Mergers of binary neutron star systems: a multi-messenger revolution

On 17 August 2017, less than two years after the direct detection of gravitational radiation from the merger of two ~30 Msun black holes, a binary neutron star merger was identified as the source of a gravitational wave signal of ~100 s duration that occurred at less than 50 Mpc from Earth. A short GRB was independently identified in the same sky area by the Fermi and INTEGRAL satellites for high energy astrophysics, which turned out to be associated with the gravitational event. Prompt follow-up observations at all wavelengths led first to the detection of an optical and infrared source located in the spheroidal galaxy NGC4993 and, with a delay of ~10 days, to the detection of radio and X-ray signals. This paper revisits these observations and focusses on the early optical/infrared source, which was thermal in nature and powered by the radioactive decay of the unstable isotopes of elements synthesized via rapid neutron capture during the merger and in the phases immediately following it. The far-reaching consequences of this event for cosmic nucleosynthesis and for the history of heavy elements formation in the Universe are also illustrated.


INTRODUCTION
Although the birth of "multi-messenger" astronomy dates back to the detection of the first solar neutrinos in the '60s and was rejuvenated by the report of MeV neutrinos from SN1987A in the Large Magellanic Cloud, the detection of gravitational radiation from the binary neutron star merger on 17 August 2017 (GW170817A) marks the transition to maturity of this approach to observational astrophysics, as it is expected to open an effective window into the study of astrophysical sources that is not limited to exceptionally close (the Sun) or rare (Galactic supernova) events. GW170817 is a textbook case for gravitational physics, because, with its accompanying short gamma-ray burst (GRB) and afterglow, and its thermal aftermath ("kilonova"), it has epitomized the different epiphanies of the coalescence of a binary system of neutron stars, and finally allowed us to unify them.
Owing its name to a typical peak luminosity of ∼ 10 42 erg s −1 , i.e. 1000 times larger than that of a typical nova outburst, kilonova is the characteristic optical and infrared source accompanying a binary neutron star merger and due to the radioactive decay of the many unstable isotopes of large atomic weight elements synthesized via rapid neutron capture in the promptly formed dynamical ejecta and in the delayed post-merger ejecta. Its evolution, as well as that of the GRB afterglow, was recorded with exquisite detail, thanks to its closeness (40 Mpc). Scope of this paper is to review the electromagnetic multi-wavelength observations of GW170817 with particular attention to the kilonova phenomenon.
The outline of the paper is as follows: Section 2 sets the context of binary systems of neutrons stars and describes the predicted outcomes of their coalescences; Section 3 presents the case of GW170817, the only so far confirmed example of double neutron star merger and the multi-wavelength features of its electromagnetic counterpart (short GRB and kilonova); Section 4 focusses on the kilonova, elaborates on its observed optical and near-infrared light curves and spectra, draws the link with nucleosynthesis of heavy elements, and outlines the theoretical framework that is necessary to describe the kilonova properties and implications; Section 5 summarizes the results and provides an outlook of this line of research in the near future.

BINARY NEUTRON STAR MERGERS
Neutron stars are the endpoints of massive stars evolution and therefore ubiquitous in the Universe: on average, they represent about 0.1% of the total stellar content of a galaxy. Since massive stars are mostly in binary systems (Sana et al., 2012), neutron star binaries should form readily, if the supernova explosion of either progenitor massive star does not disrupt the system (Renzo et al., 2019). Alternatively, binary neutron star systems can form dynamically in dense environments like stellar clusters (see Ye et al., 2020, and references therein). Binary systems composed by a neutron star and a black hole are also viable, but rare (Pfahl et al., 2005), which may account for the fact that none has so far been detected in our Galaxy.
The prototype binary neutron star system in our Galaxy is PSR B1913+16, where one member was detected as a pulsar in a radio survey carried out at the Arecibo observatory (Hulse & Taylor, 1974), and the presence of its companion was inferred from the periodic changes in the observed pulsation period of 59 ms (Hulse & Taylor, 1975). Among various tests of strong general relativity enabled by the radio monitoring of this binary system, which earned the Nobel Prize for Physics to the discoverers in 1993, was the measurement of the shrinking of the binary system orbit, signalled by the secular decrease of the 7.75 hours orbital period, that could be entirely attributed to energy loss via gravitational radiation (Taylor & Weisberg, 1982;Weisberg & Huang, 2016, and references therein).
With an orbital decay rate ofṖ = −2.4 × 10 −12 s s −1 , the merging time of the PSR B1913+16 system is ∼300 Myr. Following the detection of PSR B1913+16, another dozen binary neutron stars systems were detected in our Galaxy (e.g. Wolszczan, 1991;Burgay et al., 2003;Tauris et al., 2017;Martinez et al., 2017). Almost half of these have estimated merging times significantly shorter than a Hubble time. The campaigns conducted by the LIGO interferometers in Sep 2015-Jan 2016 (first observing run) and, together with Virgo, in Nov 2016-Aug 2017 (second observing run), the latter leading to the first detection of gravitational waves from a merging double neutron star system (see Section 3), constrained the local merger rate density to be 110-3840 Gpc −3 yr −1 (Abbott et al., 2019). This is consistent with previous estimates (see e.g. Burgay et al., 2003), and, under a series of assumptions, marginally consistent with independent estimates based on double neutron star system formation in the classical binary evolution scenario (Chruslinska et al., 2018). Ye et al. (2020) have estimated that the fraction of merging binary neutron stars that have formed dynamically in globular clusters is negligible. Under the assumption that the event detected by LIGO on 25 April 2019 was produced by a binary neutron star coalescence the local rate of neutron star mergers would be updated to 250 − 2810 Gpc −3 yr −1 (Abbott et al., 2020a).
The merger of a binary neutron star system has three predicted outcomes: 1) a gravitational wave signal that is mildly isotropic, with a stronger intensity in the polar direction than in the equatorial plane; 2) a relativistic outflow, which is highly anisotropic and can produce an observable high energy transient;

Binary neutron star mergers
3) a thermal, radioactive source emitting most of its energy at ultraviolet, optical and near-infrared wavelengths, as detailed in the next three sub-sections.

Gravitational waves
Coalescing binary systems of degenerate stars and stellar mass black holes are optimal candidates for the generation of gravitational waves detectable from ground-based interferometers as the strong gravity conditions lead to huge velocities and energy losses (Shapiro & Teukolsky, 1983), and the frequency of the emitted gravitational waves reaches several kHz, where the sensitivity of the advanced LIGO, Virgo and KAGRA interferometers is designed to be maximal (Abbott et al., 2018).
The time behavior of binary systems of compact stars consists of three phases: a first inspiral phase in a close orbit that shrinks as gravitational radiation of frequency proportional to the orbital frequency is emitted, a merger phase where a remnant compact body is produced as a result of the coalescence of the two stars, and a post-merger, or ringdown, phase where the remnant still emits gravitational radiation while settling to its new stable configuration. During the inspiral, the amplitude of the sinusoidal gravitational signal rapidly increases as the distance between the two bodies decreases and the frequency increases (chirp), while in the ringdown phase the signal is an exponentially damped sinusoid. This final phase may encode critical information on the equation of state of the newly formed remnant (a black hole or, in the case of light neutron stars, a massive neutron star or metastable supramassive neutron star). The mathematical tool that is used to describe this evolution is the waveform model, that aims at reproducing the dynamics of the system through the application of post-Newtonian corrections of increasing order and at providing the essential parameters that can then be compared with the interferometric observations (Blanchet, 2014;Nakano et al., 2019).
Since the amplitude of gravitational waves depends on the masses of the binary member stars, the signal will be louder, and thus detectable from larger distances, for binary systems that involve black holes than those with neutron stars. The current horizon for binary neutron star merger detection with LIGO is ∼200 Mpc, and 25-30% smaller with Virgo and KAGRA (Abbott et al. 2018). The dependence of the gravitational waves amplitude on the physical parameters of the system implies that gravitational wave sources are standard sirens (Schutz, 1986), provided account is taken of the correlation between the luminosity distance and the inclination of the orbital plane with respect to the line of sight (Nissanke et al., 2010;Abbott et al., 2016).

Short gamma-ray bursts
GRBs, flashes of radiation of 100-1000 keV that outshine the entire Universe in this band, have durations between a fraction of a second and hundreds or even thousands of seconds. However, the duration distribution is bimodal, with a peak around 0.2 seconds (short or sub-second GRBs) and one around 20 seconds (long GRBs; Kouveliotou et al., 1993). This bimodality is reflected in the spectral hardness, which is on average larger in short GRBs, and in a physical difference between the two groups. While most long GRBs are associated with core-collapse supernovae (Galama et al., 1998;Woosley & Bloom, 2006;Levan et al., 2016), sub-second GRBs are produced by the merger of two neutron stars or a neutron star and a black hole, as long predicted based on circumstantial evidence (Eichler et al., 1989;Fong et al., 2010;Berger et al., 2013;Tanvir et al., 2013) and then proven by the detection of GW170817 and of its high energy counterpart GRB170817A (Section 3). The observed relative ratio of long versus short GRBs depends on the detector sensitivity and effective energy band (e.g. Burns et al., 2016). However, the Pian Binary neutron star mergers duration overlap of the two populations is very large, so that the minimum of the distribution has to be regarded as a rather vaguely defined value (Bromberg et al., 2013).
About 140 short GRBs were localized so far to a precision that is better than 10 arc-minutes 1 ; of these, ∼100, ∼40 and ∼10 have a detected afterglow in X-rays, optical and radio wavelengths, respectively, and ∼30 have measured redshifts (these range between z = 0.111 and z = 2.211, excluding the nearby GRB170817A, see Section 3.1.1, and GRB090426, z = 2.61, whose identification as a short GRB is not robust, Antonelli et al., 2009). Short GRBs are located at projected distances of a fraction of, to several kiloparsecs from the centers of their host galaxies, which are of both early and late type, reflecting the long time delay between the formation of the short GRB progenitor binary systems and their mergers (Berger, 2014).
According to the classical fireball model, both prompt event and multi-wavelength afterglow of short GRBs are produced in a highly relativistic jet directed at a small angle with respect to the line of sight, whose aperture can be derived from the achromatic steepening (or "jet break") of the observed afterglow light curve (Nakar, 2007). In principle, this could be used to reconstruct the collimation-corrected rate of short GRBs, to be compared with predictions of binary neutron star merger rates. However, these estimates proved to be very uncertain, owing to the difficulty of measuring accurately the jet breaks in short GRB afterglows (Fong et al., 2015;Jin et al., 2018;Lamb et al., 2019;Pandey et al., 2019).

r-process nucleosynthesis
Elements heavier than iron cannot form via stellar nucleosynthesis, as not enough neutrons are available for the formation of nuclei and temperatures are not sufficiently high to overcome the repulsive Coulomb barrier that prevents acquisition of further baryons into nuclei (Burbidge et al., 1957). Supernovae (especially the thermonuclear ones) produce large amounts of iron via decay (through 56 Co) of radioactive 56 Ni synthesized in the explosion. Heavier nuclei form via four neutron capture processes (Thielemann et al., 2011), the dominant ones being slow and rapid neutron capture, in brief s-and rprocess, respectively, where "slow" and "rapid" refer to the timescale of neutron accretion into the nucleus with respect to that of the competing process of β − decay. In the s-process, neutron captures occur with timescales of hundreds to thousands of years, making β − decay highly probable, while r-process neutron capture occurs on a timescale of ∼0.01 seconds, leading to acquisition of many neutrons before β − decay can set on. As a consequence, the s-process produces less unstable, longer-lived isotopes, close to the so-called valley of β-stability (the decay time of a radioactive nucleus correlates inversely with its number of neutrons), while the r-process produces the heaviest, neutron-richest and most unstable isotopes of heavy nuclei, up to uranium (Sneden et al., 2008;Côté et al., 2018;Horowitz et al., 2019;Kajino et al., 2019;Cowan et al., 2020). Among both s-process and r-process elements, some are particularly stable owing to their larger binding energies per nucleon, which causes their abundances to be relatively higher than others. In the abundances distribution in the solar neighbourhood these are seen as maxima ("peaks") centered around atomic numbers Z = 39 (Sr-Y-Zr) , 57 (Ba-La-Ce-Nd), 82 (Pb) for the s-process and, correspondingly somewhat lower atomic numbers Z = 35 (Se-Br-Kr), 53 (Te-I-Xe), 78 (Ir-Pt-Au) for the r-process (e.g. Cowan et al., 2020).
Both s-process and r-process naturally occur in environments that are adequately supplied with large neutron fluxes. For the s-process, these are eminently asymptotic giant branch stars, where neutron captures are driven by the 13 C(α, n) 16 O and 22 Ne(α, n) 25 Mg reactions (Busso et al., 1999). The r-process requires much higher energy and neutron densities, which are only realized in most physically extreme environments. While it can be excluded that big-bang nucleosynthesis can accommodate heavy elements formation in any significant amount (Rauscher et al., 1994), there is currently no consensus on the relative amounts of nucleosynthetic yields in the prime r-process candidate sites: core-collapse supernovae and mergers of binary systems composed by neutron stars or a neutron star and a black hole.
Core-collapse supernovae have been proposed starting many decades ago as sites of r-process nucleosynthesis through various mechanisms and in different parts of the explosion, including dynamical ejecta of prompt explosions of O-Ne-Mg cores (Hillebrandt et al., 1976;Wheeler et al., 1998;Wanajo et al., 2002); C+O layer of O-Ne-Mg-core supernovae (Ning et al., 2007); He-shell exposed to intense neutrino flux (Epstein et al., 1988;Banerjee et al., 2011); re-ejection of fallback material (Fryer et al., 2006); neutrino-driven wind from proto-neutron stars (Woosley et al., 1994;Takahashi et al., 1994); magnetohydrodynamic jets of rare core-collapse SNe (Nishimura et al., 2006;Winteler et al., 2012). Similarly old is the first proposal that the tidal disruption of neutron stars by black holes in close binaries (Lattimer & Schramm, 1974, 1976Symbalisty & Schramm, 1982;Davies et al., 1994) and coalescences of binary neutron star systems (Eichler et al., 1989) could be at the origin of r-process nucleosynthesis. This should manifest as a thermal optical-infrared source of radioactive nature of much lower luminosity (a factor of 1000) and shorter duration (rise time of a few days) than supernova (Li & Paczyński, 1998).
The models for r-process elements production in core-collapse supernova all have problems inherent their physics (mostly related to energy budget and neutron flux density). On the other hand, the binary compact star merger origin may fail to explain observed r-process element abundances in very low metallicities stars, i.e. at very early cosmological epochs, owing to the non-negligible binary evolution times (see Cowan et al., 2020 for an accurate review of all arguments in favour and against either channel). While the event of 17 August 2017 (Section 3) has now provided incontrovertible evidence that binary neutron star mergers host r-process nucleosynthesis, the role of core-collapse supernovae cannot be dismissed although their relative contribution with respect to the binary compact star channel must be assessed (Ji et al., 2016;Ramirez-Ruiz et al., 2015;Côté et al., 2019;Safarzadeh et al., 2019;Simonetti et al., 2019). It cannot be excluded that both "weak" and "strong" r-process nucleosynthesis takes place, with the former occurring mainly in supernova and possibly failing to produce atoms up to the third peak of r-process elemental abundance distribution (Cowan et al., 2020). The hint that heavy elements may be produced in low-rate events with high yields (Sneden et al., 2008;Wallner et al., 2015;Macias & Ramirez-Ruiz, 2019) points to binary compact star mergers or very energetic (i.e. expansion velocities larger than 20000 km s −1 ) core-collapse supernovae as progenitors, rather than regular corecollapse supernovae. Along these lines, Siegel et al. (2019) have proposed that the accretion disks of collapsars (the powerful core-collapse supernovae that accompany long GRBs, Woosley & Bloom, 2006) produce neutron-rich outflows that synthesize heavy r-process nuclei.
Neutrons are tightly packed together in neutrons stars, but during coalescence of a binary neutron star system the tidal forces disrupt them and the released material forms promptly a disk-like rotating structure (dynamical ejecta, Rosswog et al., 1999;Shibata & Hotokezaka, 2019) where the neutron density rapidly drops to optimal values for r-process occurrence (∼ 10 24−32 neutrons cm −3 , Freiburghaus et al., 1999) and for copious formation of neutron-rich stable and unstable isotopes of large atomic number elements (Fernández & Metzger, 2016;Tanaka, 2016;Tanaka et al., 2018;Wollaeger et al., 2018;Metzger, 2019).

THE BINARY NEUTRON STAR MERGER OF 17 AUGUST 2017
On 17 August 2017, the LIGO and Virgo interferometers detected for the first time a gravitational signal that corresponds to the final inspiral and coalescence of a binary neutron star system (Abbott et al., 2017a). The sky uncertainty area associated with the event was 28 square degrees, in principle too large for a uniform search for an electromagnetic counterpart with ground-based and orbiting telescopes. However, its small distance (40 +8 −14 Mpc), estimated via the "standard siren" property of gravitational wave signals associated with binary neutron star mergers, suggested that the aftermath could be rather bright, and motivated a large-scale campaign at all wavelengths from radio to very high energy gamma-rays, which was promptly and largely rewarded by success and then timely followed by a long and intensive monitoring (Abbott et al., 2017b,c), as described in Section 3.1. Searches of MeV-to-EeV neutrinos directionally coincident with the source using data from the Super-Kamiokande, ANTARES, IceCube, and Pierre Auger Observatories between 500 seconds before and 14 days after the merger returned no detections (Albert et al., 2017;Abe et al., 2018).
Based on the very detection of electromagnetic radiation, Bauswein et al. (2017) have argued that the merger remnant may not be a black hole or at least the post-merger collapse to a black hole may be delayed. Since the post-merger phase ("ring-down") signal of GW170817 was not detected (Abbott et al., 2017e), this hypothesis cannot be tested directly with gravitational data. Bauswein et al. (2017) also derived lower limits on the radii of the neutron stars.

The electromagnetic counterpart of GW170817
Independent of LIGO-Virgo detection of the gravitational wave signal, the Gamma-ray Burst Monitor (GBM) onboard the NASA Fermi satellite and the Anticoincidence Shield for the gamma-ray Spectrometer (SPI) of the International Gamma-Ray Astrophysics Laboratory (INTEGRAL) satellite were triggered by a faint short GRB (duration of ∼2 seconds), named GRB170817A (Abbott et al., 2017b;Goldstein et al., 2017;Savchenko et al., 2017). This gamma-ray transient, whose large error box was compatible with that determined by LIGO-Virgo, lags the gravitational merger by 1.7 seconds, a delay that may be dominated by the propagation time of the jet to the gamma-ray production site (Beniamini et al., 2020; see however Salafia et al., 2018). The preliminary estimate of the source distance provided a crucial constraint on the maximum distance of the galaxy that could plausibly have hosted the merger, so that the searching strategy was based on targeting galaxies within a ∼50 Mpc cosmic volume (see e.g. Gehrels et al., 2016) with telescopes equipped with large (i.e. several square degrees) field-of-view cameras.
About 70 ground-based optical telescopes participated to the hunt and each of them adopted a different pointing sequence. This systematic approach enabled many groups to identify the optical counterpart candidate in a timely manner (with optical magnitude V ≃ 17), i.e. within ∼12 hours of the merger (Arcavi et al., 2017;Lipunov et al., 2017;Soares-Santos et al., 2017;Valenti et al., 2017;Tominaga et al., 2018). Coulter et al. (2017) were the first to report a detection with the optical 1m telescope Swope at Las Campanas Observatory. The optical source lies at 10 arc-seconds angular separation, corresponding to a projected distance of ∼2 kpc, from the center of the spheroidal galaxy NGC 4993 at 40 Mpc (Blanchard et al., 2017;Im et al., 2017;Levan et al., 2017;Pan et al., 2017;Tanvir et al., 2017).
Rapid follow-up of the gravitational wave and GRB signal in X-rays did not show any source comparable to, or brighter than a typical afterglow of a short GRB. Since both the gravitational data and the faintness of the prompt GRB emission suggested a jet viewed significantly off-axis, this could be expected, as the afterglows from misaligned GRB jets have longer rise-times than those of jets observed at small viewing angles (Van Eerten & MacFadyen, 2011). Therefore, X-ray monitoring with Swift/XRT, Chandra and Nustar continued, and ∼10 days after merger led to the detection with Chandra of a faint source (L X ≃ 10 40 erg s −1 ) (Evans et al., 2017;Margutti et al., 2017;Troja et al., 2017), whose intensity continued to rise up to ∼100 days (D'Avanzo et al., 2018;Troja et al., 2020). Similarly, observations at cm and mm wavelengths at various arrays, including VLA and ALMA, failed to detect the source before ∼16 days after the gravitational signal, which was interpreted as evidence that a jetted source accompanying the binary neutron star merger must be directed at a significant angle (≥20 degrees) with respect to the line of sight (Alexander et al., 2017;Andreoni et al., 2017;Hallinan et al., 2017;Kim et al., 2017;Pozanenko et al., 2018).
The Fermi Large Area Telescope covered the sky region of GW170817 starting only 20 minutes after the merger, and did not detect any emission in the energy range 0.1-1 GeV to a limiting flux of 4.5 × 10 −10 erg s −1 cm −2 in the interval 1153-2027 seconds after the merger (Ajello et al., 2018). Follow-up observation with the atmospheric Cherenkov experiment H.E.S.S. (0.3-8 TeV) from a few hours to ∼5 days after merger returned no detection to a limit of a few 10 −12 erg s −1 cm −2 (Abdalla et al., 2017). A summary of the results of the multi-wavelength observing campaign within the first month of gravitational wave signal detection are reported in Abbott et al. (2017c).
While the radio and X-ray detections are attributed to the afterglow of the short GRB, the ultraviolet, optical and near-infrared data are dominated by the kilonova at early epochs (with a possible contribution at 4 days at blue wavelengths from cooling of shock-heated material around the neutron star merger, Piro & Kollmeier, 2018), and later on by the afterglow, as described in the next two Sections.

The gamma-ray burst and its multi-wavelength afterglow
The short GRB170817A, with an energy output of ∼ 10 46 erg, was orders of magnitude dimmer than most short GRBs (Berger, 2014). Together with a viewing angle of ∼30 deg estimated from the gravitational wave signal (Abbott et al., 2017a), this led to the hypothesis that the GRB was produced by a relativistic jet viewed at a comparable angle. However, the early light curve of the radio afterglow is not consistent with the behavior predicted for an off-axis collimated jet and rather suggests a quasi-spherical geometry, possibly with two components, a more collimated one and a nearly isotropic and mildly relativistic one, which is responsible also for producing the gamma-rays (Mooley et al., 2018a). This confirms numerous predictions whereby the shocked cloud surrounding a binary neutron star merger forms a mildly relativistic cocoon that carries an energy comparable to that of the jet and is responsible for the prompt emission and the early multi-wavelength afterglow (Lazzati et al., 2017a,b;Nakar & Piran, 2017;Bromberg et al., 2018;Xie et al., 2018), and is supported by detailed numerical simulations (Lazzati et al., 2018;Gottlieb et al., 2018).
Using milli-arcsecond resolution radio VLBI observations at 75 and 230 days Mooley et al. (2018b) detected superluminal motion with β = 3 − 5, while Ghirlanda et al. (2019) determine that, at 207 days, the source is still angularly smaller than 2 milli-arcseconds at the 90% confidence, which excludes that a nearly isotropic, mildly relativistic outflow is responsible for the radio emission, as in this case the source apparent size, after more than six months of expansion, should be significantly larger and resolved by the VLBI observation. These observations point to a structured jet as the source of GRB170817A, with a narrow opening angle (θ op ≃ 3.4 degrees) and an energetic core (∼ 3 × 10 52 erg) seen under a viewing angle of ∼15 degrees (Ghirlanda et al., 2019). This is further confirmed by later radio observations Pian Binary neutron star mergers extending up to 300 days after merger, that show a sharp downturn of the radio light curve, suggestive of a jet rather than a spherical source (Mooley et al., 2018c).
The optical/near-infrared kilonova component subsided rapidly (see Section 3.1.2) leaving room to the afterglow emission: the HST observations at ∼100 days after the explosion show a much brighter source than inferred from the extrapolation of the early kilonova curve to that epoch (Lyman et al., 2018). This late-epoch flux is thus not consistent with kilonova emission and rather due to the afterglow produced within an off-axis structured jet (Fong et al., 2019). At X-ray energies, the GRB counterpart is still detected with Chandra three years after explosion (Troja et al., 2020), but its decay is not fully compatible with a structured jet, indicating that the physical conditions have changed or that an extra component is possibly emerging (e.g. a non-thermal aftermath of the kilonova ejecta, see next Section).

The kilonova
The early ground-based optical and near-infrared and space-based (with Swift/UVOT) near-ultraviolet follow-up observations started immediately after identification of the optical counterpart of GW170817, detected a rapid rise (∼1 day timescale, Arcavi et al. 2017) and wavelength-dependent time decay, quicker at shorter wavelengths (Andreoni et al., 2017;Cowperthwaite et al., 2017;Díaz et al., 2017;Drout et al., 2017;Evans et al., 2017;McCully et al., 2017;Nicholl et al., 2017;Tanvir et al., 2017;Utsumi et al., 2017;Villar et al., 2017). The optical light is polarized at the very low level of (0.50±0.07)% at 1.46 days, consistent with intrinsically unpolarized emission scattered by Galactic dust, indicating that no significant alignment effect in the emission or geometric preferential direction is present in the source at this epoch, consistent with expectation for kilonova emission (Covino et al., 2017).
Starting the same night when the optical counterpart was detected, low resolution spectroscopy was carried out at the Magellan telescope (Shappee et al., 2017). This spectrum shows that the source is not yet transparent as it is emitting black body radiation, whose maximum lies however blue-ward of the sampled wavelength range, suggesting that the initial temperature may have been larger than ∼10000 K. The following night (1.5 days after merger) the spectrum is still described by an almost perfect black body law whose maximum at ∼5000 K was fully resolved by spectroscopy at the Very Large Telescope (VLT) with the X-Shooter spectrograph over the wavelength range 3500-24000Å (Pian et al., 2017). At this epoch, the expansion velocity of the expelled ejecta, whose total mass was estimated to be 0.02-0.05 M ⊙ (Pian et al., 2017;Smartt et al., 2017;Waxman et al., 2018), was ∼20% of the light speed, which is only mildly relativistic and therefore much less extreme than the ultra-relativistic kinematic regime of the GRB and of its early afterglow, analogous to the observed difference between the afterglows and the supernovae accompanying long GRBs. At 2.5 days after merger the spectrum starts deviating from a black body as the ejecta become increasingly transparent and absorption lines are being imprinted on the spectral continuum by the atomic species present in the ejecta (Chornock et al., 2017;Pian et al., 2017;Smartt et al., 2017). In the following days these features become prominent and they evolve as the ejecta decelerate and the photosphere recedes (Figure 1). Early attempts at spectral modelling used best guesses and empirical approximate formulae, that are affected by a high level of uncertainty, and a very limited number of individual r-process elements may have been identified (e.g. strontium, Watson et al., 2019).
At ∼10 days after merger, the kilonova spectrum fades out of the reach of the largest telescopes. The radioactive source could still be monitored photometrically for another week in optical and near-infrared (Cowperthwaite et al., 2017;Drout et al., 2017;Kasliwal et al., 2017;Pian et al., 2017;Smartt et al., 2017;Tanvir et al., 2017); it was last detected at 4.5 µm with the Spitzer satellite 74 days post merger (Villar et al., 2018). The kilonova ejecta are also expected to interact with the circum-binary medium

Binary neutron star mergers
and produce low-level radio and X-ray emission that peaks years after the merger (Kathirgamaraju et al., 2019). Search for this component has not returned (yet) a detection at radio wavelengths (Hajela et al., 2019), but it may start to be revealed at X-rays (Troja et al., 2020).

The host galaxy of GW170817
HST and Chandra images, combined with VLT MUSE integral field spectroscopy of the optical counterpart of GW170817, show that its host galaxy, NGC 4993, is a lenticular (S0) galaxy at z = 0.009783 that has undergone a recent (∼1 Gyr) galactic merger. This merger may be responsible for igniting weak nuclear activity (Levan et al., 2017). No globular or young stellar cluster is detected at the location of GW170817, with a limit of a few thousand solar masses for any young system. The population in the vicinity is predominantly old and the extinction from local interstellar medium low. Based on these data, the distance of NGC4993 was determined to be (41.0 ± 3.1) Mpc (Hjorth et al., 2017). The HST imaging made it also possible to establish the distance of NGC4993 through the surface brightness fluctuation method with an uncertainty of ∼6% (40.7 ± 1.4 ± 1.9 Mpc, random and systematic errors, respectively), making it the most precise distance measurement for this galaxy (Cantiello et al., 2018). Combining this with the recession velocity measured from optical spectroscopy of the galaxy, corrected for peculiar motions, returns a Hubble constant H 0 = 71.9 ± 7.1 km s −1 Mpc −1 .
Based only on the gravitational data and the standard siren argument, and assuming that the optical counterpart represents the true sky location of the gravitational-wave source instead of marginalizing over a range of potential sky locations, Abbott et al. (2017d) determined a "gravitational" distance of 43.8 +2.9 −6.9 Mpc that is refined with respect to the one previously reported in Abbott et al. (2017a). Together with the corrected recession velocity of NGC4993 this yields a Hubble constant H 0 = 70 +12 −8 km s −1 Mpc −1 , comparable to, but less precise than that obtained from the superluminal motion of the radio counterpart core, H 0 = 70.3 +5.3 −5.0 km s −1 Mpc −1 .

KILONOVA LIGHT CURVE AND SPECTRUM
The unstable isotopes formed during coalescence of a binary neutron star system decay radioactively and the emitted gamma-ray photons are down-scattered to the ultraviolet, optical and infrared thermal radiation that constitutes the kilonova source (Section 3.1.2). Its time decline is determined by the convolution of radioactive decay chain curves of all present unstable nuclei. This is analogous to the supernova phenomenon, where however the vastly dominant radioactive chain is 56 Ni decaying into 56 Co, and then into 56 Fe.
While radioactive nuclei decay, atoms recombine, as the source is cooling, and absorption features are imprinted in the kilonova spectra. Among neutron-rich nuclei, the lanthanides (atomic numbers 57 to 71) series have full f-shells and therefore numerous atomic transitions that suppress the spectrum at shorter wavelengths ( 8000Å). Spectra of dynamical ejecta of kilonova may therefore be heavily intrinsically reddened, depending on the relative abundance of lanthanides (Barnes & Kasen, 2013;Kasen et al., 2013;Tanaka & Hotokezaka, 2013). Prior to the clear detection of kilonova accompanying GW170817 (Section 3), such a source may have been detected in HST images in near-infrared H band of the afterglow of GRB130603B (Berger et al., 2013;Tanvir et al., 2013). Successive claims for association with short GRBs and kilonova radiation were similarly uncertain (Jin et al., 2015(Jin et al., , 2016. If the neutron stars coalescence does not produce instantaneously a black hole, and a hypermassive neutron star is formed as a transitory remnant, a neutrino wind is emitted, that may inhibit the formation of neutrons and reduce the amount of neutron-rich elements (Fernández & Metzger, 2013;Kiuchi et al., 2014;Metzger & Fernández, 2014;Perego et al., 2014;Kasen et al., 2015;Lippuner et al., 2017). This "post-merger" kilonova component, of preferentially polar direction, is thus relatively poor in lanthanides and gives rise to a less reddened spectrum (Tanaka et al., 2017;Kasen et al., 2017).
The optical/near-infrared spectral behavior of kilonova is analogous to that of supernovae with the largest kinetic energies (> 10 52 erg), like those associated with GRBs: the large photospheric velocities broaden the absorption lines and blueshift them in the direction of the observer. Furthermore, broadening causes the lines to blend, which makes it difficult to isolate and identify individual atomic species (Iwamoto et al., 1998;Mazzali et al., 2000;Nakamura et al., 2001). While these effects can be controlled and de-convolved with the aid of a radiation transport model as it has been done for supernovae of all types (Mazzali et al., 2016;Ergon et al., 2018;Hillier & Dessart, 2019;Shingles et al., 2020;Ashall & Mazzali, 2020), a more fundamental hurdle in modelling kilonova spectra consists in the much larger number of electronic transitions occurring in r-process element atoms than in the lighter ones that populate supernova ejecta, and in our extremely limited knowledge of individual atomic opacities of these neutron-rich elements, owing to the lack of suitable atomic data. First systematic atomic structure calculations for lanthanides and for all r-process elements were presented by Fontes et al. (2020) and Tanaka et al. (2020), respectively.

SUMMARY AND FUTURE PROSPECTS
The gravitational and electromagnetic event of 17 August 2017 provided the long-awaited confirmation that binary neutron star mergers are responsible for well identifiable gravitational signals at kHz frequencies, for short GRBs, and for thermal sources, a.k.a. kilonovae or macronovae, produced by the radioactive decay of unstable heavy elements synthesized via r-process during the coalescence. The intensive and longterm electromagnetic monitoring from ground and space allowed clear detection of the counterpart at all wavelengths. Brief (∼2 s) gamma-ray emission, peaking at ∼200 keV and lagging the gravitational signal by 1.7 seconds, is consistent with a weak short GRB. At ultraviolet-to-near-infrared wavelengths, the kilonova component -never before detected to this level of accuracy and robustness -dominates during the first 10 days and decays rapidly under detection threshold thereafter, while an afterglow component emerges around day ∼100. Up to the most recent epochs of observation (day ∼1000 at X-rays), the kilonova does not add significantly to the bright radio and X-ray afterglow component. Multi-epoch VLBI observations measured -for the first time in a GRB -superluminal motion of the radio source, thus providing evidence of late-epoch emergence of a collimated off-axis relativistic jet.
Doubtlessly, this series of breakthroughs were made possible by the closeness of the source (40 Mpc), almost unprecedented for GRBs, and by the availability of first-class ground-based and space-borne instruments. The many findings and exceptional new physical insight afforded by GW170817/GRB170817A make it a rosetta stone for future similar events. When a sizeable group of sources with good gravitational and electromagnetic detections will be available, the properties of binary systems containing at least one neutron star, of their mergers and their aftermaths can be mapped. It will then become possible to clarify how the dynamically ejected mass depends on the binary system parameters, mass asymmetry and neutron stars equation of state (Ruffert & Janka, 2001;Hotokezaka et al., 2013), how the jet forms and evolves, which kinematic regimes and geometry it takes up in time, and how can the GRB and afterglow observed phenomenologies help distinguish the intrinsic properties from viewing angle effects (Janka et al., 2006;Lamb & Kobayashi, 2018;Ioka & Nakamura, 2019), what is the detailed chemical content of kilonova ejecta and how the r-process abundance pattern inferred from Pian Binary neutron star mergers kilonova spectra compares with the history of heavy elements cosmic enrichment (Rosswog et al., 2018), how can kilonovae help constrain the binary neutron star rates and how does the parent population of short GRBs evolve (Guetta & Stella, 2009;Yang et al., 2017;Belczynski et al., 2018;Artale et al., 2019;Matsumoto & Piran, 2020), how gravitational and electromagnetic data can be used jointly to determine the cosmological parameters (Schutz, 1986;Del Pozzo, 2012;Abbott et al., 2017d), to mention only some fundamental open problems. Comparison of the optical and near-infrared light curves of GW170817 kilonova with those of short GRBs with known redshift suggests infact significant diversity in the kilonova component luminosities (Gompertz et al., 2018;Rossi et al., 2020).
Regrettably, short GRBs viewed at random angles, and not pole on, are relativistically beamed away from the observer direction and kilonovae are intrinsically weak. These circumstances make electromagnetic detections very difficult if the sources lie at more than ∼100 Mpc, as proven during the third and latest observing run (Apr 2019 -Mar 2020) of the gravitational interferometers network. In this observing period, two merger events possibly involving neutron stars were reported by the LIGO-Virgo consortium: GW190425, caused by the coalescence of two compact objects of masses each in the range 1.12-2.52 M ⊙ , at ∼160 Mpc (Abbott et al., 2020a), and GW190814, caused by a 23 M ⊙ black hole merging with a compact object of 2.6 M ⊙ at 240 Mpc (Abbott et al., 2020b). In neither case did the search for an electromagnetic counterpart return a positive result (Coughlin et al., 2019;Gomez et al., 2019;Ackley et al., 2020;Andreoni et al., 2020;Antier et al., 2020;Kasliwal et al., 2020), owing presumably to the large distance and sky error areas, but also possibly to the fact that all coalescing stars may have been black holes, as the neutron star nature of the binary members lighter than 3 M ⊙ could not be confirmed.
The search for electromagnetic counterparts of gravitational radiation signals is currently thwarted primarily by the large uncertainty of their localization in the sky, which is usually no more accurate than several dozens of square degrees. Much smaller error boxes are expected to be available when the KAGRA (which had already joined LIGO-Virgo in the last months of the 2019-2020 observing run) and the INDIGO interferometers will operate at full regime as part of the network during the next observing run (Abbott et al., 2018). Observing modes, strategies, and simulations are being implemented to optimize the electromagnetic multi-wavelength search and follow-up (Bartos et al., 2016;Patricelli et al., 2018;Cowperthwaite et al., 2019;Graham et al., 2019;Artale et al., 2020), and new dedicated space-based facilities are designed with critical capabilities of large sky area coverage and rapid turnaround (e.g. ULTRASAT, Sagiv et al., 2014;THESEUS, Amati et al., 2018, Stratta et al., 2018DORADO, Cenko et al. 2019), to maximize the chance of detection of dim, fast-declining transients.
Finally, the possible detection of MeV and >GeV neutrinos associated with the kilonova (Kyutoku & Kashiyama, 2018) and with the GRB (Bartos et al., 2019;Aartsen et al., 2020), respectively, will bring an extra carrier of information into play, and thus complete the multi-messenger picture associated with the binary neutron star merger phenomenon.

CONFLICT OF INTEREST STATEMENT
The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

AUTHOR CONTRIBUTIONS
The author is fully accountable for the content of this work.  Smartt et al. (2017), at phases indicated in days after merger time, corrected for Galactic extinction E(B-V) = 0.1 mag, de-redshifted, and offset in flux by multiples of a 5×10 −17 erg s −1 cm −2Å−1 additive constant with respect to the 10.5d spectrum. Wavelength ranges of poor atmospheric transmission were blanked out. Some spectra have been re-calibrated with respect to the originally published version, courtesy of J. Gillanders, J. Selsing and S. Smartt.