Kilonovae

The coalescence of double neutron star (NS–NS) and black hole (BH)–NS binaries are prime sources of gravitational waves (GW) for Advanced LIGO/Virgo and future ground-based detectors. Neutron-rich matter released from such events undergoes rapid neutron capture (r-process) nucleosynthesis as it decompresses into space, enriching our universe with rare heavy elements like gold and platinum. Radioactive decay of these unstable nuclei powers a rapidly evolving, approximately isotropic thermal transient known as a “kilonova”, which probes the physical conditions during the merger and its aftermath. Here I review the history and physics of kilonovae, leading to the current paradigm of day-timescale emission at optical wavelengths from lanthanide-free components of the ejecta, followed by week-long emission with a spectral peak in the near-infrared (NIR). These theoretical predictions, as compiled in the original version of this review, were largely confirmed by the transient optical/NIR counterpart discovered to the first NS–NS merger, GW170817, discovered by LIGO/Virgo. Using a simple light curve model to illustrate the essential physical processes and their application to GW170817, I then introduce important variations about the standard picture which may be observable in future mergers. These include \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document}∼hour-long UV precursor emission, powered by the decay of free neutrons in the outermost ejecta layers or shock-heating of the ejecta by a delayed ultra-relativistic outflow; and enhancement of the luminosity from a long-lived central engine, such as an accreting BH or millisecond magnetar. Joint GW and kilonova observations of GW170817 and future events provide a new avenue to constrain the astrophysical origin of the r-process elements and the equation of state of dense nuclear matter.

sample 1 of BH-BH merger events (Abbott et al. 2019b) is already being used to place constraints on the formation channels of compact binary systems (e.g., Abbott et al. 2016), as well as fundamental tests of general relativity in the previously inaccessible strong-field regime (e.g., Miller 2016; Abbott et al. 2019c). We are fortunate witnesses to the birth of a new field of research: Gravitational-Wave Astronomy.
On August 17, 2017, near the end of their second observing run, Advanced LIGO/Virgo detected its first merger of a double NS binary (Abbott et al. 2017b). This event, like other GW detections, was dubbed GW170817 based on its date of discovery. The individual masses of the binary components measured by the GW signal, M 1 , M 2 ≈ 1.16-1.60 M (under the assumption of low NS spin) and the precisely-measured chirp mass are fully consistent with being drawn from the known population of Galactic binary NSs (Abbott et al. 2019a), particularly those with GW merger times less than the age of the universe (Zhao and Lattimer 2018). The lack of evidence for tidal interaction between the merging objects in the inspiral GW signal allowed for stringent upper limits to be placed on the tidal deformability and radii of NSs (Abbott et al. 2018;De et al. 2018), properties closely related to the pressure of neutron-rich matter above nuclear saturation density (see Horowitz et al. 2019 for a recent review). Beyond information encoded in the GW strain data, the discovery of electromagnetic (EM) emission accompanying the GW chirp has the potential to reveal a much richer picture of these events (Bloom et al. 2009). By identifying the host galaxies of the merging systems, and their precise locations within or around their hosts, we obtain valuable information on the binary formation channels, age of the stellar population, evidence for dynamical formation channels in dense stellar systems, or displacement due to supernova (SN) birth kicks, in a manner similar to techniques long applied to gamma-ray bursts (GRBs) and supernovae (e.g., Fruchter et al. 2006;Fong and Berger 2013). From the host galaxy redshifts, we obtain independent distance estimates to the sources, thus reducing degeneracies in the GW parameter estimation, especially of the binary inclination with respect to the line of sight (e.g., Cantiello et al. 2018;Chen et al. 2019). Redshift measurements also enable the use of GW events as standard rulers to measure the Hubble constant, or more generally, probe the cosmic expansion history (Schutz 1986;Holz and Hughes 2005;Nissanke et al. 2013). Remarkably, all of these these opportunities, and many others to be discussed later in this review, became reality with GW170817.
Except in rare circumstances, the mergers of stellar-mass BH-BH binaries are not expected to produce luminous EM emission due to the absence of significant matter surrounding these systems at the time of coalescence. Fruitful synthesis of the GW and EM skies will therefore most likely first be achieved from NS-NS and BH-NS mergers. Given the discovery of a single NS-NS merger in the O1 and O2 observing runs, LIGO/Virgo infer a volumetric rate of 110-3840 Gpc −3 yr −1 (Abbott et al. 2019b), corresponding to an expected NS-NS rate of ≈ 6-120 yr −1 once Advanced LIGO/Virgo reach design sensitivity by the early 2020s (Abbott et al. 2017d). The O1/O2 upper limit on the NS-BH merger rate is 600 Gpc −3 yr −1 (Abbott et al. 2019b). This range is broadly consistent with theoretical expectations on the rates (e.g., population synthesis models of field binaries; e.g., Dominik et al. 2015) as well as those derived empirically from the known population of Galactic double NS systems (Phinney 1991;Kalogera et al. 2004;Kim et al. 2015; see Abadie et al. 2010 for a review of rate predictions).
Among the greatest challenge to the joint EM/GW endeavor are the large uncertainties in the measured sky positions of the GW sources, which are primarily determined by triangulating the GW arrival times with an array of interferometers. When a detection is made by just the two North American LIGO facilities, sky error regions are very large (e.g., ≈ 850 deg 2 for the BH-BH merger GW150914, though later improved to ≈ 250 deg 2 ; Abbott et al. 2016a, b). However, with the addition of the Virgo detector in Italy, and eventually KAGRA in Japan (Somiya 2012) and LIGO-India, these can be reduced to more manageable values of ∼ 10 -100 deg 2 or less (Fairhurst 2011;Nissanke et al. 2013;Rodriguez et al. 2014). Indeed, information from Virgo proved crucial in reducing the sky error region of GW170817 to 30 deg 2 (Abbott et al. 2017c), greatly facilitating the discovery of its optical counterpart. Nevertheless, even in the best of cases, the GW-measured sky areas still greatly exceed that covered in a single pointing by most radio, optical, and X-ray telescopes, especially those with the required sensitivity to detect the potentially dim EM counterparts of NS-NS and BH-NS mergers (Metzger and Berger 2012). Figure 1 summarizes the EM counterparts of NS-NS and BH-NS mergers as a function of the observer viewing angle relative to the binary axis. Multiple lines of evidence, both observational (e.g., Fong et al. 2014) and theoretical 2 (Eichler et al. 1989;Narayan et al. 1992), support an association between NS-NS/BH-NS mergers and the "short duration" class of GRBs. The latter are those bursts with durations in the gamma-ray band less than about 2 s (Nakar 2007;Berger 2014), in contrast to the longer lasting bursts of duration 2 s which are instead associated with the core collapse of very massive stars. For a typical LIGO/Virgo source distance of 200 Mpc, any gammaray transient with properties matching those of the well-characterized cosmological population of GRBs would easily be detected by the Fermi, Swift or Integral satellites within their fields of view, or even with the less sensitive but all-sky Interplanetary Network of gamma-ray telescopes (Hurley 2013).

Gamma-ray bursts
The tightly collimated, relativistic outflows responsible for short GRBs are commonly believed to be powered by the accretion of a massive remnant disk onto the Fig. 1 Summary of the electromagnetic counterparts of NS-NS and BH-NS mergers and their dependence on the viewing angle with respect to the axis of the GRB jet. The kilonova, in contrast to the GRB and its afterglow, is relatively isotropic and thus represents the most promising counterpart for the majority of GW-detected mergers. Image reproduced with permission from Metzger and Berger (2012), copyright by AAS compact BH or NS remnant following the merger (e.g., Narayan et al. 1992). This is expected to occur within seconds of the merger, making their temporal association with the termination of the GW chirp unambiguous (the gamma-ray sky is otherwise quiet). Once a GRB is detected, its associated afterglow can in many cases be identified by promptly slewing a sensitive X-ray telescope to the location of the burst. This exercise is now routine with Swift, but may become less so in the future without a suitable replacement mission. Although gamma-ray detectors themselves typically provide poor sky localization, the higher angular resolution of the X-ray telescope allows for the discovery of the optical or radio afterglow; this in turn provides an even more precise position, which can help to identify the host galaxy.
A prompt burst of gamma-ray emission was detected from GW170817 by the Fermi and Integral satellites with a delay of ≈ 1.7 s from the end of the inspiral (Abbott et al. 2017d;Goldstein et al. 2017;Savchenko et al. 2017). However, rapid localization of the event was not possible, for two reasons: (1) the merger was outside the field-of-view of the Swift BAT gamma-ray detector and therefore a relatively precise sky position was not immediately available; (2) even if rapidly slewing of the X-ray telescope had been made, the X-ray afterglow may not have been detectable at such early times. Deep upper limits on the X-ray luminosity of GW170817 at t = 2.3 days ) reveal a much dimmer event than expected for a cosmological GRB placed at the same distance at a similar epoch. As we discuss below, the delayed rise and low luminosity of the synchrotron afterglow were the result of our viewing angle being far outside the core of the ultra-relativistic GRB jet, unlike the nearly on-axis orientation from which cosmological GRBs are typically viewed (e.g., Ryan et al. 2015).
Although short GRBs are probably the cleanest EM counterparts, their measured rate within the Advanced LIGO detection volume, based on observations prior to GW170817, was expected to low, probably less than once per year to decade (Metzger and Berger 2012). The measured volumetric rate of short GRBs in the local universe of R SGRB ∼ 5 Gpc −3 yr −1 (Wanderman and Piran 2015) can be reconciled with the much higher NS-NS merger rate R BNS ∼ 10 3 Gpc −3 yr −1 (Abbott et al. 2019b) if the gamma-ray emission is beamed into a narrow solid angle 4π by the bulk relativistic motion of the GRB jet (Fong et al. 2015;Troja et al. 2016). Given a typical GRB jet opening angle of θ jet ≈ 0.1 radians, the resulting beaming fraction of f −1 b = θ 2 jet /2 ∼ 1/200 ∼ R SGRB /R BNS , is consistent with most or all short GRBs arising from NS-NS mergers (though uncertainties remain large enough that a contribution from other channels, such as BH-NS mergers, cannot yet be excluded).
While the discovery of gamma-rays from the first GW-detected NS-NS merger came as a surprise to most, its properties were also highly unusual. The isotropic luminosity of the burst was ∼ 10 3 times smaller than the known population of cosmological short GRBs on which prior rate estimates had been based. A similar burst would have been challenging to detect at even twice the distance of GW170817. Several different theoretical models were proposed to explain the origin of gamma-ray signal from GW170817 (e.g., Granot et al. 2017;Fraija et al. 2019;Gottlieb et al. 2018;Beloborodov et al. 2018), but in all cases its low luminosity and unusual spectral properties are related to our large viewing angle θ obs ≈ 0.4 relative to the binary angular momentum Abbott et al. 2019a) being 4 times larger than the opening angle of the jet core, θ jet 0.1 (e.g., Fong et al. 2017;Mooley et al. 2018;Beniamini et al. 2019). Given that the majority of future GW-detections will occur at greater distances and from even larger θ obs (for which the prompt jet emission is likely to be less luminous) than for GW170817, it remains true that only a fraction of NS-NS mergers detected are likely to be accompanied by detectable gamma-rays (e.g., Mandhai et al. 2018;Howell et al. 2019). Nevertheless, given the unique information encoded in the prompt gamma-ray emission when available (e.g., on the properties of the earliest ejecta and the timing of jet formation relative to binary coalescence), every effort should be made to guarantee the presence of a wide-field gamma-ray telescopes in space throughout the next decades.
For the majority of GW-detected mergers viewed at θ obs θ jet , the most luminous GRB emission will be beamed away from our line of sight by the bulk relativistic motion of the emitting jet material. However, as the relativistic jet slows down by colliding with and shocking the interstellar medium, even off-axis viewers enter the causal emission region of the synchrotron afterglow (e.g., Totani and Panaitescu 2002). Such a delayed off-axis non-thermal afterglow emission was observed from GW170817 at X-ray (e.g., Troja et al. 2017;Margutti et al. 2017;Haggard et al. 2017), radio (e.g., Hallinan et al. 2017;Alexander et al. 2017), and optical frequencies (after the kilonova faded; e.g., Lyman et al. 2018). The afterglow light curve, which showed a gradual rise to a peak at t ∼ 200 days followed by extremely rapid fading, reveals details on the angular structure of the jet (e.g., Perna et al. 2003;Lamb and Kobayashi 2017;Lazzati et al. 2018;Xie and MacFadyen 2019). The latter structured may have been imprinted by the relativistic GRB jet as it pierced through the merger ejecta (Nakar and Piran 2017;Lazzati et al. 2017). We return later to an additional possible signature of the ejecta being shock-heated by the jet on the early-time kilonova emission.

Kilonovae
In addition to the beamed GRB and its afterglow, the merger of NS-NS and BH-NS binaries are also accompanied by a more isotropic counterpart, commonly known as a 'kilonova ' (or, less commonly, 'macronova'). Kilonovae are thermal supernova-like transients lasting days to weeks, which are powered by the radioactive decay of heavy neutron-rich elements synthesized in the expanding merger ejecta (Li and Paczyński 1998). They provide both a robust EM counterpart to the GW chirp, which is expected to accompany a fraction of BH-NS mergers and essentially all NS-NS mergers, as well as a direct probe of the unknown astrophysical origin of the heaviest elements (Metzger et al. 2010b).
This article provides a pedagogical review of kilonovae, including a brief historical background and recent developments in this rapidly evolving field (Sect. 2). Section 3 describes the basic physical ingredients, including the key input from numerical relativity simulations of the merger and its aftermath. For pedagogical reasons, the discussion is organized around a simple toy model for the kilonova light curve (Sect. 4), which synthesizes most of the key ingredients in a common and easy-to-interpret framework. My goal is to make the basic results accessible to anyone with the ability to solve a set of coupled ordinary differential equations.
I begin by introducing the simplest model of lanthanide-rich ejecta heated by radioactivity, which produces a week-long near-infrared (NIR) transient ('red kilonova'; Sect. 4.1.1) which is proceeded in at least some cases by ∼ day-long UV/optical-wavelength emission ('blue kilonova') arising from lanthanide-free components of the ejecta (Sect. 4.1.2). Section 5 describes observations of the thermal UVOIR kilonova emission observed following GW170817 and its theoretical interpretation within this largely pre-existing theoretical framework. I also summarize the lessons GW170817 has provided, about the origin of r -process elements, the equation of state of neutron stars, and the final fate of the merger remnant.
Section 6 explores several variations on this canonical picture, some of which were not possible to test in the case of GW170817 and some which are ruled out in that event but could be relevant to future mergers, e.g., with different ingoing binary parameters. These include early (∼ hours-long) 'precursor' emission at ultra-violet wavelengths (UV), which is powered either by the decay of free neutrons in the outermost layers of the ejecta (Sect. 6.1.1) or prompt shock heating of the ejecta by a relativistic outflow such as the GRB jet (Sect. 6.1.2). In Sect. 6.2 we consider the impact on the kilonova signal of energy input from a long-lived accreting BH or magnetar central engine. In Sect. 7, I assess the prospects for discovering kilonovae following short GRBs and for future GW-triggers of NS-NS/BH-NS mergers. I use this opportunity to make predictions for the diversity of kilonova signals with GW-measured properties of the binary, which will become testable once EM observations are routinely conducted in Freiburghaus et al.: NS-NS dynamical ejecta ⇒ r-process abundances 2005 Kulkarni: kilonova powered by free neutron-decay ("macronova"), central engine coincidence with a large sample of GW-detected merger events. I conclude with some personal thoughts in Sect. 8. Although I have attempted to make this review self-contained, the material covered is necessarily limited in scope and reflects my own opinions and biases. I refer the reader to a number of other excellent recent reviews, which cover some of the topics discussed briefly here in greater detail: Nakar (2007), Faber and Rasio (2012), Berger (2014), Rosswog (2015), Fan and Hendry (2015), Baiotti and Rezzolla (2017), Baiotti (2019), including other short reviews dedicated exclusively to kilonovae (Tanaka 2016;Yu 2019). I encourage the reader to consult Fernández and Metzger (2016) for a review of the broader range of EM counterparts of NS-NS/BH-NS mergers. A few complementary reviews have appeared since GW170817 overviewing the interpretation of this event (e.g., Miller 2017; Bloom and Sigurdsson 2017;Metzger 2017b;Siegel 2019) or its constraints on the nuclear EOS (e.g., Raithel 2019). I also encourage the reader to consult the initial version of this review, written the year prior to the discovery of GW170817 (Metzger 2017a). Table 1 summarizes key events in the historical development of kilonovae. Burbidge et al. (1957) and Cameron (1957) realized that approximately half of the elements heavier than iron are synthesized via the capture of neutrons onto lighter seed nuclei like iron) in a dense neutron-rich environment in which the timescale for neutron capture is shorter than the β-decay timescale. This 'rapid neutron-capture process', or r -process, occurs along a nuclear path which resides far on the neutron-rich side of the valley of stable isotopes. Despite these works occurring over 70 years ago, the astrophysical environments giving rise to the r -process remains an enduring mystery, among the greatest in nuclear astrophysics (e.g., Qian and Wasserburg 2007;Arnould et al. 2007;Thielemann et al. 2011;Cowan et al. 2019, for contemporary reviews).

NS mergers as sources of the r-process
Among the most critical quantities which characterize the viability of a potential r -process event is the electron fraction of the ejecta, where n p and n n are the densities of protons and neutrons, respectively. Ordinary stellar material usually has more protons than neutrons (Y e ≥ 0.5), while matter with a neutron excess (Y e < 0.5) is typically required for the r -process. Core collapse supernovae have long been considered promising r -process sources. This is in part due to their short delays following star formation, which allows even the earliest generations of metal-poor stars in our Galaxy to be polluted with r -process elements (e.g., Mathews et al. 1992;Sneden et al. 2008). Throughout the 1990s, the high entropy 3 neutrino-heated winds from proto-neutron stars (Duncan et al. 1986;Qian and Woosley 1996), which emerge seconds after a successful explosion, were considered the most likely r -process site 4 within the core collapse environment (Woosley et al. 1994;Takahashi et al. 1994). However, more detailed calculations of the wind properties (Thompson et al. 2001;Arcones et al. 2007;Fischer et al. 2010;Hüdepohl et al. 2010;Roberts et al. 2010;Martínez-Pinedo et al. 2012;Roberts et al. 2012) later showed that the requisite combination of neutron-rich conditions (Y e 0.5) and high entropy were unlikely to obtain. Possible exceptions include the rare case of a very massive proto-NS (Cardall and Fuller 1997), or in the presence of non-standard physics such as an eV-mass sterile neutrino (Tamborra et al. 2012;Wu et al. 2014).
Another exception to this canonical picture may occur if the NS is formed rotating rapidly and is endowed with an ultra-strong ordered magnetic field B 10 14 -10 15 G, similar to those which characterize Galactic magnetars. Magneto-centrifugal acceleration within the wind of such a "millisecond magnetar" to relativistic velocities can act to lower its electron fraction or reduce the number of seed nuclei formed through rapid expansion ). This could occur during the early supernova explosion phase Nishimura et al. 2015) as well as during the subsequent cooling phase of the proto-NS over several seconds (Thompson 2003;Metzger et al. 2007;Vlasov et al. 2014). Despite the promise of such models, numerical simulations of MHD supernovae are still in a preliminary state, especially when it comes to the accurate neutrino transport need to determine the ejecta Y e and the high resolution three-dimensional grid needed to capturing the growth of non-axisymmetric (magnetic 3 A high entropy (low density) results in an α-rich freeze-out of the 3-body and effective 4-body reactions responsible for forming seed nuclei in the wind, similar to big bang nucleosynthesis. The resulting higher ratio of neutrons to seed nuclei (because the protons are trapped in α particles) then allows the r -process to proceed to heavier elements. 4 Another r -process mechanism in the core collapse environment results from ν−induced spallation in the He shell (Banerjee et al. 2011). This channel is limited to very low metallicity Z 10 −3 and thus cannot represent the dominant r -process source over the age of the galaxy (though it could be important for the first generations of stars). kink or sausage mode) instabilities. The latter can disrupt and slow the expansion rate of jet-like structures (Mösta et al. 2014), rendering the creation of the heaviest (third abundance-peak) r -process elements challenging to obtain (Halevi and Mösta 2018).
The observed rate of hyper-energetic supernovae (the only bona fide MHD-powered explosions largely agreed upon to exist in nature) is only ∼ 1/1000 of the total core collapse supernova rate (e.g., Podsiadlowski et al. 2004). Therefore, a higher r -process yield per event 10 −2 M is required to explain a significant fraction of the Galactic abundances through this channel. However, in scenarios where the r -process takes place in a prompt jet during the supernova explosion, it is inevitable that the r -process material will mix into the outer layers of the supernova ejecta along with the shocksynthesized 56 Ni, the latter being responsible for powering the supernova's optical luminosity. As we shall discuss later in the context of kilonovae, such a large abundance of lanthanide elements mixed into the outer ejecta layers would substantially redden the observed colors of the supernova light in a way incompatible with observed hyperenergetic (MHD) supernovae (Siegel et al. 2019). 5 Contemporaneously with the discovery of the first binary pulsar (Hulse and Taylor 1975), Schramm (1974, 1976) proposed that the merger of compact star binaries-in particular the collision of BH-NS systems-could give rise to the r -process by the decompression of highly neutron-rich ejecta (Meyer 1989). Symbalisty and Schramm (1982) were the first to suggest NS-NS mergers as the site of the r -process. Blinnikov et al. (1984) and Paczyński (1986) first suggested a connection between NS-NS mergers and GRBs. Eichler et al. (1989) presented a more detailed model for how this environment could give rise to a GRB, albeit one which differs significantly from the current view. Davies et al. (1994) performed the first numerical simulations of mass ejection from merging neutron stars, finding that ∼ 2% of the binary mass was unbound during the process. Freiburghaus et al. (1999) presented the first explicit calculations showing that the ejecta properties extracted from a hydrodynamical simulation of a NS-NS merger  indeed produces abundance patterns in basic accord with the solar system r -process.
The neutrino-driven wind following a supernova explosion accelerates matter from the proto-NS surface relatively gradually, in which case neutrino absorption reactions on nucleons (particularly ν e +n → p +e − ) have time to appreciably raise the electron fraction of the wind from its initial low value near the NS surface. By contrast, in NS-NS/BH-NS mergers the different geometry and more dynamical nature of the system allows at least a fraction of the unbound ejecta (tidal tails and disk winds) to avoid strong neutrino irradiation, maintaining a much lower value of Y e 0.2 (Sect. 3.1).
When averaged over the age of the Galaxy, the required production rate of heavy r -process nuclei of mass number A > 140 is ∼ 2 × 10 −7 M yr −1 (Qian 2000). Given a rate R NS−NS of detection of NS-NS mergers by Advanced LIGO/Virgo at design sensitivity (defined here as a horizon distance of 200 Mpc for NS-NS mergers), the required r -process mass yield per merger event to explain the entire Galactic abundances is very approximately given by (e.g., Metzger et al. 2009;Vangioni et al. 2016) As described in Sect. 3.1, numerical simulations of NS-NS/BH-NS mergers find a range of total ejecta masses of M r ∼ 10 −3 −10 −1 M , while M r ≈ 0.03−0.06 M was inferred from the kilonova of GW170817 (Sect. 5). Although large uncertainties remain, it is safe to conclude that NS mergers are likely major, if not the dominant, sources of the r -process in the universe. Several additional lines of evidence support 'high yield' r -process events consistent with NS-NS/BH-NS mergers being common in our Galaxy, both now and in its early history. These include the detection of 244 Pu on the ocean floor at abundances roughly two orders lower than that expected if the currently active r -process source were frequent, low-yield events like ordinary core collapse supernovae (Wallner et al. 2015;Hotokezaka et al. 2015). Similar arguments show that actinide abundances in the primordial solar system require a rare source for the heaviest r -process elements (Bartos and Marka 2019;Côté et al. 2019b). A fraction of the stars in the ultra-faint dwarf galaxy Reticulum II are highly enriched in r -process elements, indicating that this galaxy was polluted early in its history by a single rare r -process event (Ji et al. 2016). Similarly, the large spread seen in the r -process abundances in many Globular Clusters (e.g., Roederer 2011) may also indicate a rare source acting at low metallicity. Reasonable variations in the ejecta properties of NS mergers could in principle explain the observed variability in the actinide abundances of metal-poor stars (Holmbeck et al. 2019).
Nevertheless, NS mergers are challenged by some observations, which may point to alternative r -process sites. Given the low escape speeds of dwarf galaxies of ∼ 10 km s −1 , even moderate velocity kicks to binaries from the process of NS formation would remove the binaries from the galaxy prior to merger (and thus preventing the merger ejecta from polluting the next generation of stars). Although a sizable fraction of the Galactic NS-NS binaries have low proper motions and are indeed inferred to have experienced very low supernova kicks (Beniamini et al. 2016), even relatively modest spatial offsets of the merger events from the core of the galaxies make it challenging to retain enough r -process enhanced gas (Bonetti et al. 2019). Another challenge to NS mergers are the short delay times 10-100 Myr between star formation and merger which are required to explain stellar populations in the low-metallicity Galactic halo (Safarzadeh et al. 2019) and Globular Clusters (Zevin et al. 2019). Depending on the efficiency of compositional mixing between the merger ejecta and the ISM of the Galaxy, realistic delay time distributions for NS-NS/NS-BH mergers within a consistent picture of structure formation via hierarchical growth (Kelley et al. 2010) were argued to produce chemical evolution histories consistent with observations of the abundances of r -process elements in metal-poor halo stars as a function of their iron abundance (Shen et al. 2015;Ramirez-Ruiz et al. 2015;van de Voort et al. 2015). However, Safarzadeh et al. 2019 comes to a different conclusion, while van de Voort et al. (2019) found that r -process production by rare supernovae better fit the abundances of metal-poor stars than NS mergers. Due to the fact that the delay time distribution of NS mergers at late times after star formation is expected to similar to that of Type Ia supernova (which generate most of the iron in the Galaxy), NS mergers are also challenged to explain the observed decrease of [Eu/Fe] with increasing [Fe] at later times in the chemical evolution history of the Galaxy Côté et al. 2019a).
Together, these deficiencies may hint at the existence of an additional high-yield r -process channel beyond NS-NS mergers which can operate at low metallicities with short delay times. The collapse of massive, rotating stars ("collapsars"), the central engines of long-duration gamma-ray bursts (which are directly observed to occur in dwarf low-metallicity galaxies), are among the most promising contenders (Pruet et al. 2004;Fryer et al. 2006;Siegel et al. 2019). As we shall discuss, the physical conditions of hyper-accreting disks and their outflows following NS mergers, as probed by their kilonova emission, offers an indirect probe of the broadly similar physical conditions which characterize the outflows generated in collapsars. Evidence for r -process in the outflows of NS merger accretion flows would indirectly support an r -process occurring in collapsars as well. Li and Paczyński (1998) first argued that the radioactive ejecta from a NS-NS or BH-NS merger provides a source for powering thermal transient emission, in analogy with supernovae. They developed a toy model for the light curve, similar to that we describe in Sect. 4. Given the low mass and high velocity of the ejecta from a NS-NS/BH-NS merger, they concluded that the ejecta will become transparent to its own radiation quickly, producing emission which peaks on a timescale of about one day, much faster than for normal supernovae (which instead peak on a timescale of weeks or longer).

A brief history of kilonovae
Lacking a model for the nucleosynthesis (the word "r -process" does not appear in their work), Li and Paczyński (1998) parametrized the radioactive heating rate of the ejecta at time t after the merger according to the following prescription, where M is the ejecta mass and f was a free parameter. The ∝ 1/t time dependence was motivated by the total heating rate which results from the sum of the radioactive decay heating rateQ i ∝ exp(−t/τ i ) of a large number of isotopes i, under the assumption that their half-lives τ i are distributed equally per logarithmic time (at any time t, the heating rate is dominated by isotopes with half-lives τ i ∼ t). Contemporary models, which process the thermodynamic history of the expanding ejecta based on numerical simulations of the merger through a detailed nuclear reaction network, show that the heating rate at late times actually approaches a steeper power law decay ∝ t −α , with α ≈ 1.1-1.4 (Metzger et al. 2010b;Roberts et al. 2011;Korobkin et al. 2012), similar to that found for the decay rate of terrestrial radioactive waste (Way and Wigner 1948). Metzger et al. (2010b) and Hotokezaka et al. (2017) describe how this power-law decay can be understood from the basic physics of β-decay and the properties of nuclei on the neutron-rich valley of stability. Li and Paczyński (1998) left the normalization of the heating rate f , to which the peak luminosity of the kilonova is linearly proportional, as a free parameter, considering a range of models with different values of f = 10 −5 -10 −3 . More recent calculations, described below, show that such high heating rates are extremely optimistic, leading to predicted peak luminosities 10 43 -10 44 erg s −1 (Li and Paczyński 1998, their Fig. 2) which exceed even those of supernovae. These overpredictions leaked to other works throughout the next decade; for instance, Rosswog (2005) predicted that BH-NS mergers are accompanied by transients of luminosity 10 44 erg s −1 , which would rival the most luminous transients ever discovered. This unclear theoretical situation led to observational searches for kilonovae following short GRBs which were inconclusive since they were forced to parametrized their results (usually non-detections) in terms of the allowed range of f (Bloom et al. 2006;Kocevski et al. 2010) instead of in terms of more meaningful constraints on the ejecta properties such as its mass. Metzger et al. (2010b) were the first to determine the true luminosity scale of the radioactively-powered transients of NS mergers by calculating light curve models using radioactive heating rates derived self-consistently from a nuclear reaction network calculation of the r -process, based on the dynamical ejecta trajectories of Freiburghaus et al. (1999). 6 Based on their derived peak luminosities being approximately one thousand times brighter than a nova, Metzger et al. (2010b) introduced the term 'kilonova' to describe the EM counterparts of NS mergers powered by the decay of r -process nuclei. They showed that the radioactive heating rate was relatively insensitive to the precise electron fraction of the ejecta and to the assumed nuclear mass model, and they were the first to consider how efficiently the decay products thermalize their energy in the ejecta. Their work highlighted the critical four-way connection, now taken for granted, between kilonovae, short GRBs, GWs from NS-NS/BH-NS mergers, and the astrophysical origin of the r -process.
Prior to Metzger et al. (2010b), it was commonly believed that kilonovae were in fact brighter, or much brighter, than supernovae (Li and Paczyński 1998;Rosswog 2005). One exception is Kulkarni (2005), who assumed that the radioactive power was supplied by the decay of 56 Ni or free neutrons. However, 56 Ni cannot be produced in the neutron-rich ejecta of a NS merger, while all initially free neutrons are captured into seed nuclei during the r -process (except perhaps in the very outermost, fastest expanding layers of the ejecta; see Sect. 6.1.1). Kulkarni introduced the term "macronovae" for such Nickel/neutron-powered events. Despite its inauspicious physical motivation and limited use in the literature until well after the term kilonova was already in use, many authors continue to use the macronova terminology, in part because this name is not tied to a particular luminosity scale (which may change as our physical models evolve). 6 As a student entering this field in the mid 2000s, it was clear to me that the optical transients proposed by Li and Paczyński (1998) were not connected in most people's mind with the r -process. Rosswog (2005) in principle had all the information needed to calculate the radioactive heating rate of the ejecta based on the earlier Freiburghaus et al. (1999) calculations, and thus to determine the true luminosity scale of these merger transients well before Metzger et al. (2010b). I make this point not to cast blame, but simply to point out that the concept, now taken for granted, that the radioactive heating rate was something that could actually be calculated with any precision, came as a revelation, at least to a student of the available literature.
Once the radioactive heating rate was determined, attention turned to the yet thornier issue of the ejecta opacity. The latter is crucial since it determines at what time and wavelength the ejecta becomes transparent and the light curve peaks. Given the general lack of experimental data or theoretical models for the opacity of heavy r -process elements, especially in the first and second ionization states of greatest relevance, Metzger et al. (2010b), Roberts et al. (2011) adopted grey opacities appropriate to the Fe-rich ejecta in Type Ia supernovae. However, Kasen et al. (2013) showed that the opacity of r -process elements can be significantly higher than that of Fe, due to the high density of line transitions associated with the complex atomic structures of some lanthanide and actinide elements (Sect. 3.2). This finding was subsequently confirmed by Tanaka and Hotokezaka (2013). As compared to the earlier predictions (Metzger et al. 2010b), these higher opacities push the bolometric light curve to peak later in time (∼ 1 week instead of a ∼ 1 day timescale), and at a lower luminosity . More importantly, the enormous optical wavelength opacity caused by line blanketing moved the spectral peak from optical/UV frequencies to the nearinfrared (NIR). Later that year, Tanvir et al. (2013) and Berger et al. (2013) presented evidence for excess infrared emission following the short GRB 130603B on a timescale of about one week using the Hubble Space Telescope.
However, not all of the merger ejecta necessarily will contain lanthanide elements with such a high optical opacity (e.g., Metzger et al. 2008a). While ejecta with a relatively high electron fraction 0.25 Y e 0.4 has enough neutrons to synthesize radioactive r -process nuclei, the ratio of neutrons to lighter seed nuclei is insufficient to reach the relatively heavy lanthanide elements of atomic mass number A 140 which (if not blocked by high-opacity low-Y e material further out) produces emission with more rapid evolution and bluer colors, similar to those predicted by the original models (Metzger et al. 2010b;Roberts et al. 2011). Metzger and Fernández (2014) dubbed the emission from high-Y e , lanthanide-poor ejecta a "blue" kilonova, in contrast to "red" kilonova emission originating from low-Y e , lanthanide-rich portions of the ejecta . They further argued that both "blue" and "red" kilonova emission, arising from different components of the merger ejecta, could be seen in the same merger event, at least for some observing geometries. As will be discussed in §5, such hybrid "blue" + "red" kilonova models came to play an important role in the interpretation of GW170817. Figure 2 is a timeline of theoretical predictions for the peak luminosities, timescales, and spectral peak of the kilonova emission.

Basic ingredients
The physics of kilonovae can be understood from basic considerations. Consider the merger ejecta of total mass M, which is expanding at a constant mean velocity v, such that its mean radius is R ≈ vt after a time t following the merger. Perhaps surprisingly, it is not unreasonable to assume spherical symmetry to first order because the ejecta will have a chance to expand laterally over the many orders of magnitude in scale from the merging binary (R 0 ∼ 10 6 cm) to the much larger radius (R peak ∼ 10 15 cm) at which the kilonova emission peaks (Roberts et al. 2011;Grossman et al. 2014;Rosswog et al. 2014).

Fig. 2
Schematic timeline of the development kilonova models in the space of peak luminosity and peak timescale. The wavelength of the predicted spectral peak are indicated by color as marked in the figure. Shown for comparison are the approximate properties of the "red" and "blue" kilonova emission components observed following GW170817 (e.g., Cowperthwaite et al. 2017;Villar et al. 2017) The ejecta is extremely hot immediately after being ejected from the viscinity of the merger (Sect. 3.1). This thermal energy cannot, however, initially escape as radiation because of its high optical depth at early times, and the correspondingly long photon diffusion timescale through the ejecta, where ρ = 3M/(4π R 3 ) is the mean density and κ is the opacity (cross section per unit mass). As the ejecta expands, the diffusion time decreases with time t diff ∝ t −1 , until eventually radiation can escape on the expansion timescale, as occurs once t diff = t (Arnett 1982). This condition determines the characteristic timescale at which the light curve peaks, where the constant β ≈ 3 depends on the precise density profile of the ejecta (see Sect. 4). For values of the opacity κ ∼ 0.5-30 cm 2 g −1 which characterize the range from lanthanide-free and lanthanide-rich matter ; Table 4), respectively, Eq. (7) predicts characteristic durations ∼ 1 day -1 week.
The temperature of matter freshly ejected at the radius of the merger R 0 10 6 cm generally exceed 10 9 −10 10 K. However, absent a source of persistent heating, this Left: sources of radioactive heating include the decay of ∼ 10 −2 M of r-process nuclei, as first modeled in a parametrized way by Li and Paczyński (1998) (Eq. 4, grey band) and then by Metzger et al. (2010b) using a full reaction network, plotted here using the analytic fit of Korobkin et al. (2012) (Eq. 22, black line) and including the thermalization efficiency of Barnes et al. (2016) (Eq. 25). The outermost layers of the ejecta may contain ∼ 10 −4 M free neutrons (red line), which due to their comparatively long half-life can enhance the kilonova emission during the first few hours if present in the outermost layers of the ejecta due to premature freeze-out of the r -process (Sect. 6.1.1). Right: heating sources from a central engine. These include fall-back accretion (blue lines), shown separately for NS-NS (solid line) and BH-NS (dashed line) mergers, based on results by Rosswog (2007) for an assumed jet efficiency j = 0.1 (Eq. 34). Also shown is the rotational energy input from the magnetic dipole spin-down of a stable magnetar remnant with an initial spin period of P = 0.7 ms and dipole field strengths of B = 10 15 G (brown lines) and 10 16 G (orange lines). Dashed lines show the total spin-down luminosity L sd (Eq. 36), while solid lines show the effective luminosity available to power optical/X-ray emission once accounting for suppression of the efficiency of thermalization due to the high scattering opacity of e ± pairs in the nebula (Eq. 39; Metzger and Piro, 2014). The isotropic luminosity of the temporally-extended X-ray emission observed following the short GRB 080503 is shown with a green line (for an assumed source redshift z = 0.3; Perley et al. 2009) matter will cool through adiabatic expansion, losing all but a fraction ∼ R 0 /R peak ∼ 10 −9 of its initial thermal energy before reaching the radius R peak = vt peak at which the ejecta becomes transparent (Eq. 7). Such adiabatic losses would leave the ejecta so cold as to be effectively invisible at large distances.
In a realistic situation, the ejecta will be continuously heated, by a combination of sources, at a total rateQ(t) (Fig. 3). At a minimum, this heating includes contributions from radioactivity due to r -process nuclei and, possibly at early times, free neutrons. More speculatively, the ejecta can also be heated from within by a central engine, such as the emergence of the GRB jet or over longer timescales by the rotational energy of a magnetar remnant. In most cases of relevance,Q(t) is constant or decreasing with time less steeply than ∝ t −2 . The peak luminosity of the observed emission then equals the heating rate at the peak time (t = t peak ), i.e., a result commonly known as "Arnett's Law" (Arnett 1982). Equations (7) and (8) make clear that, in order to quantify the key observables of kilonovae (peak timescale, luminosity, and effective temperature), we must understand three key ingredients: -The mass and velocity of the ejecta from NS-NS/BH-NS mergers, as comprised by several distinct components. -The opacity κ of expanding neutron-rich matter.
-The variety of sources contributing to the ejecta heatingQ(t), particularly on timescales of t peak , when the ejecta first becomes transparent.
The remainder of this section addresses the first two issues. The range of different heating sources, which can give rise to different types of kilonovae, are covered in Sect. 4.

Sources of neutron-rich ejecta
Two broad sources of ejecta characterize NS-NS and BH-NS mergers (see Fernández and Metzger 2016; Shibata and Hotokezaka 2019, for recent reviews). First, there is matter ejected on the dynamical timescale of milliseconds, either by tidal forces or due to compression-induced heating at the interface between merging bodies (Sect. 3.1.1). Debris from the merger, which is not immediately unbound or incorporated into the central compact object, can possess enough angular momentum to circularize into an accretion disk around the central remnant. A disk can also be generated by outwards transport of angular momentum and mass during the post-merger evolution of the central NS remnant prior to BH formation. Outflows from this remnant disk, taking place on longer timescales of up to seconds, provide a second important source of ejecta (Sect. 3.1.2) ( Table 2). In BH-NS mergers, significant mass ejection and disk formation occurs only if the BH has a low mass M • and is rapidly spinning; in such cases, the NS is tidally disrupted during the very final stages of the inspiral instead of being swallowed whole (giving effectively zero mass ejection). Roughly speaking, the condition for the latter is that the tidal radius of the NS, R t ∝ M 1/3 • , exceed the innermost stable circular orbit of the BH, R isco ∝ M • (see Foucart 2012;Foucart et al. 2018 for a more precise criterion for mass ejection, calibrated to GR merger simulations). For a NS of radius 12 km and mass 1.4 M , this requires a BH of mass 4(12) M for a BH Kerr spin parameter of χ BH = 0.7(0.95). For slowly-spinning BHs (as appears to characterize most of LIGO/Virgo's BH-BH systems), the BH mass range giving rise to tidal disruptionand hence a kilonova or GRB-could be very small.
In the case of a NS-NS merger, the ejecta properties depend sensitively on the fate of the massive NS remnant which is created by the coalescence event. The latter  Köppel et al. 2019), which they find is insensitive to the binary mass ratio q = M 2 /M 1 for q 0.7 (however, see Kiuchi et al. 2019).
Mergers that do not undergo prompt collapse (M tot < M crit ) typically result in the formation of rapidly-spinning NS remnant of mass ∼ M tot (after subtracting mass lost through neutrino and GW emission and in the dynamical ejecta), which is at least temporarily stable against gravitational collapse to a BH. The maximum stable mass of a NS exceeds its non-rotating value, M TOV , if the NS is rapidly spinning close to the break-up velocity Özel et al. 2010;Kaplan et al. 2014).
A massive NS remnant, which is supported exclusively by its differential rotation, is known as a hypermassive NS (HMNS). A somewhat less massive NS, which can be supported even by its solid body rotation (i.e., after differential rotation has been removed), is known as a supramassive NS (SMNS). A HMNS is unlikely to survive for more than a few tens to hundreds of milliseconds after the merger, before collapsing to a BH due to the loss of differential rotation and accretion of mass by internal hydromagnetic torques and gravitational wave radiation (Shibata and Taniguchi 2006;Duez et al. 2006;Siegel et al. 2013). In contrast, SMNS remnants must spin-down to the point of collapse through the global loss of angular momentum. The latter must take place through less efficient processes, such as magnetic dipole radiation or GW emission arising from small non-axisymmetric distortions of the NS, and hence such objects can in principle survive for much longer before collapsing. Finally, the merger of a particularly low mass binary, which leaves a remnant mass less than M TOV , will produce an indefinitely stable remnant (Metzger et al. 2008b;Giacomazzo and Perna 2013), from which a BH can never form, even once its angular momentum has been entirely removed. Such cases are likely very rare.

Fig. 4
Left: properties of the merger ejecta which affect the EM emission as a function of the binary chirp mass M c (Eq. 1), taken here as a proxy for the total binary mass M tot . Vertical dashed lines delineate the threshold masses for different merger remnants as marked, for an example EOS with M TOV = 2.1 M and radius R 1.6 = 12 km of a 1.6 M NS. The top panel shows the ejecta kinetic energy, which we take to be the sum of the initial kinetic energy of the ejecta (estimated using fits to numerical relativity simulations; Coughlin et al. 2018Coughlin et al. , 2019 and, in the case of stable or SMNSs, the rotational energy which can be extracted from the remnant before forming a BH (Margalit and Metzger 2017). The bottom panel shows the ejecta mass, both dynamical and disk wind ejecta, estimated as in Coughlin et al. (2019), where 50% of the disk mass is assumed to be ejected at v = 0.15 c (e.g., Siegel and Metzger 2017). The finite width of the lines results from a range of binary mass ratio q = 0.7-1, to which the tidal dynamical ejecta is most sensitive. The ejecta mass line is colored qualitatively according to the dominant color of the kilonova emission, which becomes redder for more massive binaries (with shorter-lived remnants) due to their more neutron-rich ejecta (Metzger and Fernández 2014). Right: distribution of BNS merger chirp masses drawn from a NS population representative of Galactic double NSs (Kiziltan et al. 2013). Dashed vertical curves separate the M c parameter space based on the possible merger outcomes in each region. The fraction of mergers expected to occur in each region (the integral over the PDF within this region) is stated above the region in red (see also Table 3) . Image reproduced with permission from Margalit and Metzger (2019), copyright by the authors Table 3  The right panel of Fig. 4 shows the chirp mass distribution of known Galactic NS-NS binaries (Kiziltan et al. 2013) compared to the allowed ranges in the binary mass thresholds separating different remnant classes (stable NS, SMNS, HMNS, prompt collapse) given current EOS constraints. The chirp mass of GW170817 is fully consistent with being drawn from the Galactic NS-NS population, while indications from the EM observations suggest that a HMNS remnant formed in this event (Sect. 5). If the extra-galactic population of merging NS-NS binaries is indeed similar to the known Galactic population, Margalit and Metzger (2019) predict that 18% − 65% of mergers would result in SMNS remnants, while only a small fraction < 3% would produce indefinitely stable NS remnants (Table 3). As we discuss in Sect. 6.2.2, additional energy input from a long-lived magnetar remnant could substantially boost the kilonova emission. The fractions of mergers leading to a prompt BH collapse, and relatively little ejecta or disk, ranges from tens of percent to extremely infrequently.

Dynamical ejecta
NS-NS mergers eject unbound matter through processes that operate on the dynamical time, and which depend primarily on the total binary mass, the mass ratio, and the EOS. Total dynamical ejecta masses typically lie in the range 10 −4 -10 −2 M for NS-NS mergers (e.g., Hotokezaka et al. 2013a;Radice et al. 2016a;Bovard et al. 2017), with velocities 0.1-0.3 c. For BH-NS mergers, the ejecta mass can be up to ∼ 0.1M with similar velocities as in the NS-NS case (Kyutoku et al. 2013Foucart et al. 2017). The ejecta mass is typically greater for eccentric binaries (East et al. 2012;Gold et al. 2012), although the dynamical interactions giving rise to eccentric mergers require high stellar densities, probably making them rare events compared to circular inspirals (Tsang 2013). Very high NS spin can also enhance the quantity of dynamical ejecta (e.g., Dietrich et al. 2017a;East et al. 2019;Most et al. 2019).
Two main ejection processes operate in NS-NS mergers. First, material at the contact interface between the merging stars is squeezed out by hydrodynamic forces and is subsequently expelled by quasi-radial pulsations of the remnant (Oechslin et al. 2007;Bauswein et al. 2013b;Hotokezaka et al. 2013a), ejecting shock-heated matter in a broad range of angular directions. The second process involves spiral arms from tidal interactions during the merger, which expand outwards in the equatorial plane due to angular momentum transport by hydrodynamic processes. The relative importance of these mechanisms depends on the EOS and the binary mass ratio q, with lower values of q 1 (asymmetric) binaries ejecting greater quantities of mass (Bauswein et al. 2013b;Lehner et al. 2016). The ejecta mass also depends on the BH formation timescale; for the prompt collapses which characterize massive binaries, mass ejection from the contact interface is suppressed due to prompt swallowing of this region. Figure 5 shows the total ejecta mass and mean velocity of the dynamical ejecta inferred from a range of NS-NS simulations compiled from the literature (Bauswein et al. 2013b;Hotokezaka et al. 2013a;Radice et al. 2018b;Sekiguchi et al. 2016;Ciolfi et al. 2017).
In BH-NS mergers, mass is ejected primarily by tidal forces that disrupt the NS, with the matter emerging primarily in the equatorial plane (Kawaguchi et al. 2015). The ejecta from BH-NS mergers also often covers only part of the azimuthal range , which may introduce a stronger viewing angle dependence on the kilonova emission than for NS-NS mergers.
Another key property of the dynamical ejecta, in addition to the mass and velocity, is its electron fraction, Y e . Simulations that do not account for weak interactions find the ejecta from NS-NS mergers to be highly neutron-rich, with an electron fraction . More recent merger calculations that include the effects of e ± captures and neutrino irradiation in full general-relativity have shown that the dynamical ejecta may have a wider electron fraction distribution (Y e ∼ 0.1−0.4) than models which neglect weak interactions Radice et al. 2016a). As a result, lighter r -process elements with 90 A 130 are synthesized in addition to third-peak elements ). These high-Y e ejecta components are distributed in a relatively spherically-symmetric geometry, while the primarily tidallyejected, lower-Y e matter is concentrated closer to the equatorial plane (Fig. 7).

Disk outflow ejecta
All NS-NS mergers, and those BH-NS mergers which end in NS tidal disruption outside the BH horizon, result in the formation of an accretion disk around the central NS or BH remnant. The disk mass is typically ∼ 0.01-0.3 M , depending on the total mass and mass ratio of the binary, the spins of the binary components, and the NS EOS (e.g., Oechslin and Janka 2006). Relatively low disk masses are expected in the case of massive binaries that undergo prompt collapse to a BH, because the process of massive disk formation is intimately related to the internal redistribution of mass and angular momentum of the remnant as it evolves from a differentially rotating to solid body state (which has no time to occur in a prompt collapse). Outflows from this disk, over a timescales of seconds or longer, represent an important source of ejecta mass which can often dominate that of the dynamical ejecta.
At early times after the disk forms, its mass accretion rate is high and the disk is a copious source of thermal neutrinos (Popham et al. 1999). During this phase, mass loss is driven from the disk surface by neutrino heating, in a manner analogous to neutrinodriven proto-NS winds in core collapse supernovae (Surman et al. 2008;Metzger et al. 2008c). Spiral density waves, which are excited in the disk by the oscillations of the central NS remnant, may also play a role in outwards angular momentum transport and mass ejection during this early phase (Nedora et al. 2019). Time dependent models of the long-term evolution of these remnant tori, which include neutrino emission and absorption, indicate that when BH formation is prompt, the amount of mass ejected through this channel is small, contributing at most a few percent of the outflow, because the neutrino luminosity decreases rapidly in time (Fernández and Metzger 2013;Just et al. 2015). However, if the central NS remnant survives for longer than ∼ 50 ms (as a HMNS or SMNS), then the larger neutrino luminosity from the NS remnant ejects a non-negligible amount of mass (∼ 10 −3 M , primarily from the NS itself instead of the disk; Dessart et al. 2009;Perego et al. 2014;Martin et al. 2015;Richers et al. 2015). As we discuss below, ejecta from the star could be substantially enhanced if the central remnant has a strong ordered magnetic field .
The disk evolves in time due to the outwards transport of angular momentum, as mediated e.g., by spiral density waves or (more generically) magnetic stresses created by MHD turbulence generated by the magneto-rotational instability. Initial time-dependent calculations of this 'viscous spreading' followed the disk evolution over several viscous times using one-zone (Metzger et al. 2008a) and one-dimensional height-integrated ) models. These works showed that, as the disk evolves and its accretion rate decreases, the disk transitions from a neutrino-cooled state to a radiatively inefficient (geometrically thick disk) state as the temperature, and hence the neutrino cooling rate, decreases over a timescale of seconds (see also Lee et al. 2009;Beloborodov 2008). Significant outflows occur once during the radiatively inefficient phase, because viscous turbulent heating and nuclear recombination are unbalanced by neutrino cooling (Kohri et al. 2005). This state transition is also accompanied by "freeze-out" 8 of weak interactions, leading to the winds being neutron-rich (Metzger et al. 2008a). Neutron-rich mater is shielded within the degenerate disk midplane, being ejected only once the disk radius has become large enough, and the neutrino luminosity low enough, that weak interactions no longer appreciably raise Y e in the outflow.
These early estimates were followed by two-dimensional, axisymmetric hydrodynamical models of the disk evolution, which show that, in the case of prompt BH formation, the electron fraction of the disk outflows lies in the range Y e ∼ 0.2-0.4 (Fernández and Metzger 2013;Just et al. 2015), sufficient to produce the entire mass range of r -process elements (Just et al. 2015;Wu et al. 2016). The total fraction of the disk mass which is unbound by these "viscously-driven" winds ranges from ∼ 5% for a slowly spinning BH, to ∼ 30% for high BH spin χ BH 0.95 (Just et al.

Fig. 6
Longer-lived remnants produce higher Y e disk wind ejecta and bluer kilonovae. Shown here is the mass distribution by electron fraction Y e of the disk wind ejecta, calculated for different assumptions about the lifetime, t collapse , of the central NS remnant prior to BH formation, from the axisymmetric α-viscosity hydrodynamical calculations of Metzger and Fernández (2014). A vertical line approximately delineates the ejecta with enough neutrons to synthesize lanthanide elements (Y e 0.25) generate a red kilonova from that with Y e 0.25 which is lanthanide-poor and will generate blue kilonova emission. The NS lifetime has a strong effect on the ejecta composition because it is a strong source of electron neutrinos, which convert neutrons in the disk to protons via the process ν e + n → p + e − . This figure is modified from a version in Lippuner et al. (2017) 2015; Fernández et al. 2015a); see also Kiuchi et al. (2015), who simulated the longterm evolution of BH-NS disks but without following the electron fraction evolution. These large disk ejecta fractions and neutron-rich ejecta were confirmed by the first 3D GRMHD simulations of the long-term disk evolution (Siegel andMetzger 2017, 2018;Fernández et al. 2019), with Siegel and Metzger (2017) finding that up to 40% of the initial torus may be unbound. The velocity and composition of magnetized disk outflows appears to be sensitive to the strength and geometry of the large-scale net magnetic flux threading the accretion disk Christie et al. 2019).
An even larger fraction of the disk mass (up to ∼ 90%) is unbound when the central remnant is a long-lived hypermassive or supramassive NS instead of a BH, due to the presence of a hard surface and the higher level of neutrino irradiation from the central remnant (Metzger and Fernández 2014;Fahlman and Fernández 2018). A longer-lived remnant also increases the electron fraction of the ejecta, which increases monotonically with the lifetime of the HMNS (Fig. 6). Most of the ejecta is lanthanidefree (Y e 0.3) if the NS survives longer than about 300 ms (Metzger and Fernández 2014;Kasen et al. 2015;Lippuner et al. 2017). Even when BH formation is prompt, simulations with Monte Carlo radiation transport included find that the earliest phases of disk evolution can produce at least a modest quantity of high-Y e material .
The mass ejected by the late disk wind can easily be comparable to, or larger than, that in the dynamical ejecta (e.g., Wu et al. 2016, their Fig. 1). Indeed, the total ejecta mass inferred for GW170817 greatly exceeds that of the dynamical ejecta Fig. 7 Different components of the ejecta from NS-NS mergers and the possible dependence of their kilonova emission on the observer viewing angle, θ obs , relative to the binary axis, in the case of a relatively prompt BH formation (left panel) and a long-lived magnetar remnant (right panel). In both cases, the dynamical ejecta in the equatorial plane is highly neutron-rich (Y e 0.1), producing lanthanides and correspondingly "red" kilonova emission peaking at NIR wavelengths. Mass ejected dynamically in the polar directions may be sufficiently neutron-poor (Y e 0.3) to preclude lanthanide production, powering "blue" kilonova emission at optical wavelengths (although this component may be suppressed if BH formation is extremely prompt). The outermost layers of the polar ejecta may contain free neutrons, the decay of which powers a UV transient lasting a few hours following the merger (Sect. 6.1.1). Re-heating of the ejecta by a delayed relativistic outflow (e.g., the GRB jet or a wind from the magnetar remnant) may also contribute to early blue emission (Sect. 6.1.2). The innermost ejecta layers originate from accretion disk outflows, which may emerge more isotropically. When BH formation is prompt, the disk wind ejecta is mainly neutron-rich, powering red kilonova emission (Fernández and Metzger 2013;Just et al. 2015;Wu et al. 2016;Siegel and Metzger 2017). If the NS remnant is instead long-lived relative to the disk lifetime, then neutrino emission can increase Y e sufficiently to suppress lanthanide production and result in blue disk wind emission ( Fig. 6; e.g., Metzger and Fernández 2014;Perego et al. 2014). Energy input from the central accreting BH or magnetar remnant enhance the kilonova luminosity compared to that exclusively from radioactivity (Sect. 6.2) found in merger simulations (Fig. 5), but is consistent in both its mass and velocity with originating from a disk wind (e.g., Siegel and Metzger 2017). As the disk outflows emerge after the dynamical ejecta, the disk outflow material will be physically located behind the latter (Fig. 7).
Beyond the dynamical and disk wind ejecta, other ejecta sources have been proposed, though these remain more speculative because the physical processes at work are less robust. Mass loss may occur from the differentially rotating NS during the process of angular momentum redistribution (Fujibayashi et al. 2018;Radice et al. 2018b). However, the details of this mechanism and its predictions for the ejecta properties depend sensitively on the uncertain physical source and operation of the "viscosity" currently put into the simulations by hand; unlike the quasi-Keplerian accretion disk on larger radial scales, the inner regions of the NS remnant possess a positive shear profile d /dr > 0 are therefore not unstable to the magneto-rotational instability.
Outflows can also occur from the HMNS/SMNS or stable NS remnant as it undergoes Kelvin-Helmholtz contraction and neutrino cooling over a timescale of seconds. At a minimum there will be outflows driven from the NS surface by neutrino heating (Dessart et al. 2009), which typically will possess a relatively low mass-loss rate 10 −3 M s −1 and low asymptotic velocity ∼ 0.1 c. However, the NS remnant possesses an ordered magnetic field of strength ∼ 10 14 −10 15 G, then the mass-loss rate  Kasen et al. (2013). Bound-free opacities are estimated as that of neutral iron (Verner et al. 1996), which should crudely approximate the behavior of heavier r -process elements. Electron scattering opacity accounts for the Klein-Nishina suppression at energies m e c 2 and (very schematically) for the rise in opacity that occurs above the keV energy scale due to all electrons (including those bound in atoms) contributing to the scattering opacity when the photon wavelength is smaller than the atomic scale. At the highest energies, opacity is dominated by pair creation by gamma-rays interacting with the electric fields of nuclei in the ejecta (shown schematically for Xenon, A = 131, Z = 54). Not included are possible contributions from r -process dust; or γ −γ pair creation opacity at photon energies m e c 2 ∼ 10 6 eV (see Eq. 9) and velocity of such an outflow is substantially enhanced by the centrifugal force along the open magnetic field lines (e.g., Thompson et al. 2004). ) argue that such a magnetar wind, from a HMNS of lifetime ∼ 0.1−1 s, was responsible for the fastest ejecta in GW170817 (Sect. 5). While the presence of an ordered magnetic field of this strength is physically reasonable, its generation from the smaller scale magnetic field generated during the merger process has yet to be conclusively demonstrated by numerical simulations.

Ejecta opacity
It is no coincidence that kilonova emission is centered in the optical/IR band, as this is the first spectral window through which the expanding merger ejecta becomes transparent. Figure 8 illustrates semi-quantitatively the opacity of NS merger ejecta near peak light as a function of photon energy. At the lowest frequencies (radio and far-IR), free-free absorption from ionized gas dominates, as shown with a red line in Fig. 8, and calculated for the approximate ejecta conditions three days post merger. As the ejecta expands, the free-free opacity will decrease rapidly due to the decreasing density ∝ ρ ∝ t −3 and the fewer number of free electrons as the ejecta cools and recombines.
At near-IR/optical frequencies, the dominant source of opacity is a dense forest of line (bound-bound) transitions. The magnitude of this effective continuum opacity is determined by the strengths and wavelength density of the lines, which in turn depend sensitively on the ejecta composition. If the ejecta contains elements with relatively simple valence electron shell structures, such as iron, then the resulting opacity is comparatively low (dashed brown line), only moderately higher than the Fe-rich ejecta in Type Ia supernovae (Pinto and Eastman 2000). On the other hand, if the ejecta also contains even a modest fraction of elements with partially-filled f-shell valence shells, such as those in the lanthanide and actinide group, then the opacity can be an order of magnitude or more higher Tanaka and Hotokezaka 2013;Fontes et al. 2015Fontes et al. , 2017Even et al. 2019). In both cases, the opacity rises steeply from the optical into the UV, due to the increasing line density moving to higher frequencies.
Considerable uncertainty remains in current calculations of the La/Ac opacities because the atomic states and line strengths of these complex elements are not measured experimentally. Theoretically, such high-Z atoms represent an unsolved problem in N-body quantum mechanics, with statistical models that must be calibrated to experimental data. Beyond identifying the line transitions themselves, there is considerably uncertainty in how to translate these data into an effective opacity. The commonly employed "line expansion opacity" formalism (Pinto and Eastman 2000), based on the Sobolev approximation and applied to kilonovae by Barnes and Kasen (2013) and Tanaka and Hotokezaka (2013), may break down if the line density is sufficiently high that the wavelength spacing of strong lines becomes comparable to the intrinsic thermal) width of the lines Fontes et al. 2015Fontes et al. , 2017. Nevertheless, the qualitative dichotomy between the opacity of La/Ac-free and La/Ac-bearing ejecta is likely to be robust and will imprint diversity in the kilonova color evolution (Sect. 4.1.2). 9 Despite the strong time-and wavelength-dependence of the opacity, for purposes of an estimate it is reasonable to model the kilonova using a constant effective "grey" (wavelength-independent) opacity, κ. Including a large range of r -process nuclei in their analysis, Tanaka et al. (2019) found (for temperatures T = 5−10 × 10 3 K characteristic of the ejecta near the time of peak emission) values of κ 20−30 cm 2 g −1 for Y e 0.2 (sufficient neutrons for the r -process to extend up to or beyond the third abundance peak at A ∼ 195 with a large lanthanide mass fraction), κ ∼ 3−5 cm 2 g −1 for Y e ≈ 0.25 − 0.35 (r -process extending only to the second abundance peak A ∼ 130 with a small or zero lanthanide mass fraction) and κ ∼ 1 cm 2 g −1 for Y e ≈ 0.40 (mainly neutron-rich Fe-group nuclei and a weak r -process). The approximate opacity range corresponding to different ejecta composition is summarized in Table 4.
Throughout the far UV and X-ray bands, bound-free transitions of the partially neutral ejecta dominates the opacity (blue line in Fig. 8). This prevents radiation from escaping the ejecta at these frequencies, unless non-thermal radiation from the central magnetar or BH remnant remains luminous enough to re-ionize the ejecta at late times (Sect. 6.2.2). Margalit et al. (2018) find that an X-ray luminosity L X 10 42 erg/s would be required to ionize M ej ∼ 10 −3 M of Fe-like ejecta expanding at ≈ 0.2 c at t ∼ 1 day after the merger. However, a substantially higher luminosity L X 10 44 −10 45 erg s −1 would be needed 10 to ionize the greater expected quantity of disk wind ejecta M ej ∼ 0.01-0.1 M , especially considering that the bound-free opacity of r -process nuclei will be even higher than Fe. Such high X-ray luminosities are too large to be powered by fall-back accretion onto the central remnant from the merger ejecta, but could be achievable from the rotational energy input from a long-lived magnetar remnant (Fig. 3, right panel).
At hard X-rays and gamma-ray energies, electron scattering, with Klein-Nishina corrections, provides an important opacity (which becomes highly inelastic at energies m e c 2 ). For gamma-ray energies m e c 2 , the opacity approaches a constant value κ Aγ ≈ α fs κ T (Z 2 /A) due to electron/positron pair creation on nuclei, where α fs 1/137, and A and Z are the nuclear mass and charge, respectively (e.g., Zdziarski and Svensson 1989). For r -process nuclei with Z 2 /A 10 − 20 this dominates inelastic scattering at energies 10 MeV. The low opacity 0.1 cm 2 g −1 in the ∼ MeV energy range implies that gamma-rays released by radioactive decay of r -process elements largely escape the ejecta prior to the optical/NIR peak without their energy being thermalized (Sect. 4.1).
Gamma-rays with very high energies (m e c 2 ) 2 /hν s ∼ 0.3TeV(hν s /1eV) −1 can also create electron/positron pairs by interacting with (more abundant) lower energy optical or X-ray photons of energy hν s m e c 2 . The γ − γ pair creation optical depth through the ejecta of radius R = vt is roughly given by where U rad L/(4π R 2 c) is the energy density of the seed photons of luminosity L. The fact that τ γ −γ 1 for characteristic thermal luminosities ∼ 10 40 −10 42 erg s −1 shows that ∼ GeV-TeV photons from a putative long-lived central engine (e.g., millisecond magnetar; Sect. 6.2.2) will be trapped for days to weeks after the merger. Prompt TeV emission from a NS-NS merger is thus unlikely to come from the merger remnant, but still could be generated within the relativistically-beamed GRB jet on much larger physical scales.

Unified toy model
Thermal emission following a NS-NS or BH-NS merger (a "kilonova", broadly defined) can be powered by a variety of different energy sources, including radioactivity and central engine activity (see Fig. 3 for a summary of different heating sources). This section describes a simple model for the evolution of the ejecta and its radiation, which we use to motivate the potential diversity of kilonova emission properties.
Though ultimately no substitute for full, multi-dimensional, multi-group radiative transfer, this 1D toy model does a reasonable job at the factor of a few level. Some sacrifice in accuracy may be justified in order to facilitate a qualitative understanding, given the other uncertainties on the mass, heating rate, composition, and opacity of the ejecta. Following the merger, the ejecta velocity structure approaches one of homologous expansion, with the faster matter lying ahead of slower matter . We approximate the distribution of mass with velocity greater than a value v as a power-law, where M is the total mass, v 0 ≈ 0.1 c is the average (∼ minimum) velocity. We adopt a fiducial value of β ≈ 3, motivated by a power-law fit to the dynamical ejecta in the numerical simulations of (Bauswein et al. 2013b). In general the velocity distribution derived from numerical simulations cannot be fit by a single power-law (e.g., Fig. 3 of Piran et al. 2013), but the following analysis can be readily extended to the case of an arbitrary velocity distribution.
In analogy with Eq. (6), radiation escapes from the mass layer M v on the diffusion timescale where κ v is the opacity of the mass layer v and in the second equality makes use of Eq. (10) with β = 3. Equating t d,v = t gives the mass depth from which radiation peaks for each time t, where t peak is the peak time for diffusion out of the whole ejecta mass, e.g., Eq. (7) evaluated for v = v 0 . Emission from the outer layers (mass M v < M) peaks first, while the luminosity of the innermost shell of mass ∼ M peaks at t = t peak . The deepest layers usually set the peak luminosity of the total light curve, except when the heating rate and/or opacity are not constant with depth (e.g., if the outer layers are free neutrons instead of r -process nuclei; Sect. 6.1.1).
As the ejecta expands, the radius of each layer of depth M v of mass δ M v evolves according to The thermal energy δ E v of the layer evolves according to where the first term accounts for losses due to PdV expansion in the radiationdominated ejecta. The second term in Eq. (14), accounts for radiative losses (the observed luminosity) and t lc,v = R v /c limits the energy loss time to the light crossing time (this becomes important at late times when the layer is optically thin). The third term in Eq. (14), accounts for sources of heating, including radioactivity (Q r ,v ; Sect. 4.1), a millisecond magnetar (Q mag ; Sect. 6.2.2) or fall-back accretion (Q fb ; Sect. 6.2.1). The radioactive heating rate, being intrinsic to the ejecta, will in general vary between different mass layers v. In the case of magnetar or accretion heating, radiation must diffuse from the central cavity through the entire ejecta shell (Fig. 7, right panel). One must in general also account for the evolution of the ejecta velocity (Eq. 10) due to acceleration by pressure forces. For radioactive heating, the total energy input Q r ,v dt is less than the initial kinetic energy of the ejecta (Metzger et al. 2010a;Rosswog et al. 2013;Desai et al. 2019), in which case changes to the initial velocity distribution (Eq. 10) are safely ignored. However, free expansion is not a good assumption when there is substantial energy input from a central engine. In such cases, the velocity v 0 of the central shell (of mass M and thermal energy E v 0 ) is evolved separately according to d dt where the source term on the right hand side balances the PdV loss term in the thermal energy equation (14), and R 0 is the radius of the inner mass shell. Equation (17) neglects two details: (1) special relativistic effects, which become important for low ejecta mass 10 −2 M and the most energetic magnetar engines (Zhang 2013;Gao et al. 2013;Siegel and Ciolfi 2016a, b); (2) the secondary shock driven through the outer ejecta layers by the nebula inflated by a long-lived central engine (e.g., Kasen et al. 2016) and its effects on the outer velocity distribution (e.g., Suzuki and Maeda 2017). Under the idealization of blackbody emission, the temperature of the thermal emission is is the total luminosity (summed over all mass shells). The radius of the photosphere R ph (t) is defined as that of the mass shell at which the optical depth The flux density of the source at photon frequency ν is given by where D is the source luminosity distance. We have neglected cosmological effects such as K-corrections, but these can be readily included.
For simplicity in what follows, we assume that the opacity κ v of each mass layer depends entirely on its composition, i.e. we adopt a temperature-independent grey opacity. The most relevant feature of the composition is the mass fraction of lanthanide or actinide elements, which in turn depends most sensitively on the ejecta Y e (Sect. 3.2). Following Tanaka et al. (2019, Table 4), one is motivated to take smoothly interpolating when necessary. We caution that these values were calculated by  for ejecta temperatures in the range 5−10 × 10 3 K, i.e. similar to those obtained close to peak light, and therefore may not be appropriate much earlier (the first hours) or later (nebular phase). The full emission properties are determined by solving Eq. (14) for δ E v , and hence L v , for a densely sampled distribution of shells of mass δ M v and velocity v > v 0 . When considering radioactive heating acting alone, one can fix the velocity distribution (Eq. 10). For an energetic engine, the velocity of the central shell is evolved simultaneously using Eq. (17).
As initial conditions at the ejection radius R(t = 0) ≈ 10−100 km, it is reasonable to assume that the initial thermal energy is comparable to its final kinetic energy, . If the ejecta expands freely from the site of ejection, the predicted light curves are largely insensitive to the details of this assumption because the initial thermal energy is anyways quickly removed by adiabatic expansion. Section 6.1.2 explores an exception to this rule when the ejecta is re-heated by shock interaction with a delayed outflow from the central engine.

R-process heating
At a minimum, the ejecta is heated by the radioactive decay of heavy r -process nuclei. This occurs at a rateQ where X r ,v is the r -process mass fraction in mass layer M v and e r is the specific heating rate. For neutron-rich ejecta (Y e 0.2), the latter can be reasonably well approximated by the fitting formula (Korobkin et al. 2012) where t 0 = 1.3 s and σ = 0.11 s are constants, and th,m is the thermalization efficiency (see below). Equation (22) Using Eqs. (7) and (8) the peak luminosity can be estimated as Given the reasonably large range in values allowed in NS-NS or BH-NS mergers, M ∼ 10 −3 −0.1 M , κ ∼ 0.5 − 30 cm 2 g −1 , v ∼ 0.1−0.3 c, one can have L peak ≈ 10 39 −10 42 erg s −1 .
The time dependence ofq r is more complicated for higher 0.2 Y e 0.4, with 'wiggles' caused by the heating rate being dominated by a few discrete nuclei instead of the large statistical ensemble present at low Y e (Korobkin et al. 2012;Martin et al. 2015). However, when averaged over a realistic Y e distribution, the heating rate on timescales of days-weeks (of greatest relevance to the peak luminosity; Eq. 8), is constant to within a factor of a few for Y e 0.4 (Lippuner and Roberts 2015;Wu et al. 2019b; see Fig. 10). The radioactive decay power is sensitive to various uncertainties in the assumed nuclear physics (nuclear masses, cross sections, and fission fragment distribution) at the factor of a few level (e.g., Wu et al. 2019b) 11 , a point we shall return to when discussing GW170817 (Sect. 5).
Radioactive heating occurs through a combination of β-decays, α-decays, and fission (Metzger et al. 2010b;Barnes et al. 2016;Hotokezaka et al. 2016). The thermalization efficiency, th,v , depends on how these decay products share their energies with the thermal plasma. Neutrinos escape from the ejecta without interacting; ∼ MeV gamma-rays are trapped at early times ( 1 day), but leak out at later times given the low Klein-Nishina-suppressed opacity ( Fig. 8; Hotokezaka et al. 2016;Barnes et al. 2016). β-decay electrons, α-particles, and fission fragments share their kinetic energy effectively with the ejecta via Coulomb collisions (Metzger et al. 2010b) and by ionizing atoms (Barnes et al. 2016). For a fixed energy release rate, the thermalization efficiency is smallest for β-decay, higher for α-decay, and the highest for fission fragments. The thermalization efficiency of charged particles also depends on the magnetic field orientation within the ejecta, since the particle Larmor radius is generally shorter than the mean free path for Coulomb interactions. Because the actinide yield around mass number A ∼ 230 varies significantly with the assumed nuclear mass model, Barnes et al. (2016) finds that the effective heating rate including thermalization efficiency can vary by a factor of 2 -6, depending on time. Barnes et al. (2016) find that the combined efficiency from all of these processes typically decreases from th,v ∼ 0.5 on a timescale of 1 day to ∼ 0.1 at t ∼ 1 week (their Fig. 13). In what follows, we adopt the fit provided in their Table 1, where t day = t/1 day, and {a v , b v , d v } are constants that will in general depend on the mass and velocity of the layer under consideration. For simplicity, we adopt fixed values of a v = 0.56, b v = 0.17, c v = 0.74, corresponding to a layer with M = 10 −2 M and v 0 = 0.1 c. As we shall discuss, the luminosity and color evolution of kilonovae encode information on the total quantity of r -process ejecta and, in principle, the abundance of lanthanide/actinide elements. However, insofar as the lanthanides cover the atomic mass range A ∼ 140-175, kilonova observations at peak light do not readily probe the creation of the heaviest elements, those near the third r -process peak (A ∼ 195) and the transuranic elements ( A 240).
One avenue for probing the formation of ultra-heavy elements is by the light curve's decay at late times, weeks to months after maximum light. At such late times the radioactive heating is often dominated by a few discrete isotopes with well-measured half-lifes (e.g., 223 Ra [t 1/2 = 11.4 days], 225 Ac [t 1/2 = 10.0 days], 225 Ra [t 1/2 = 14.9 days], 254 Cf [t 1/2 = 60.5 days]) which could produce distinctive features (e.g., bumps or exponential decay-like features) in the bolometric light curve of characteristic timescale ∼ t 1/2 (Zhu et al. 2018;Wu et al. 2019b), much in the way that the half-life of 56 Co is imprinted in the decay of Type Ia supernovae. The ability in practice to identify individual isotopes through this method will depend on accurate models for the ejecta thermalization entering the nebular phase (Kasen and Barnes 2019;Waxman et al. 2019) as well as dedicated broad-band, multi-epoch follow-up of nearby kilonovae in the NIR (where most the nebular emission likely emerges) with sensitive facilities like the James Webb Telescope (Kasliwal et al. 2019;Villar et al. 2018). Table 5 compiles all r -process isotopes with half-lives in the range 10-100 day (Wu et al. 2019b).
Gamma-ray decay lines from the r -process element decays escape the ejecta within days or less of the merger and could in principle be directly observed from an extremely nearby event 3-10 Mpc with future gamma-ray satellites Korobkin et al. 2019). A related, but potentially more promising near-term strategy is a gamma-ray search for remnants of past NS mergers in our Galaxy (Wu et al. 2019a;Korobkin et al. 2019). Among the most promising isotopes for this purpose is 126 Sn, which has several lines in the energy range 415-695 keV and resides close to the second r -process peak, because of its half-life t 1/2 = 2.3 × 10 5 yr is comparable to the ages of the most recent Galactic merger(s). Wu et al. (2019a) estimate that multiple remnants are detectable as individual sources by next-generation gamma-ray satellites with line sensitivities ∼ 10 −6 −10 −8 γ cm −2 s −1 .

Red kilonova: lanthanide-bearing ejecta
All NS-NS mergers, and the fraction of BH-NS mergers in which the NS is tidally disrupted before being swallowed by the BH, will unbind at least some highly neutronrich matter (Y e 0.25) capable of forming heavy r -process nuclei. This lanthanidebearing high-opacity material resides within the equatorially-focused tidal tail, or in more spherical outflows from the accretion disk (Fig. 7, top panel). The disk outflows will contain a greater abundance of low-Y e 0.25 material in NS-NS mergers if the BH formation is prompt or the HMNS phase short-lived (Fig. 7, top panel).
The left panel of Fig. 9 shows an example light curve of such a 'red' kilonova, calculated using the toy model assuming an ejecta mass M = 10 −2 M , opacity κ = 20 cm 2 g −1 , minimum velocity v 0 = 0.1 c, and velocity index β = 3, at an assumed distance of 100 Mpc. For comparison, dashed lines show light curves calculated from Barnes et al. (2016), based on a full one-dimensional radiative transfer calculation, for similar parameters. The emission is seen to peak at NIR wavelengths on a timescale of several days to a week at J and K bands (1.2 and 2.2 µm, respectively).
One notable feature of the light curves calculated using full radiative transfer is the significant suppression of the emission in the UV/optical wavebands URV due to the high lanthanide opacity. Here, the assumption of a gray opacity made in the toy model results in an overestimation of the UV flux relative to that found by the full radiative transfer calculation (Barnes et al. 2016). This difference results in part because the true line opacity increases strongly moving to higher frequencies due to the higher density of lines in the UV (Fig. 8).

Blue kilonova: lanthanide-free ejecta
In addition to the highly neutron-rich ejecta (Y e 0.30), some of the matter unbound from a NS-NS merger may contain a lower neutron abundance (Y e 0.30) and thus will be free of lanthanide group elements (e.g., Metzger and Fernández 2014;Perego et al. 2014;Wanajo et al. 2014). This low-opacity ejecta can reside either in the polar regions, due to dynamical ejection from the NS-NS merger interface, or in more isotropic outflows from the accretion disk (e.g., Miller et al. 2019). The quantity of high-Y e matter will be greatest in cases when BH formation is significantly delayed relative to the lifetime of the accretion disk due to the strong neutrino luminosity of the NS remnant (Fig. 7, right panel). The right panel of Fig. 9 shows a model otherwise identical to that in the left panel, but which assumes a lower opacity κ = 1 cm 2 g −1 more appropriate to lanthanide-free Fig. 9 Kilonova light curves in AB magnitudes for a source at 100 Mpc, calculated using the toy model presented in Sect. 4, assuming a total ejecta mass M = 10 −2 M , minimum velocity v 0 = 0.1 c, and gray opacity κ = 20 cm 2 g −1 . The left panel shows a standard "red" kilonova, corresponding to ejecta bearing lanthanide elements, while the right panel shows a "blue" kilonova poor in lanthanides (κ = 1 cm 2 g −1 ). Shown for comparison in the red kilonova case with dashed lines are models from Barnes et al. (2016) for v = 0.1 c and M = 10 −2 M . Depending on the relative speeds of the two components and the viewing angle of the observer, both red and blue emission components can be present in a single merger, originating from distinct portions of the ejecta (Fig. 7) ejecta. The emission now peaks at the visual bands R and I, on a timescale of about 1 day at a level 2-3 magnitudes brighter than the lanthanide-rich case. This luminous, fast-evolving visual signal was key to the discovery of the kilonova counterpart of GW170817 (Sect. 5).

Mixed blue + red kilonova
In general, the total kilonova emission can be thought of as a combination of distinct 'blue' and 'red' components. This is because both high-and low-Y e ejecta components can be simultaneously visible following a merger, particularly for viewing angles close to the binary rotation axis (Fig. 7). For viewers closer to the equatorial plane, the blue emission may in some cases be blocked by the high-opacity lanthanide-rich tidal ejecta . Thus, although the presence of a days to week-long NIR transient is probably a generic feature of all mergers, the early blue kilonova phase might only be visible or prominent in a fraction of events. On the other hand, if the blue component expands faster than the tidal ejecta (or the latter is negligibly low in mass-e.g., for an equal-mass merger), the early blue emission may be visible from a greater range of angles than just pole-on (e.g., Christie et al. 2019).
It has become common practice following GW170817 to model the total kilonova light curve by adding independent 1D blue (low κ) and red (high κ) models on top of one another (e.g., Villar et al. 2017), i.e. neglecting any interaction between the ejecta components. While extremely useful for obtaining qualitative insight, in detail this assumption will result in quantitative errors in the inferred ejecta properties (e.g., Kasen et al. 2017;Wollaeger et al. 2018). With well-sampled photometry (in both time and frequency), the total ejecta mass should be reasonably well-measured: once the ejecta has become effectively transparent at late times the bolometric luminosity directly traces the radioactive energy input. However, at early times when the ejecta is still opaque, the radial and angular structure of the opacity (i.e., lanthanide abundance X La ) can couple distinct ejecta components in a way not captured by combining two independent 1D models (e.g., Kasen et al. 2015). Radial dependence of X La (e.g., due to a Y e gradient) is straightforward to implement in the toy model through a mass shell-dependent value of κ v (Eq. 11). If X La increases with radius (i.e. if the "red" ejecta resides physically outside of the "blue"), then in principle even a low amount of red (high-opacity) fast material can 'reprocess' the radioactive luminosity generated from a much greater mass of blue (low-opacity) slower material residing behind it. Kawaguchi et al. (2018) found that this could in principle lead to an over-estimate the quantity of blue ejecta if one models the light curve by simply adding independent red and blue components. For these reasons, caution must be taken in naively adding 'blue' and 'red' models and a detailed analysis must take into account not just photometric light curve information, but also spectral features (in the above example, for instance, reprocessing by an outer thin layer of lanthanide-rich matter would generate strong UV line blanketing). The geometry of the blue and red components of the ejecta, and the observer viewer angle, are also in principle distinguishable by their relative levels of polarization (Covino et al. 2017;Matsumoto 2018;Bulla et al. 2019).  Table 6 summaries a few key properties of GW170817 as inferred from its GW/EM emission.

GW170817: the first LIGO NS-NS merger
The timeline of the discovery was recounted in the capstone paper written jointly by LIGO/Virgo and astronomers involved in the EM follow-up (Abbott et al. 2017c) and will not be recounted here. In limiting the scope of our discussion, we also do not address the host galaxy and environment of the merger and its implication for NS-NS merger formation channels Hjorth et al. 2017;Im et al. 2017;Levan et al. 2017;Pan et al. 2017), nor shall we discuss the plethora of other science opportunities the kilonova enabled (e.g., H0 cosmology; Abbott et al. 2017a). We also do not touch upon inferences about the GRB jet and its connection to the observed prompt gamma-ray and non-thermal afterglow emission (e.g., Bromberg et al. 2018;Kasliwal et al. 2017;Gottlieb et al. 2018;Murguia-Berthier et al. 2017a;Salafia et al. 2018;Xiao et al. 2017), though possible connections between the jet and the early kilonova emission will be discussed in Sect. 6.1.2.

The kilonova
AT2017gfo started out blue in color, with a featureless thermal spectrum that peaked at UV frequencies (e.g., Nicholl et al. 2017;McCully et al. 2017;Evans et al. 2017), before rapidly evolving over the course of a few days to become dominated by emission with a spectral peak in the near-infrared (NIR) Pian et al. 2017;Tanvir et al. 2017). Although early blue colors are not uncommon among astrophysical transients (most explosions start hot and thereafter cool from expansion), the very fast evolution of AT2017gfo was completely unlike that seen in an previously known extragalactic event, making its connection to GW170817 of high significance (even before folding in theoretical priors on the expected properties of kilonovae). Simultaneous optical ( In discussing the interpretation of AT2017gfo, we start with the most basic and robust inferences that can be made, before moving onto areas where there is less universal agreement. Perhaps the first question one might ask is: What evidence exists that AT2017gfo was powered by r -process heating? and, if so, How much radioactive material was synthesized? Figure 10 from Wu et al. (2019b) shows the bolometric luminosity L bol (t) compiled from observations in the literature (Smartt et al. 2017;Cowperthwaite et al. 2017;Waxman et al. 2018;Arcavi 2018) compared to several distinct models for the time-dependent heating rate of r -process decayQ r (Eq. 21), in which the authors have varied the mean Y e of the ejecta contributing to the heating and the nuclear mass model, the latter being one of the biggest nuclear physics uncertainties.
A first takeaway point is the broad similarity between the observed L bol (t) evolution and the power-law-like decay predicted by the decay of a large ensemble of r -process isotopes (Metzger et al. 2010b). Furthermore, the total ejecta mass one requires to match the normalization of L bol varies with the assumptions, ranging from ≈ 0.02 M (Y e = 0.15; DZ31 mass model) to 0.06 M (Y e = 0.35; FRDM mass model). This range broadly agrees with that reported by independent groups modeling GW170817 (see Côté et al. 2018 for a compilation). It is also entirely consistent with the range of ejecta masses predicted from NS-NS mergers (Sect. 3.1), as we elaborate further below. Although non-r -process powered explanations for AT2017gfo can be constructed (e.g., invoking magnetar power; Sect. 6.2.2), they require several additional assumptions and thus are disfavored by Ockham's razor.
If the yield of r -process elements in GW170817 is at all representative of that of NS-NS mergers in the Universe as a whole (as the similarity of its GW-inferred properties compared to the Galactic binary NS population support it being; e.g., Zhao and Lattimer 2018), then, even adopting the lowest NS-NS merger rate currently allowed from LIGO/Virgo of ∼ 100 Gpc −1 yr −1 , an order-of-magnitude estimate (Eq. 3) leads to the conclusion that NS-NS mergers are major sources of r -process elements in the universe (e.g., Kasen et al. 2017;Côté et al. 2018). However, given large current uncertainties on the Galactic rate of NS-NS mergers and the precise abundance distribution synthesized in GW170817, it cannot yet be established that mergers are the exclusive, or even dominant, r -process site (see discussion in Sect. 2.1).
With the production and ejection of at least a few hundredths of a solar mass of neutron-rich elements established, the next question is the detailed nature of its composition. Specifically, which r -process elements were formed? Figure 11 shows a compilation of photometric data from the literature on AT2017gfo by Villar et al. (2017). The blue/UV bands (e.g., F225W, F275W) fade rapidly from the first observation at 11 hours, while the NIR bands (e.g., JHK) show a much flatter decay over the first week. The early-time blue emission suggests that the outermost layers of the merger ejecta (at least those dominating the observed emission) are composed of light r -process material with a low opacity (blue kilonova; Sect. 4.1.2) synthesized from merger ejecta with a relatively high 12 electron fraction, Y e 0.25. The more persistent late NIR emission instead requires matter with higher opacity, consistent with the inner ejecta layers containing at least a moderate amount of lanthanide or actinide elements (red kilonova; Sect 4.1.1).
Motivated by the theoretical prediction of distinct lanthanide-free and lanthaniderich ejecta components (e.g., Metzger and Fernández 2014), many groups interpreted AT2017gfo using mixed models described in the previous section, comprised of 2 or 3 separate ejecta components with different lanthanide abundances (e.g., Kasen et al. 2017;Tanaka et al. 2017;Drout et al. 2017;Kasliwal et al. 2017;Perego et al. 2017; however, see Waxman et al. 2018). As one example, the solid lines in Fig. 11 show a best-fit model from Villar et al. (2018) based on the sum of three spherical gray-opacity kilonova models ("blue", "purple", "red") with respective opacities κ = (0.5, 3, 10) cm 2 g −1 (similar to those given in Table 4) and from which they infer for the respective components ejecta masses M ej ≈ (0.02, 0.047, 0.011) M and mean velocities v ej ≈ (0.27, 0.15, 0.14)c. Mapping the opacities back to electron fractions using e.g., Table 4, one infers that most of the ejecta possessed intermediate values of Y e ≈ 0.25 − 0.35 which generated elements up to the second r -process peak. Smaller quantities of the ejecta had Y e 0.4 or Y e 0.25, the latter containing a sufficient neutron abundance to produce some lanthanide elements (A 140) if not nuclei 12 Rosswog et al. (2018) argue that the ejecta must have possessed Y e 0.3−0.35 to produce a smooth radioactive heating rate consistent with the bolometric light curve (due to the fact that discrete isotopes dominate the heating rate for higher Y e ; e.g., Lippuner and Roberts 2015). However, this makes the overrestrictive assumption that all the ejecta contains a single precise value of Y e . A small but finite spread in Y e about the mean valueȲ e results in a smooth light curve decay consistent with observations even for   Villar et al. (2017), copyright by the authors extending up to third r -process peak (A ∼ 195) or beyond. Despite the unprecedented data set available for AT2017gfo, it is unfortunately not possible to reconstruct the detailed abundance pattern synthesized, for instance to test its consistency with that observed in metal-poor stars or in our solar system (e.g., Hotokezaka et al. 2018).
From the inferred ejecta mass, velocity, and Y e distribution, the next question is: What phase or phases during or following the merger was this material released? One thing is clear: the dynamical ejecta alone is insufficient. Fig. 5 shows that the total ejecta mass 0.02 M exceeds the predicted dynamical ejecta from essentially all NS-NS merger simulations published to date, while the average velocity of the bulk of the lower-Y e ejecta ≈ 0.1 c is also significantly less than predicted for the dynamical ejecta. The bulk of the ejecta, particularly the redder low-Y e component, is instead most naturally explained as an outflow from the remnant accretion torus created around the central compact object following the merger (Sect. 3.1.2). GRMHD simulation of the post-merger disk evolution demonstrate that ≈ 40% of the initial mass of the torus (i.e. up to ∼ 0.08 M in wind ejecta for initial disk masses up to ≈ 0.2 M predicted by simulations) is unbound at an average velocity of v ≈ 0.1 c (e.g., Siegel and Metzger 2017;Fernández et al. 2019). The disk wind ejecta can contain a range of electron fractions (and thus produce blue or red emission), depending e.g., on the lifetime of the central NS remnant prior BH formation (see Fig. 6).
The physical source of the lanthanide-poor ejecta (Y e 0.35) responsible for powering the early-time blue emission is more open to debate. Table 7 summarizes  (Dietrich et al. 2017b;Dietrich and Ujevic 2017;Gao et al. 2017) b However, a small NS radius may be in tension with the creation of a large accretion disk needed to produce the red KN ejecta (Radice et al. 2018c) several possible origins, along with some of their pros and cons. The high velocities v blue ≈ 0.2−0.3 c and composition (Y e 0.25) broadly agree with predictions for the shock-heated dynamical ejecta (e.g., Oechslin and Janka 2006;Sekiguchi et al. 2016;Radice et al. 2016b). However, one concern is the large quantity M blue 10 −2 M , which again is higher than the total predicted by most merger simulations (Fig. 5), especially considering that only a fraction of the dynamical ejecta will possess a high Y e . If the blue ejecta is dynamical in origin, this could provide evidence for a small value of the NS radius (Nicholl et al. 2017, Sect. 5.2) because the quantity of shock-heated ejecta appears to grow with the NS compactness (Bauswein et al. 2013b).
The highest velocity tail of the kilonova ejecta might not be an intrinsic property, but instead the result of shock-heating of an originally slower ejecta cloud by a relativistic jet created following some delay after the merger (e.g., Bucciantini et al. 2012;Duffell et al. 2015;Gottlieb et al. 2018;Kasliwal et al. 2017;Bromberg et al. 2018;Piro and Kollmeier 2018;Sect. 6.1.2). However, even models which invoke an early component of cocoon emission ) require a radioactive-powered component of emission which dominates after the first 11 hours that contains a large mass ∼ 10 −2 M of low-opacity (high Y e ) matter (stated another way, the observations do not require an early extra emission component beyond radioactivity; see Fig. 12). Further disfavoring a jet-related origin for the blue kilonova ejecta is that the kinetic energy of the latter M blue v 2 blue /2 ≈ 10 51 ergs exceeds the kinetic energies of cosmological short gammaray bursts (≈ 10 49 −10 50 erg ;Nakar 2007;Berger 2014), and that of the off-axis jet specifically required to fit the off-axis afterglow of GW170817 (e.g., Margutti et al. 2018), by a large factor ∼ 10 − 100. Metzger et al. (2018) proposed an alternative source for the blue kilonova ejecta: a magnetized wind which emerges from the HMNS remnant ≈ 0.1−1 s prior to its collapse to a BH. While the HMNS remnant was proposed as a potential ejecta source for GW170817, (e.g., Evans et al. 2017), the velocity and mass-loss rate of such purely neutrino-powered winds (Dessart et al. 2009;Perego et al. 2014) are insufficient to explain that of the observed kilonova. Metzger et al. (2018) emphasize the role that Fig. 12 Bolometric kilonova light curve during the first few hours of a NS-NS merger, calculated for several model assumptions that can reproduce the measured luminosity L bol ≈ 10 42 erg s −1 of AT2017gfo at t ≈ 11 hr (blue uncertainty bar; e.g., Arcavi et al. 2017b;Cowperthwaite et al. 2017;Drout et al. 2017). Black solid lines show how r -process only models change with the assumed timescale t 0 = 0.01, 0.1, 1 s at which the outer ejecta was last "thermalized", i.e. endowed with an internal thermal energy comparable to its asymptotic kinetic energy (at t t 0 , the ejecta is heated is solely by r -process radioactivity in these models). A small value of t 0 ∼ 0.01 s corresponds to a dynamical ejecta origin with no additional heating, while a large value of t 0 ∼ 0.01−1 s represents the case of a long-lived engine (GRB jet, magnetar wind or accretion disk outflow) which re-heats the ejecta on a timescale ∼ t 0 . We adopt parameters β = 3, v 0 = 0.25 c, M = 0.025 M , κ = 0.5 cm 2 g −1 (except for the t 0 = 1 s case, for which M = 0.015 M ). Red dashed lines show models with t 0 = 0.01 s but for which the outer layer of mass M n is assumed to contains free neutrons instead of r-process nuclei (a model similar to those shown in Fig. 14). Note that the early-time signatures of neutron decay are largely degenerate with late-time shock re-heating of the ejecta. Image reproduced with permission from Metzger et al. (2018), copyright by the authors strong magnetic fields have on increasing the mass-loss rate and velocity of the wind through centrifugal slinging (similar to models of magnetized winds from ordinary stars; Belcher and MacGregor 1976). Using a series of 1D wind models, they found that a temporarily-stable magnetar remnant with a surface field strength B ≈ 1−3×10 14 G can naturally produce the mass, velocity, and composition of the blue kilonova ejecta in GW170817. We return to the role that much longer-lived magnetar remnants can play in the kilonova emission in Sect. 6.2.2.
Recently, using numerical relativity simulations which include approximate neutrino transport and a treatment of the effects of turbulent viscosity in the disk, Nedora et al. (2019) found that spiral density waves generated in the post-merger accretion disk by the central HMNS remnant can lead to the ejection of ∼ 10 −2 M in matter with Y e 0.25 and velocity ∼0.15-0.2 c. An open question is whether such spiral waves behave similarly, and can produce ejecta with sufficiently high velocities 0.2 c to explain AT2017gfo, even in the physical case in which the magneto-rotational instability operates simultaneously in the disk. If so, this mechanism would provide an additional promising source of blue kilonova ejecta from a moderately long-lived HMNS.

Inferences about the neutron star equation of state
GW170817 provided a wealth of information on a wide range of astrophysical topics. One closely connected to the focus of this review are new constraints it enabled on the equation of state (EOS) of nuclear density matter, that which is responsible for determining the internal structure of the NS and setting its key properties, such as its radius and maximum stable mass M TOV (Lattimer and Prakash 2016;Özel and Freire 2016).
Even absent a bright EM counterpart, the gravitational waveform can be used to measure or constrain the tidal deformability of the inspiraling stars prior to their disruption resulting from tidal effects on the inspiral phase evolution (e.g., Raithel et al. 2018;De et al. 2018;Abbott et al. 2018). Assuming two stars with the same EOS, observations of GW170817 were used to place limits on the radius of a 1.6 M NS of R 1.6 = 10.8 +2.0 −1.7 km (Abbott et al. 2018). Likewise, in BH-NS mergers, measurement of tidal interactions and the cut-off GW frequency at which the NS is tidally disrupted by the BH, provide an alternative method to measure NS radii (e.g., Kyutoku et al. 2011;Lackey et al. 2014;Pannarale 2013;Pannarale et al. 2015).
Unfortunately, current generation of GW detectors are far less sensitive to the postmerger signal and thus of the ultimate fate of the merger remnant, such and whether and when a BH is formed (as was true even for the high signal-to-noise event, GW170817; Abbott et al. 2017e). Here, EM observations provide a complementary view. In a BH-NS merger, the presence or absence of an EM counterpart is informative about whether the NS was tidally disrupted and thus can be used to measure its compactness (e.g., Ascenzi et al. 2019). As discussed in Sect. 3.1 and summarized in Table 3, the type of compact remnant which is created by a NS-NS merger (prompt collapse, HMNS, SMNS, or stable NS) depends sensitively on the total binary mass M tot relative to various threshold masses, which depend on unknown properties of the EOS, particularly M TOV and R 1.6 (Fig. 13). Thus, if one can infer the type of remnant produced in a given merger from the EM counterpart, e.g., the kilonova or GRB emission, then by combining this with GW measured 13 value of M tot one can constrain the values of M TOV and/or R 1.6 .
In GW170817, the large quantity of ejecta 0.02 M inferred from the kilonova, and its high electron fraction, strongly disfavored that the merger resulted in a prompt (∼ dynamical timescale) collapse to a BH. Given that the threshold for prompt collapse depends on the NS compactness (Bauswein et al. 2013a), this enabled Bauswein et al. (2017) to place a lower limit of R 1.6 10.3−10.7 km (depending on the conservativeness of their assumptions). Radice et al. (2018c) came to a physicallyrelated conclusion (that GW170817 produced a large ejecta mass not present in the case of prompt BH formation), but expressed their results as a lower limit on tidal deformability instead of R 1.6 (see also Coughlin et al. 2019 for a joint Bayesian analysis of the EM and GW data).
Going beyond the inference that GW170817 initially formed a NS remnant instead of a prompt collapse BH to infer the stability and lifetime of the remnant becomes trick-  Fig. 13 The four possible outcomes of a NS-NS merger depend on the total binary mass relative to various threshold masses, each of which is proportional to the maximum mass, M TOV , of a non-rotating NS (Table 3). Prompt BH formation or a short-lived HMNS results generates ejecta with a relatively low kinetic energy ∼ 10 50 −10 51 erg (energy stored in the differential rotation of the HMNS remnant can largely be dissipated as heat and thus lost to neutrinos). By contrast, the delayed formation of a BH through spin-down of a SMNS or stable remnant takes place over longer, secular timescales and must be accompanied by the release of substantial rotational energy ∼ 10 52 −10 53 erg. Unless effectively "hidden" through GW emission, a large fraction of this energy will be transferred to the ejecta kinetic energy (and, ultimately, the ISM forward shock), thus producing a more luminous kilonova and synchrotron afterglow than for a short-lived remnant. ier. Nevertheless, several independent arguments can be made which taken together strongly suggest that remnant was a relatively short-lived HMNS (t collapse 0.1−1 s), rather than a SMNS or indefinitely-stable NS (Margalit and Metzger 2017;Granot et al. 2017;Bauswein et al. 2017;Perego et al. 2017;Rezzolla et al. 2018;Ruiz et al. 2018;Pooley et al. 2018).
-The presence of a significant quantity of lanthanide-rich disk wind ejecta, as inferred from the presence of red kilonova emission, is in tension with the Y e distribution predicted for a merger remnant that were to have survived longer than several hundred milliseconds (Metzger and Fernández 2014;Lippuner et al. 2017; see Fig. 6). -The kinetic energies of the kilonova ejecta (∼ 10 51 erg) and the off-axis gammaray burst jet inferred from the X-ray/radio afterglow (∼ 10 49 −10 50 erg) exceed by a large factor 10−100 the rotational energy necessarily released for a SMNS or stable NS remnant to collapse to a BH ( Fig. 13; see further discussion in Sect. 6.2.2). -The formation of an ultra-relativistic GRB jet on a timescale of 1 s after the merger is believed by many to require a clean polar funnel only present above a BH (e.g., Murguia-Berthier et al. 2017b; see Sect. 6.2.2 for discussion of this point).
All indications from the GW170817 afterglow point to the presence of an off-axis jet with properties consistent with the cosmological short GRB population (e.g., Wu and MacFadyen 2019). -The magnetar spin-down luminosity could power temporally-extended X-ray emission minutes to days after the merger (Sect. 6.2.2); however, the observed X-rays from GW170817 are completely explained by the GRB afterglow without excess emission from a long-lived central remnant being required Pooley et al. 2018; however, see Piro et al. 2019).
Taking the exclusion of a SMNS remnant in GW170817 for granted, and combining this inference with the measured binary mass M tot = 2.74 +0.04 −0.01 M (Abbott et al. 2017b) from the GW signal, Margalit and Metzger (2017) place an upper limit on the TOV mass of M TOV 2.17 M (see also Shibata et al. 2017;Rezzolla et al. 2018;Ruiz et al. 2018; however see , who argue for a more conservative constraint of M TOV 2.30 M ). Stated another way, if M TOV were much higher than this limit, one would expect the remnant of GW170817 to have survived longer and produced an EM signal markedly different than the one observed. If this result holds up to further scrutiny, it provides the most stringent upper limit on M TOV currently available (and one in possible tension with the high NS masses ≈ 2.4 M suggested for some so-called "black widow" pulsars; e.g., Romani et al. 2015).
The above methods for constraining the NS EOS from pure GW or joint EM/GW data come with systematic uncertainties (most yet to be quantified), albeit different ones than afflict current EM-only methods. However, there is reason to hope these will improve with additional modeling and observations. For instance, if our current theoretical understanding of the diversity of possible outcomes of NS-NS mergers with the in-going binary properties (M tot , q) is correct (Fig. 4), these predicted trends should be verifiable from a sample of future kilonova/afterglow observations (Sect. 7.3). Margalit and Metzger (2019) show that ∼ 10 joint EM-GW detections of sufficient quality to accurately ascertain the merger outcome could constrain the values of R 1.6 and M TOV to the several percent level where systematic effects are certain to dominate the uncertainties. The future is bright for EM-GW joint studies of the NS EOS in tandem with improvements in our ability to understand and even predict the diverse outcomes of NS-NS or BH-NS merger events (Sect. 7.3).

Diversity of kilonova signatures
The previous section described the most "vanilla" models of red/blue kilonovae powered exclusively by r -process heating and how the thermal UVOIR emission following GW170817 could be adequately described in this framework. This section explores additional, sources of emission which are theoretically predicted by some models but have not yet been observed (at least unambiguously). Either these emission sources were not accessible in GW170817 due to observational limitations, or they could be ruled out in this event but nevertheless may accompany future NS-NS or NS-BH mergers (e.g., for different masses of the in-going binary stars). Though some of these possibilities remain speculative, their consideration is nevertheless useful to define future observational goals and to inform search strategies regarding just how differently future mergers could appear than GW170817.

The first few hours
The UV/optical counterpart of GW170817 was discovered roughly 11 hours after the two stars merged. This delay was largely due to the event taking place over the Indian Ocean, rendering its sky position initially inaccessible to the majority of ground-based follow-up telescopes. However, roughly half of future mergers should take place in the northern hemisphere (above the LIGO detectors) which improves the chances of rapid optical follow-up, potentially within hours or less from the time of coalescence (e.g., Kasliwal and Nissanke 2014). A future wide-field UV satellite (or fleet of satellites) able to rapidly cover GW event error regions could revolutionize the early-time frontier.
Kilonovae are powered by the outwards diffusion of thermal radiation. The earliest time emission therefore probes the fastest, outermost layers of the ejecta. In the most simple-minded and conservative scenario, these layers are also heated by r -process decay, rendering the early-time emission a simple continuation of the r -process kilonova to earlier times. However, due to the higher temperatures when the ejecta is more compact, the ionization states of the outer layers (and hence which elements dominate the line opacity) could differ markedly from those at later times. Full opacity calculations which extend to thermodynamic conditions appropriate to the first few hours of the transient are a necessary ingredient to making more accurate predictions for this early phase (e.g., Tanaka et al. 2019).
However, it is also possible that the outermost layers of the ejecta are even hotterand thus more luminous-than expected from r -process heating alone. Additional sources of early-time heating include: (1) radioactive decay of free neutrons which may be preferentially present in the fast outers layers (Sect. 6.1.1) or (2) the delayed passage through the ejecta by a relativistic jet or wide-angle outflow (Sect. 6.1.2). The tail-end of such an extra emission component could in principle have contributed to the earliest epochs of optical/UV emission from GW170817 (Arcavi 2018), though the available data is fully consistent with being powered exclusively by r -process heating.

Neutron precursor emission
The majority of the ejecta from a NS-NS merger remain sufficiently dense during its decompression from nuclear densities that all neutrons are captured into nuclei during the r -process (which typically takes place seconds after matter is ejected). However, some NS-NS merger simulations find that a small fraction of the dynamical ejecta (typically a few percent, or ∼ 10 −4 M ) can expand sufficiently rapidly that the neutrons in the ejecta do not have time to be captured into nuclei (Bauswein et al. 2013b), i.e., the r -process "freezes out". In the simulations of Bauswein et al. (2013b) this fast-expanding matter, which reaches asymptotic velocities v 0.5 c, originates from the shocked interface between the merging stars and resides on the outermost layers of the polar ejecta (see also Ishii et al. 2018). Equally fast-expanding material could in principle be produced via other mechanisms which take place after the dynamical phase, e.g., the passage through the ejecta by a GRB jet or in a magnetized wind from the HMNS remnant (see next section).
Free neutrons, if present in the outer ejecta layers, provide an order of magnitude greater specific heating rate than produced by r -process nuclei on timescales of tens of minutes to hours (Fig. 3). Metzger et al. (2015) emphasized that such super-heating by such a free neutron layer could substantially enhancing the early kilonova emission (see also Kulkarni 2005).
An ejecta layer δ M v containing free neutrons experiences a radioactive heating rate of  Fig. 9. The left panel shows a calculation with an opacity appropriate to lanthanide-bearing nuclei, while the right panel shows an opacity appropriate to lanthanidefree ejecta. Models without a free neutron layer (M n = 0; Fig. 9) are shown for comparison with dashed linesQ where the initial mass fraction of neutrons, is interpolated in a smooth (but otherwise ad-hoc) manner between the neutron-free inner layers at M M n and the neutron-rich outer layers M M n , which have a maximum mass fraction of 1 − 2Y e . The specific heating rate due to neutron β-decay (accounting for energy loss to neutrinos) is given bẏ e n = 3.2 × 10 14 exp[−t/τ n ] erg s −1 g −1 , where τ n ≈ 900 s is the neutron half-life. The rising fraction of free neutrons in the outermost layers produces a corresponding decreasing fraction of r -process nuclei in the outermost layers, i.e., X r ,v = 1 − X n,v in calculating the r -process heating rate from Eq. (21). Figure 14 shows kilonova light curves, including an outer layer of neutrons of mass M n = 10 −4 M and electron fraction Y e = 0.1. In the left panel, we have assumed that the r -process nuclei which co-exist with the neutrons contain lanthanides, and hence would otherwise (absent the neutrons) produce a "red" kilonova. Neutron heating boosts the UVR luminosities on timescales of hours after the merger by a large factor compared to the otherwise identical case without free neutrons (shown for comparison with dashed lines). Even compared to the early emission predicted from otherwise lanthanide-free ejecta ("blue kilonova"), neutron decay increases the luminosity during the first few hours by a magnitude or more, as shown in the right panel of Fig. 14.
How can such a small layer of neutrons have such a large impact on the light curve? The specific heating rate due to free neutronsė n (Eq. 28) exceeds that due to r -process nucleiė r (Eq. 22) by over an order of magnitude on timescales ∼ 0.1−1 hr after the merger. This timescale is also, coincidentally, comparable to the photon diffusion depth of the inner edge of the neutron mass layer if M n 10 −5 M . Indeed, setting t d,v = t in Eq. (11), the emission from mass layer M v peaks on a timescale The total energy energy released by neutron-decay is E n ė n M n dt ≈ 6 × 10 45 (M n /10 −5 M ) erg for Y e 0.5. Following adiabatic losses, a fraction τ n /t peak,v ∼ 0.01−0.1 of this energy is available to be radiated over a timescale ∼ t peak,v . The peak luminosity of the neutron layer is thus approximately and hence is relatively insensitive to the mass of the neutron layer, M v = M n . This peak luminosity can be ∼ 10−100 times higher than that of the main r -process powered kilonova peak. The high temperature of the ejecta during the first hours of the merger will typically place the spectral peak in the UV, potentially even in cases when the free neutron-enriched outer layers contain lanthanide elements.
Additional theoretical and numerical work is needed to assess the robustness of the fast-moving ejecta and its abundance of free neutrons, which thus far has been seen in a single numerical code (Bauswein et al. 2013b). The freeze-out of the r -process, and the resulting abundance of free neutrons, is also sensitive to the expansion rate of the ejecta (Lippuner and Roberts 2015), which must currently be extrapolated from the merger simulations (running at most tens of milliseconds) to the much longer timescales of ∼ 1 second over which neutrons would nominally be captured into nuclei. Figure 14 and Eq. (30) make clear that the neutron emission is sensitive to the opacity of the ejecta at early stages, when the temperatures and ionization states of the ejecta are higher than those employed in extant kilonova opacity calculations.

Shock re-heating (short-lived engine power)
The ejecta from a NS merger is extremely hot 10 10 K immediately after becoming unbound from the central remnant or accretion disk. However, due to PdV losses, the temperature drops rapidly ∝ 1/R as the ejecta radius R = vt expands. Absent additional sources of heating, the internal energy decays in time as e 0 (t) 0.68 ρ 1/3 s 4/3 a 1/3 ≈ 4 × 10 12 erg g −1 s where the ejecta density has been set to its mean value ρ = 3M/(4π R 3 ) and the entropy s normalized to a value 20 k b per baryon typical of the shock-heated polar dynamical or disk wind ejecta (the unshocked tidal tail material can be even colder). As the ejecta expands, it receives heating from the decay of r -process nuclei at the rate given by Eq. (23). The thermal energy input from r -process decay on timescales ∼ t is thus approximately given by e r ∼ė r t ≈ 9 × 10 14 t 1 day The key point to note is that e r e 0 within minutes after the merger. This demonstrates why the initial thermal energy of the ejecta can be neglected in calculating the kilonova emission on timescales of follow-up observations of several hours to days (or, specifically, why the toy model light curves calculations presented thus far are insensitive to the precise initial value of E v ).
However, the early-time luminosity of the kilonova could be substantially boosted from this naive expectation if the ejecta is re-heated at large radii, well after its initial release (i.e. if the ejecta entropy is boosted to a substantially larger value than assumed in Eq. 31).
One way such re-heating could take place is by the passage of a relativistic GRB jet through the polar ejecta, which generates a shocked "cocoon" of hot gas (e.g., Gottlieb et al. 2018;Kasliwal et al. 2017;Piro and Kollmeier 2018). However, the efficiency of this heating process is debated. Duffell et al. (2018) found, using a large parameter study of jet parameters (jet energies ∼ 10 48 −10 51 erg and opening angles θ ∼ 0.07−0.4 covering the range thought to characterize GRBs) that the thermal energy deposited into the ejecta by the jet falls short of that produced by the r -process heating on the same timescale by an order of magnitude or more. Jet heating is particularly suppressed when the relativistic jet successfully escapes from the ejecta, evidence in GW170817 by late-time afterglow observations Mooley et al. 2018).
An alternative means to shock-heat the ejecta is by a wind from the magnetized central NS remnant prior to its collapse into a BH . Such a wind is expected to have a wide opening angle and to accelerate to trans-relativistic speeds over a characteristic timescale of seconds (e.g., Metzger et al. 2008b). 14 This delay in the wind acceleration, set by the Kelvin-Helmholtz cooling of the remnant, would naturally allow the dynamical ejecta time to reach large radii before being hit and shocked by the wind. Although this has not yet been explored in the literature, even time variability in the accretion disk outflows (Sect. 3.1.2) could generate internal shocks and re-heat the wind ejecta over timescales comparable to the disk lifetime seconds (e.g., Fernández et al. 2019).
In all of these mechanisms, re-setting of the ejecta thermal energy at large radii is key to producing luminous emission, because otherwise the freshly-deposited energy is degraded by PdV expansion before being radiated. This is illustrated by Fig. 12, where black lines show how the early-time kilonova light curve is enhanced when the ejecta is re-thermalized (its thermal energy re-set to a value comparable to its kinetic energy, i.e. E v ∼ δ M v v 2 /2) at different times, t 0 , following its initial ejection. As expected, larger t 0 (later re-thermalization) results in more luminous emission over the first few hours.
But also note that the light curve enhancement from jet/wind re-heating looks broadly similar that resulting from the outer layers being composed of free neutrons (shown for comparison as red lines in Fig. 12). This degeneracy between free neutrons and delayed shock-heating makes the two physical processes challenging to observationally distinguish (Arcavi 2018;Metzger et al. 2018). Well-sampled early-time light curves, e.g., to search for a subtle bump in the light curve on the neutron half-live τ β ≈ 10 3 s, could be necessary to make progress on the interpretation. Regardless, the first few hours of the kilonova is an important frontier for future EM follow-up efforts: the signal during this time is sensitive to the origin of the ejecta and how it interacts with the central engine (the details which are largely washed out at later times when r -process heating takes over).

Long-lived engine power
The end product of a NS-NS or BH-NS merger is a central compact remnant, either a BH or a massive NS (Sect. 3.1). Sustained energy input from this remnant can produce an additional source of ejecta heating in excess of the minimal contribution from radioactivity, thereby substantially altering the kilonova properties (e.g., Yu et al. 2013;Metzger and Piro 2014;Wollaeger et al. 2019).
Evidence exists for late-time central engine activity following short GRBs, on timescales from minutes to days. A fraction ≈ 15−25% of Swift short bursts are followed by a distinct hump of hard X-ray emission lasting for tens to hundreds of seconds following the initial prompt spike (e.g., Norris and Bonnell 2006;Perley et al. 2009;Kagawa et al. 2015). The isotropic X-ray light curve of such temporally extended emission in GRB 080503 is shown in the bottom panel of Fig. 3  ). Other GRBs exhibit a temporary flattening or "plateau" in their X-ray afterglows lasting ≈ 10 2 −10 3 s (Nousek et al. 2006), which in some cases abruptly ceases (Rowlinson et al. 2010). X-ray flares are also observed at even later timescales of ∼few days Fong et al. 2014). The power output of the central engine required to explain this emission is uncertain by several orders of magnitude because it depends on the radiative efficiency and beaming fraction of the (likely jetted) X-ray emission. Nevertheless, comparison of the left and right panels of Fig. 3 makes clear that the energy input of a central remnant can compete with, or even dominate, that of radioactivity on timescales from minutes to weeks after the merger (e.g., Kisaka et al. 2017).

Fall-back accretion
In addition to the unbound ejecta during a NS-NS/BH-NS merger (or the accretion disk after the merger), a comparable mass could remain gravitationally bound to the central remnant. Depending on the energy distribution of this matter, it will fall back to the center and enter the accretion disk over timescales ranging from seconds to days or longer after the coalescence event (Rosswog 2007;Rossi and Begelman 2009;Chawla et al. 2010;Kyutoku et al. 2015). At late times t 0.1 s, the mass fall-back rate decays as a power-laẇ where the normalizationṀ fb (t = 0.1) at the reference time t = 0.1 s can vary from ∼ 10 −3 M s −1 in NS-NS mergers, to values up to an order of magnitude larger in BH-NS mergers (Rosswog 2007;Foucart et al. 2015).
There are several caveats to the presence of fall-back accretion. Simulations show that disk outflows from the inner accretion flow in BH-NS mergers can stifle the fall-back of material, preventing it from reaching the BH on timescales t 100 ms (Fernández et al. 2015b). Heating due to the r -process over the first ∼ 1 second can also unbind matter that was originally marginally-bound, generating a cut-off in the fall-back rate after a timescale of seconds or minutes (Metzger et al. 2010a;Desai et al. 2019). Furthermore, matter which does return to the central remnant is only tenuously bound and unable to cool through neutrinos, which may drastically reduce the accretion efficiency (the fraction ofṀ fb that remains bound in the disk; Rossi and Begelman 2009). Despite these concerns, some fall-back and accretion by the central remnant is likely over the days-weeks timescales of the observed kilonovae.
If matter reaches the central compact object at the rateṀ fb (Eq. 33), then a fraction of the resulting accretion power L acc ∝Ṁ fb c 2 would be available to heat the ejecta, contributing to the kilonova luminosity. The accretion flow is still highly super-Eddington throughout this epoch (L acc L Edd ∼ 10 39 erg s −1 ) and might be expected power a collimated ultra-relativistic jet, similar but weaker than that responsible for generating the earlier GRB. At early times, the jet has sufficient power to propagate through the ejecta, producing high energy emission at larger radii (e.g., powering the short GRB or temporally-extended X-ray emission following the burst). However, as the jet power decreases in time it is more likely to become unstable (e.g., Bromberg and Tchekhovskoy 2016), in which case its Poynting flux or bulk kinetic energy would be deposited as heat behind the ejecta. A mildly-relativistic wind could be driven from the inner fall-back-fed accretion disk, which would emerge into the surroundings and collide/shock against the (potentially slower, but higher mass) ejecta shell, thermal-  Fig. 9, shown separately assuming opacities appropriate to lanthanide-bearing (κ = 20 cm 2 g −1 ; left panel) and lanthanide-free (κ = 1 cm 2 g −1 ; right panel) ejecta. We adopt ejecta heating rate following Eq. (34) for a constant efficiency j = 0.1 and have normalized the fall-back mass to an optimistic valueṀ fb (t = 0.1) = 10 −2 M s −1 izing the wind's kinetic energy and providing a heat source behind the ejecta (Dexter and Kasen 2013).
Heating by fall-back accretion can be crudely parametrized as follows, where j is a jet/disk wind efficiency factor. 15 For optimistic, but not physically unreasonable, values of j ∼ 0.01 − 0.1 andṀ fb (0.1s) ∼ 10 −3 M s −1 , Fig. 3 shows thaṫ Q fb can be comparable to radioactive heating on timescales of days to weeks. Figure 15 shows toy model light curves calculated assuming the ejecta (mass M = 10 −2 M and velocity v 0 = 0.1 c) is heated by fall-back accretion according to Eq. (34) under the optimistic assumption that jṀfb (t = 0.1) = 10 −3 M s −1 on the very high end of the values suggested by merger simulations (Rosswog 2007) and expected for BH-powered outflows ( j ∼ 1). The left and right panels show the results separately in the case of ejecta with a low (κ = 1 cm 2 g −1 ) and high (κ = 20 cm 2 g −1 ) opacity, respectively. Fall-back enhance the peak brightness, particularly those bands which peak during the first 1 day, by up to a magnitude or more, compared to the otherwise similar case with pure r -process heating (reproduced from Fig. 9 and shown for comparison with dashed lines). If the amount of fall-back, and the jet/accretion disk wind efficiency is high, we conclude that accretion power could in principle provide a moderate boost to the observed kilonova emission, particularly in cases where the ejecta mass (and thus intrinsic r -process decay power) is particularly low.
Based on the high observed X-ray luminosity from GRB 130603B simultaneous with the excess NIR emission, Kisaka et al. (2016) argue that the latter was powered by reprocessed X-ray emission rather than radioactive heating (Tanvir et al. 2013;Berger et al. 2013). However, the validity of using the observed X-rays as a proxy for the ejecta heating rely on the assumption that the X-ray emission is instrinsically isotropic (i.e., we are peering through a hole in the ejecta shell), as opposed to being geometrically or relativistically-beamed as a part of a jet-like outflow from the central engine (a substantial beaming-correction to the observed isotropic X-ray luminosity would render it too low to power the observed NIR emission). Matsumoto et al. (2018) and Li et al. (2018) made a similar argument that AT2017gfo was powered by a central engine. However, unlike in 130603B, no X-ray emission in excess of the afterglow from the external shocked ISM was observed following GW170817 at the time of the kilonova . Furthermore, the observed bolometric light curve is well explained by r -process radioactive decay without the need for an additional central energy source (Fig. 10).

Long-lived magnetar remnants
As described in Sect. 3.1, the type of compact remnant produced by a NS-NS merger prompt BH formation, hypermassive NS, supramassive NS, or indefinitely stable NS) depends sensitively on the total mass of the binary relative to the poorly constrained TOV mass, M TOV . A lower limit of M TOV 2-2.1 M is set by measured pulsar masses Antoniadis et al. 2013;Cromartie et al. 2019), while an upper limit of M TOV 2.16 M is suggested for GW170817 (Sect. 5.2). Taken as granted, these limits, combined with the assumption that the measured mass distribution of the Galactic population of binary neutron stars is representative of those in the universe as a whole, leads to the inference that ≈ 18−65% of mergers will result in a long-lived SMNS  instead of the short-lived HMNS most believe formed in GW170817 (Table 3).
The massive NS remnant created by a NS-NS merger will in general have more than sufficient angular momentum to be rotating near break-up (Radice et al. 2018a;however, see Shibata et al. 2019). A NS of mass M ns rotating near its mass-shedding limit possesses a rotational energy, where P = 2π/ is the rotational period and I is the NS moment of inertia, which we have normalized to an approximate value for a relatively wide class of nuclear equations of state I LS ≈ 1.3 × 10 45 (M ns /1.4 M ) 3/2 g cm 2 , motivated by Fig. 1 of Lattimer and Schutz (2005). This energy reservoir is enormous, both compared to the kinetic energy of the merger ejecta (≈ 10 50 −10 51 erg) and to that released by its radioactive decay. Even if only a modest fraction of E rot were to be extracted from the remnant hours to years after the merger by its electromagnetic spin-down, this would substantially enhance the EM luminosity of the merger counterparts (Yu et al. 2013;Gao et al. 2013;Metzger and Piro 2014;Gao et al. 2015;Siegel and Ciolfi 2016a).
This brings us back to a crucial qualitative difference between the formation of a HMNS and a SMNS or stable NS remnant. A HMNS can be brought to the point of collapse by the accretion of mass and redistribution of its internal angular momentum. However, energy dissipated by removing internal differential rotational support can largely be released as heat and thus will escape as neutrino emission (effectively unobservable at typical merger distances). Thus, the angular momentum of the binary can largely be trapped in the spin of the newly-formed BH upon its collapse, rendering most of E rot unavailable to power EM emission.
By contrast, a SMNS is (by definition) supported by its solid-body rotation, even once all forms of differential rotation have been removed. Thus, angular momentum must physically be removed from the system to allow the collapse, and the removal of angular momentum brings with it a concomitant amount of rotational energy. The left panel of Fig. 4 shows the "extractable" rotational energy for mergers which leave SMNS remnants and how quickly it grows with decreasing binary chirp mass (proxy for the mass M tot ). 16 This energy budget increases from 10 51 erg for remnants near the HMNS-SMNS boundary at M tot ≈ 1.2M TOV to the full rotational energy E rot ≈ 10 53 erg (Eq. 35) for the lowest mass, indefinitely stable remnants M tot M TOV .
A strong magnetic field provides an agent for extracting rotational energy from the NS remnant via electromagnetic spin-down (the same mechanism at work in ordinary pulsars). MHD simulations show that the original magnetic fields of the NS in a NS-NS merger are amplified to very large values, similar or exceeding the field strengths of 10 15 −10 16 G of Galactic magnetars (Price and Rosswog 2006;Zrake and MacFadyen 2013;Kiuchi et al. 2014). However, most of this amplification occurs on small spatial scales, and at early times when the NS is still differentially-rotating, resulting in a complex and time-dependent field geometry (Siegel et al. 2014). Nevertheless, by the time the NS enters into a state of solid body rotation (typically within hundreds of milliseconds following the merger), there are reasons to believe the remnant could possess an ordered dipole magnetic field of comparable strength, B ∼ 10 15 −10 16 G. For instance, an ordered magnetic field can be generated by an α − dynamo, driven by the combined action of its rapid millisecond rotation period and thermal and lepton gradient-driven convection in the cooling remnant (Thompson and Duncan 1993).
The spin-down luminosity of an aligned dipole 17 rotator is given by (e.g., Philippov et al. 2015) L sd = μ 2 4 c 3 = 7 × 10 50 erg s −1 I 16 The extractable energy is defined as the difference between the rotational energy at break-up and the rotational energy after the object has spun down to the point of becoming unstable and collapsing into a BH. 17 Unlike vacuum dipole spin-down, the spin-down rate is not zero for an aligned rotator in the force-free case, which is of greatest relevance to the plasma-dense, post-merger environment.
where μ = B R 3 ns is the dipole moment, R ns = 12 km is the NS radius, B is the surface equatorial dipole field, is the characteristic spin-down time over which an order unity fraction of the rotational energy is removed, where P 0 is the initial spin-period and we have assumed a remnant mass of M = 2.3 M . The latter is typically close to, or slightly exceeding, the mass-shedding limit of P = 0.7 ms. 18 The spin-down luminosity in Eq. (36) goes to zero 19 when the NS collapses to the BH at time t collapse . For a stable remnant, t collapse → ∞, but for supramassive remnants, the NS will collapse to a black hole after a finite time which can be estimated 20 by equating t collapse 0 L sd dt to the extractable rotational energy. Yu et al. (2013) suggested 21 that magnetic spin-down power, injected by the magnetar behind the merger ejecta over a timescale of days, could enhance the kilonova emission (the termed such events "merger-novae"; see also Gao et al. 2015). Their model was motivated by similar ideas applied to super-luminous supernovae (Kasen and Bildsten 2010;Woosley 2010; and is similar in spirit to the 'fall-back powered' emission described in Sect. 6.2.1. Although the spin-down luminosity implied by Eq. (36) is substantial on timescales of hours to days, the fraction of this energy which will actually be thermalized by the ejecta, and hence available to power kilonova emission, may be much smaller.
As in the Crab Nebula, pulsar winds inject a relativistic wind of electron/positron pairs. This wind is generally assumed to undergo shock dissipation or magnetic reconnection near or outside a termination shock, inflating a nebula of relativistic particles (Kennel and Coroniti 1984). Given the high energy densities of the post-NS-NS merger environment, electron/positron pairs heated at the shock cool extremely rapidly via synchrotron and inverse Compton emission inside the nebula Siegel and Ciolfi 2016a, b), producing broadband radiation from the radio to gammarays (again similar to conventional pulsar wind nebulae; e.g., Gaensler and Slane, 2006). 22 A fraction of this non-thermal radiation, in particular that at UV and soft X-ray frequencies, will be absorbed by the neutral ejecta walls and reprocessed to 18 If the remnant is born with a shorter period, mass shedding or non-axisymmetric instabilities set in which will result in much more rapid loss of angular momentum to GWs , until the NS rotates at a rate close to P 0 0.7 ms. 19 The collapse event itself has been speculated to produce a brief (sub-millisecond) electromagnetic flare (Palenzuela et al. 2013) or a fast radio burst (Falcke and Rezzolla 2014) from the detaching magnetosphere; however, no accretion disk, and hence long-lived transient, is likely to be produced (Margalit et al. 2015). 20 We also assume that dipole spin-down exceeds gravitational wave losses, as is likely valid if the nonaxisymmetric components of the interior field are 100 times weaker than the external dipole field (Dall'Osso et al. 2009). We have also neglected angular momentum losses due to f -mode instabilities (Doneva et al. 2015). 21 In fact, Kulkarni (2005) earlier had suggested energy input from a central pulsar as a power source. 22 Such magnetar nebulae may also be promising sources of high-energy neutrinos for hours to days after the merger (e.g. Fang and Metzger 2017). lower, optical/IR frequencies , where the lower opacity allows the energy to escape, powering luminous kilonova-like emission.
On the other hand, this non-thermal nebular radiation may also escape directly from the ejecta without being thermalized, e.g., through spectral windows in the opacity. This can occur for hard X-ray energies above the bound-free opacity or (within days or less) for high energy MeV gamma-rays between the decreasing Klein-Nishina cross section and the rising photo-nuclear and γ −γ opacities (Fig. 8). Furthermore, if the engine is very luminous and the ejecta mass sufficiently low, the engine can photo-ionize the ejecta shell, allowing radiation to freely escape even from the far UV and softer X-ray bands (where bound-free opacity normally dominates). While such leakage from the nebula provides a potential isotropic high energy counterpart to the merger at X-ray wavelengths (Metzger and Piro 2014;Siegel and Ciolfi 2016a, b;Wang et al. 2016), it also reduces the fraction of the magnetar spin-down luminosity which is thermalized and available to power optical-band radiation.
We parameterize the magnetar contribution to the ejecta heating aṡ where, as in the fall-back case (Eq. 34), th is the thermalization efficiency. We expect th ∼ 1 at early times when the ejecta is opaque (unless significant energy escapes in a jet), but the value of th will decrease as the optical depth of the expanding ejecta decreases, especially if the ejecta becomes ionized by the central engine. Metzger and Piro (2014) point out another inefficiency, which, unlike radiation leakage, is most severe at early times. High energy MeV gamma-rays in the nebula behind the ejecta produce copious electron/positron pairs when the compactness is high. These pairs in turn are created with enough energy to Compton upscatter additional seed photons to sufficient energies to produce another generation of pairs (and so on. . .). For high compactness 1, this process repeats multiple times, resulting in a 'pair cascade' which acts to transform a significant fraction Y ∼ 0.01−0.1 of the pulsar spin-down power L sd into the rest mass of electron/positron pairs (Svensson 1987;Lightman et al. 1987). Crucially, in order for non-thermal radiation from the central nebula to reach the ejecta and thermalize, it must diffuse radially through this pair cloud, during which time it experiences adiabatic PdV losses. Because at early times the Thomson optical depth of the pair cloud, τ n , actually exceeds the optical depth through the ejecta itself, this suppresses the fraction of the magnetar spin-down power which is available to thermalize and power the emission.
Following Metzger and Piro (2014) and Kasen et al. (2016), we account in an approximate manner for the effect of the pair cloud by suppressing the observed luminosity according to, where L is the luminosity of the kilonova, calculated as usual from the energy equation (14) using the magnetar heat source (Eq. 38), and is the ratio of the characteristic 'lifetime' of a non-thermal photon in the nebula, t life , to the ejecta expansion timescale ∼ t, where A is the (frequency-averaged) albedo of the ejecta. In what follows we assume A = 0.5.
For high spin-down power and early times (t life t), pair trapping acts to reduce the thermalization efficiency of nebular photons, reducing the effective luminosity of the magnetar-powered kilonova by several orders of magnitude compared to its value were this effect neglected. The bottom panel of Fig. 3 shows the spin-down luminosity L sd for stable magnetars with P 0 = 0.7 ms and B = 10 15 , 10 16 G. We also show the spin-down power, 'corrected' by the factor (1+t life /t) −1 , as in Eq. (39) for Y = 0.1. 23 Figures 16 and 17 shows kilonova light curves, calculated from our toy model, but including additional heating of the ejecta due to rotational energy input from an indefinitely stable magnetar with assumed dipole field strengths of B = 10 15 G and 10 16 G, respectively. We again show separately cases in which we assume a high value for the ejecta opacity κ = 20 cm 2 g −1 appropriate for lanthanide-rich ejecta compared to low opacity κ = 1 cm 2 g −1 appropriate to lanthanide-free (light r -process) elements. The latter case is probably the most physical one, because the majority of the disk wind ejecta (which typically dominates the total) is expected to possess a high-Y e in the presence of a long-lived stable merger remnant (t life ≈ ∞ in Fig. 6).
A long-lived magnetar engine has three main effects on the light curve relative to the normal (pure r -process-powered) case: (1) increase in the peak luminosity, by up 1 Page 58 of 89 B. D. Metzger   Fig. 17 Same as Fig. 16, but calculated for an ejecta opacity κ = 1 cm 2 g −1 relevant to lanthanide-free matter to 4−5 magnitudes; (2) more rapid evolution, i.e. an earlier time of peak light; (3) substantially bluer colors. Feature (1) is simply the result of additional heating from magnetar spin-down, while feature (2) results from the greater ejecta velocity due to the kinetic energy added to the ejecta by the portion of the spin-down energy released before the ejecta has become transparent (that portion going into PdV rather than escaping as radiation). Feature (3) is a simple result of the fact that the much higher luminosity of the transient increases its effective temperature, for an otherwise similar photosphere radius near peak light. Figures 16 and 17 represent close to the most "optimistic" effects a long-lived magnetar could have on the kilonova light curves. The effects will be more subtle, and closer to the radioactivity-only models, if the magnetar is a less-stable SMNS (that collapses into a BH before all of E rot is released) or if a substantial fraction of its rotational energy escapes as gamma-rays or is radiated by the magnetar through GW emission (e.g., Dall'Osso and Stella 2007; Corsi and Mészáros 2009), rather than being transferred to the ejecta through its magnetic dipole spin-down. Such effects are easy to incorporate into the toy model by cutting off the magnetar spin-down heating for t ≥ t collapse in Eq. (36), or by including additional losses due to GW radiation into the magnetar spin-down evolution (e.g., Li et al. 2018 and references therein).

Observational prospects and strategies
With the basic theory of kilonova in place, we now discuss several implications for past and present kilonova observations. Table 8 provides rough estimates for the expected range in kilonova luminosities, timescales, isotropy to accompany NS-NS and BH-NS mergers for different assumptions about the merger remnant/outcome. Figure 18 illustrates some of this diversity graphically.

Kilonova candidates following short GRBs
If short duration GRBs originate from NS-NS or NS-BH mergers, then one way to constrain kilonova models is via optical and NIR follow-up observations of nearby week Mostly red N a Estimated range in peak luminosity. Does not account for extra sources of early-time heating from free neutrons or shocks, which could enhance the peak luminosity in the first hours (Sect. 6.1) b Estimated peak timescale c Whether to expect large pole-equatorial (or azimuthal, in the NS-BH case) isotropy in the kilonova properties d Whether a BH-NS merger is accompanied by mass ejection depends on whether the NS is tidally disrupted sufficiently far outside of the BH event horizon to generate unbound tidal material and the formation of an accretion disk. Very roughly, this condition translates into a comparison between the tidal radius of the NS, R t (which depends on the BH-NS mass ratio and the NS radius), and the radius of the innermost stable circular orbit of the BH, R isco (which depends on the BH mass and spin); see futher discussion in Sect. 3.1

Fig. 18
Schematic illustration mapping different types of mergers and their outcomes to trends in their kilonova light curves. The top panel shows the progenitor system, either an NS-NS or an NS-BH binary, while the middle plane shows the final merger remnant (from left to right: an HMNS that collapses to a BH after time t collapse , a spinning magnetized NS, a non-spinning BH and a rapidly spinning BH). The bottom panel illustrates the relative amount of UV/blue emission from an neutron precursor (purple), optical emission from lanthanide-free material (blue) and IR emission from lanthanide containing ejecta (red). Note: the case of a NS-NS merger leading to a slowly spinning black hole is unlikely, given that at a minimum the remnant will acquire the angular momentum of the original binary orbit. Modified from a figure originally presented in Kasen et al. (2015), copyright by the authors short bursts on timescales of hours to a week. All else being equal, the closest GRBs provide the most stringent constraints; however, the non-thermal afterglow emissionthe strength of which can vary from burst to burst-must also be relatively weak, so that it does not outshine the thermal kilonova. Blue kilonova emission similar in luminosity to GW170817 would have been outshone by the afterglow emission in all but a small handful of observed short GRBs, but the longer-lived red kilonova emission has a better chance of sticking out above the fading, relatively blue afterglow. The NIR excess observed following GRB 130603B Tanvir et al. 2013), if powered by the radioactive decay of r -process nuclei, required a total ejecta mass of lanthanide-bearing matter of 0.1 M (Barnes et al. 2016). This is ∼ 3-5 times greater than the ejecta mass inferred for GW170817 Sect. 5). As with GW170817, the ejecta implied by kilonova models of 130603B is too high to be explained by the dynamical ejecta from a NS-NS merger, possibly implicating a BH-NS merger in which the NS was tidally disrupted well outside the BH horizon (Hotokezaka et al. 2013b;Tanaka et al. 2014;Kawaguchi et al. 2016). However, NS-NS mergers can also produce such a high ejecta mass if a large fraction of the remnant accretion disk (which possess masses up to ∼ 0.2 M ) is unbound in disk winds (Siegel and Metzger 2017 found that ≈ 40% of the disk mass could be unbound). Alternatively, the unexpectedly high luminosity of this event could attributed to energy input from a central engine rather than radioactivity (Kisaka et al. 2016), which for fallback accretion indeed produces the correct luminosity to within an order-of-magnitude (Fig. 15). Yang et al. (2015), Jin et al. (2015), Jin et al. (2016) found evidence for NIR emission in excess of the expected afterglow following the short GRBs 050709 and 060614, indicative of possible kilonova emission. The short GRB 080503 ) showed an optical peak on a timescale of ∼ 1 day, which could be explainded as blue kilonova powered by r -process heating (Metzger and Fernández 2014;Kasen et al. 2015) or a central engine (Metzger and Piro 2014;Gao et al. 2015). These possibilities unfortunately could not be distinguished because the host galaxy of GRB080503 was not identified, resulting in its distance and thus luminosity being unconstrained. 24 Gompertz et al. (2018) found three short bursts (GRBs 050509b, 061201, and 080905A) where, if the reported redshifts were correct, deep upper limits rule out the presence of a kilonova similar to AT2017gfo by several magnitudes (see also Fong et al. 2017). Given the diverse outcomes of NS-NS mergers, and how the properties of the remnant can dramatically effect the quantity and composition of the kilonova outflows, variation in the ejecta properties by an order of magnitude or more would not be unexpected. For instance, high-mass mergers that undergo a prompt collapse to a BH eject substantially less mass (particularly of the high-Y e kind capable of producing blue kilonova emission), but still could generate an accretion disk of sufficient mass to power a GRB jet. 25 24 Rebrightening in the X-ray luminosity, coincident with the optical brightening, was also observed following GRB 080503 ). Whether the optical emission is powered exclusively by r -process heating or not, this could potentially be consistent with non-thermal emission from a central engine (Metzger and Piro 2014;Gao et al. 2015;Siegel and Ciolfi 2016a, b). 25 If the jet is powered by the Blandford-Znajek mechanism, the disk mass M t only need be sufficiently massive to hold a magnetic field of the requisite strength in place. This is not a very stringent constraint: Even with deep observations of a particularly nearby burst, Fong et al. (2016a) emphasize the challenges to constrain vanilla blue/red kilonova models with ground follow-up of GRBs. This highlights the crucial role played by the Hubble Space Telescope, and in the future by the James Webb Space Telescope (JSWT) and Wide Field Infrared Survey Telescope (WFIRST), in such efforts. Fortunately, NS-NS mergers detected by Advanced LIGO at distances < 200 Mpc (redshift z < 0.045) are at least three times closer (> 2.5 mags brighter) than the nearest cosmological short GRBs. Nevertheless, the detection rate of well-localized short GRBs is currently higher than GW events. For this reason among others (e.g., the information obtained on GRB jet opening angles from afterglow jet breaks), we advocate for continued space-based observations performing late-time follow-up of short GRB afterglows to search for kilonova signatures.

Gravitational-wave follow-up
This review has hopefully made clear that kilonovae are not likely to be homogeneous in their properties, with potentially significant differences expected in their colors and luminosities, depending on the type of merging system, the properties of the ingoing binary, and, potentially, our viewing angle relative to the binary inclination (see, Figs. 4, 18 and Table 8). Here we discuss prospects and strategies of GW follow-up in the case of different merger outcomes. Prompt collapse to BH In a NS-NS merger the observed emission is expected to depend sensitively on the lifetime of the central NS remnant, which in turn will depend on the in-going binary mass (Table 3). When BH formation is prompt, the ejecta mass will in most cases be low 10 −2 M and radioactivity (and, potentially, fall-back accretion of the tidal tail) will provide the only heating sources. The lack of a HMNS remnant will also result in a greater fraction of ejecta being lanthanide-rich and thus generating red kilonova emission. While some blue ejecta could still originate from the accretion disks, its relatively low velocity ∼ 0.1 c could result in its emission being blocked for equatorial viewing angles. Little or no neutron precursor emission is expected. For a merger at ∼ 100 Mpc, a purely red r -process powered kilonova of ejecta mass ∼ 3 × 10 −3 M would peak over a timescale of a few days in the NIR at I J K ∼ 23-24 (scaling from Fig. 9, right panel). Given such relatively dim emission, only the largest aperture telescopes (e.g., DECam, Subaru HSC, or LSST) are capable of detecting the low-M ej red kilonova of a prompt collapse (Fig. 9, right panel). HMNS remnant (GW170817-like) The situation is more promising for lower-mass mergers in which at least a moderately long-lived HMNS remnant forms. Shock-heated matter from the merger interface can generate high-Y e dynamical ejecta comparable of exceeding that of the lanthanide-rich tidal tail. Likewise, outflows from the magnetized HMNS, or from the accretion disk prior to BH formtion (Fig. 6), will produce a greater quantity of high-Y e lanthanide free material than in the prompt collapse case. For a merger generating ∼ 10 −2 M of high-Y e ejecta (similar to GW170817) at ∼ 100 Footnote 25 continued GRB jets are weak in energy E GRB ∼ 10 49 erg compared to the maximum accretion power available, ∼ M t c 2 ∼ 10 52 (M t /10 −2 M ) erg.
Mpc, the resulting blue kilonova emission could peak at UVR ∼ 19-20 (Fig. 9, right panel) on a timescale of several hours to days. Even if the blue component is somehow not present or is blocked by lanthanide-rich matter, a source at 100 Mpc could still reach U ∼ 20 on a timescale of hours if the outer layers of the ejecta contain free neutrons (Fig. 14) or if the ejecta has been shock heated within a second after being first ejected (Fig. 12). Although not much brighter in magnitude than the NIR peak at later times, the blue kilonova may be the most promising counterpart for the majority of follow-up telescopes, for which the greatest sensitivity at optical wavelengths (a fact that the discovery of AT2017gfo has now made obvious). It is thus essential that follow-up observations begin within hours to one day following the GW trigger. SMNS/stable remnant Although potentially rare, the brightest counterparts may arise from the mergers of low-mass NS-NS binaries which generate long-lived SMNS or stable NS remnants (t collapse 300 ms). Even ignoring the possibility of additional energy input from magnetar spin-down, the quantity of blue disk wind ejecta in this case is substantially enhanced (t collapse = ∞ in Fig. 6) and could approach ∼ 0.1 M , boosting the peak luminosity of the blue kilonova by a magnitude or more from what was observed in GW170817. Allowing also for energy input from the magnetar rotational energy (in addition to radioactivity), the transient at 100 Mpc could reach U V I ≈ 16 − 17 (Figs. 16, 17). However, the precise luminosity is highly uncertain, as it depends on several unknown factors: the dipole magnetic field of the magnetar remnant, the thermalization efficiency of the magnetar nebula by the ejecta (Eq. 38), and the NS collapse time (which in turn depends on the binary mass and the magnetic field strength). Shallower follow-up observations, such as those used in the discovery and follow-up of GW170817, are thus still relevant to kilonova follow-up even for more distant events (they could also be sufficient to detect the on-axis GRB afterglow in the rare case of a face-on merger). BH-NS mergers BH-NS mergers are "all or nothing" events. If the NS is swallowed whole prior to being tidally disrupted, then little or no kilonova emission is expected. However, in the potentially rare cases when the BH is low in mass and rapidly spinning (in the prograde orbital sense), then the NS is tidally disrupted well outside of the horizon and the quantity of dynamical ejecta can be larger than in NS-NS mergers, by a typical factor of ∼ 10 (Sect. 3.1). All else being equal, this results in the kilonova peaking one magnitude brighter in BH-NS mergers. Likewise, the mass fall-back rate in BH-NS mergers can be up to ∼ 10 times higher than in NS-NS mergers (Rosswog 2007), enhancing potential accretion-powered contributions to the kilonova emission (Fig. 15, bottom panel).
However, the amount of high-Y e ejecta is potentially less than in NS-NS mergers due to the lack in BH-NS mergers of shock-heated ejecta or a magnetar remnant, and for the same reason no neutron precursor is anticipated, unless it can be somehow generated by the GRB jet. The accretion disk outflows could still produce a small quantity of blue ejecta, but its velocity is likely to be sufficiently low ∼ 0.1 c that it will be blocked by the (faster, more massive) tidal tail, at least for equatorial viewing angles. Taken together, the kilonova emission from BH-NS mergers is more likely to dominated by the red component, although moderate amounts of high-Y e matter and blue emission could still be produced by the disk winds (Just et al. 2015;Fernández et al. 2015b). Unfortunately for purposes of follow-up, any benefits of the higher dynamical ejecta mass on the light curve luminosity may be more than offset by the larger expected source distance, which will typically be ≈ 2-3 times greater than the 200 Mpc horizon characteristic of NS-NS mergers for an otherwise equal GW event detection rates. See, e.g., Bhattacharya et al. (2019) for further discussion of the diverse EM counterparts of BH-NS mergers. Search strategies Several works have explored the optimal EM follow-up strategies of GW sources, or ways to achieve lower latency GW triggers (Metzger and Berger 2012;Cowperthwaite and Berger 2015;Gehrels et al. 2016;Ghosh et al. 2016;Howell et al. 2016;Rana et al. 2017). Extremely low latency (Cannon et al. 2012;Chen and Holz 2017), though crucial to searching for a potential low-frequency radio burst (Kaplan et al. 2016), is generally not essential for kilonova follow-up. One possible exception is the speculative neutron precursor (Sect. 6.1.1), which peaks hours after the merger. However, in this case, the greatest advantage is arguably to instead locate the follow-up telescope in North America (Kasliwal and Nissanke 2014), producing a better chance of the source occuring near zenith (since the LIGO detectors are most sensitive directly above or below the detectors).
Future EM follow-up efforts would also be aided if LIGO were to provide more information on the properties of its binaries to the wider astronomy community at the time of the GW trigger . The predicted signal from a NS-NS is expected to depend on the binary inclination and the total binary mass, M tot . The inclination cannot be measured with high precision because it is largely degenerate with the (initially unknown) source distance, though it could be determined once the host galaxy was identified. However, the chirp mass and total binary mass can be reasonably accurately determined in low latency (e.g., Biscoveanu et al. 2019). Once the mapping between EM counterparts and the binary mass is better established, providing the binary mass in low latency could provide crucial information for informing search strategies, or even prioritizing, EM follow-up. This is especially important given the cost/scarcity of follow-up resources capable of performing these challenging deep searches over large sky areas.
The generally greater sensitivity of telescopes at optical wavelengths, as compared to the infrared, motivates a general strategy by which candidate targets are first identified by wide-field optical telescopes on a timescale of days, and then followed-up with spectroscopy or photometry in the NIR over a longer timescale of ∼ 1 week. Cowperthwaite and  show that no other known or predicted astrophysical transients are as red and evolve as quickly as kilonovae, thus reducing the number of optical false positives to a manageable level. Follow-up observations of candidates at wavelengths of a few microns could be accomplished, for instance, by the James Webb Space Telescope (Bartos et al. 2016), WFIRST (Gehrels et al. 2015), or a dedicated GW follow-up telescope with better target-of-opportunity capabilities.
Another goal for future kilonova observations would be a spectroscopic measurement of absorption lines from individual r -process elements (however, see Watson et al. 2019, for a possible detection of Sr ii in GW170817). Individual lines are challenging to identify for the simple reason that most of their specific wavelengths cannot be predicted theoretically with sufficient precision and have not been measured experimentally. Furthermore, at early times, the absorption lines are Doppler-broadened near peak light due to the substantial velocities v 0.1 c of the ejecta. Broad absorp-tion features were seen in AT2017gfo (e.g., Chornock et al. 2017), but these likely represented a blend of multiple lines. Fortunately, line-widths become narrower postmaximum as the photosphere recedes to lower velocity coordinates through the ejecta and nebular lines appear (ejecta velocities as low as ∼ 0.03 c are predicted for the disk wind ejecta). Unfortunately, emission becomes significantly dimmer at these late times and line blending could remain an issue. Spectroscopic IR observations of such dim targets is a compelling science case for future 30-meter telescopes. For instance, the planned Infrared Imaging Spectrograph (IRIS) on the Thirty Meter Telescope (Skidmore et al. 2015) will obtain a signal to noise ratio of 10 per wavelength channel (spectral resolution R = 4000) for a K = 25 mag point source.

Summary of predictions
We conclude with a summary of predictions for a future large samples of kilonovae. Once Advanced LIGO/Virgo reach design sensitivity within a few years, NS-NS mergers could be detected as frequently as once per week to once per month. This rate will increase by roughly another factor of ∼ 8 with the planned LIGO A+ upgrades in the mid 2020s (Reitze et al. 2019b), and yet further with proposed third generation GW detectors, such as Einstein Telescope (Punturo et al. 2010) and Cosmic Explorer (Reitze et al. 2019a), which could come online in the 2030s.
As mentioned earlier in this review, high-fidelity first-principles kilonova models are not currently available. Reasons for this include uncertainties in: (1) the predicted ejecta properties (mass, velocity, Y e ) due to numerical limitations combined with present ignorance about the NS equation of state; (2) the properties of unstable neutronrich nuclei, which determine the details of the r -process and the radioactive heating rate; (3) radiative transfer in face of the complex atomic structure of heavy r -process elements and the potential break-down of the approximations which are more safely applied to modeling supernova with less exotic ejecta composition.
Nevertheless, we can still highlight a few trends which a future sample of joint GW/EM events should bear out if our understanding of these events is even qualitatively correct. These predictions include: -For a total binary mass (∼ chirp mass) significantly higher than GW170817, the remnants of NS-NS merger will undergo a prompt collapse to a BH, resulting in a kilonova which is dimmer and redder than AT2017gfo. The blue component of the ejecta, if present at all, will arise from the relatively low-velocity disk wind and is likely be blocked for equatorial viewers by the lanthanide-rich tidal tail ejecta. A GRB jet can still be produced because a low-mass accretion disk can still form, but it could be less energetic than the off-axis jet inferred for GW170817. Such events could be relatively rare if the prompt collapse threshold mass is high. -For binary masses below the (uncertain) prompt collapse threshold mass, the merger will form a HMNS remnant, perhaps similar to that believed to be generated in GW170817. The strength of the blue kilonova relative to the red kilonova emission will increase with decreasing total binary mass, as the ejecta is dominated by the disk wind and the average Y e rises with the collapse time (Fig. 6). BH formation is likely to still be relatively prompt, providing no hindrance to the production of an ultra-relativistic GRB jet and afterglow emission.
-Below another uncertain critical binary mass threshold, the merger remnant will survive for minutes or longer as a quasi-stable SMNS or indefinitely stable NS remnant. The disk wind ejecta will be almost entirely blue but a weak red component could still be present from the tidal tail ejecta (even though such an emission component was likely swamped in GW170817). The mean velocity of the ejecta will increase with decreasing remnant mass, as the amount of rotational energy extracted from the remnant prior to BH formation increases (Fig. 4). Less certain is whether a powerful ultra-relativistic GRB jet will form in this case, due to baryon pollution from the wind of the remnant (though exceptions to such an expectation would be highly informative). However, the afterglow produced by the interaction of the jet/kilonova ejecta with the ISM (particularly the radio emission, which is roughly isotropic and generated by even trans-relativistic ejecta) will be even more luminous than in the HMNS case due to the addition of magnetar rotational energy (Metzger and Bower 2014;Horesh et al. 2016;Fong et al. 2016b). -For BH-NS mergers, in the (possibly rare) subset of low-mass/high-spin BHs which disrupt the NS well outside the horizon, the red kilonova will be more luminous, and extend to higher velocities, than in the NS-NS case due to the greater quantity of tidal tail ejecta. As in the prompt collapse of a NS-NS, the blue component-if present-will arise from the relatively low-velocity disk wind and thus could be blocked for equatorial viewers by the lanthanide-rich tidal tail ejecta, which is likely to be more massive than in the NS-NS case.

Final thoughts
The first edition of this review was written in the year prior to the discovery of GW170817. In the final section of that article I pondered: "Given the rapid evolution of this field in recent years, it is natural to question the robustness of current kilonova models. What would it mean if kilonova emission is ruled out following a NS-NS merger, even to stringent limits?" I concluded this was untenable: "First, it should be recognized that-unlike, for instance, a GRB afterglow-kilonovae are largely thermal phenomena. The ejection of neutron-rich matter during a NS-NS merger at about ten percent of the speed of light appears to be a robust consequence of the hydrodynamics of such events, which all modern simulations agree upon. Likewise, the fact that decompressing nuclear-density matter will synthesize heavy neutron rich isotopes is also robust. The properties of individual nuclei well off of the stable valley are not well understood, although that will improve soon due to measurements with the new Facility for Rare Isotope Beams (e.g., Horowitz et al. 2019). However, the combined radioactive heating rate from a large ensemble of decaying nuclei is largely statistical in nature and hence is also relatively robust, even if individual isotopes are not; furthermore, most of the isotopes which contribute to the heating on the timescale of days to weeks most relevant to the kilonova peak are stable enough that their masses and half-lives are experimentally measured." Despite the success that astrophysics theorists and numerical relativists had in anticipating many of the properties of GW170817, the tables may soon be turned as we struggle to catch up to the rich phenomenology which is likely to explode from the new GW/EM science. In particular, much additional work is required on the theory side to transform kilonovae into quantitative probes of the diversity of merger outcomes and nuclear physics. Among the largest remaining uncertainty in kilonova emission relates to the wavelength-dependent opacity of the ejecta, in particular when it includes lanthanide/actinides isotopes with partially-filled f-shell valence shells Tanaka and Hotokezaka 2013;Fontes et al. 2015). As discussed in Sect. 3.2, the wavelengths and strengths of the enormous number of lines of these elements and ionization states are not experimentally measured and are impossible to calculate from first principles from multi-body quantum mechanics with current computational capabilities. Furthermore, how to handle radiative transport in cases when the density of strong lines becomes so large that the usual expansion opacity formalism breaks down requires further consideration and simulation work. Another theoretical issue which deserves prompt attention is the robustness of the presence of free neutrons in the outermost layers of the ejecta, given their potentially large impact on the very early-time kilonova optical emission (Sect. 6.1.1). The issue of large-scale magnetic field generation, and its impact on the GRB jet and post-merger outflows (e.g., from the remnant NS or accretion disk), is likely to remain a challenging issue for man years to come, one which is hardly unique to this field of astrophysics. Here, I suspect that nature will teach us more than we can deduce ourselves.
With ongoing dedicated effort, as more detections or constraints on kilonovae become possible over the next few years, we will be in an excellent position use these observations to probe the physics of binary NS mergers, their remnants, and their role as an origin of the r -process.