Gravitational waves from binary neutron stars

I review the current global status of research on gravitational waves emitted from mergers of binary neutron star systems, focusing on general-relativistic simulations and their use to interpret data from the gravitational-wave detectors, especially in relation to the equation of state of compact stars.


Introduction
The composition of neutron stars (a better term would be compact stars) is still largely unknown. Binaryneutron-star (BNS) mergers are prime observational targets for gravitational-wave detectors and one of the main reasons they are being studied is to obtain information on the equation of state of ultrahigh-density matter [13,46], in addition to information on the origin of the abundance of heavy elements in the Universe [104,187,191], on formation mechanisms for magnetars [71,89,105,141] and short gamma-ray bursts [29], and with the goal to test the theory of gravity [185].
The separation of the two neutron stars in a binary decreases in time because of the emission of gravitational radiation. During the merger, when the outer parts of the stars come into contact, a shear surface forms and a Kelvin-Helmholtz instability develops, curling the interface forming a series of vortices [45]. Later, also the cores of the two neutron stars coalesce. The merged object then continues to oscillate, until reaching a configuration that is nearly stationary: a rotating neutron star or, after collapsing, a Kerr black hole. The current estimate of minimum and maximum masses for neutron stars leads to think that most BNS merger would eventually end in a collapse to a black hole. Depending on the initial masses of the stars and on the equation of state (which is still largely unknown for rest-mass densities above twice the nuclear saturation density), this can happen promptly after the merger (i.e., without any bounce of the stellar cores), or after an interval of time. The merger of BNSs with higher masses and/or softer equations of state ends in a prompt collapse to a black hole, while binaries with lower rest masses and/or stiffer equations of state produce a merged object that does not collapse immediately.
In the latter case, the remnant would initially be an unstable, differentially rotating stellar object undergoing violent oscillations. Differential rotation and thermal gradients [17] prevent prompt collapse also for remnants with masses above the upper limit for static neutron stars (or, more precisely, for uniformly rotating neutron stars, whose mass limit is slightly higher than that of nonrotating neutron stars). Such configurations resulting from BNS mergers are often called hypermassive neutron stars (HMNS), even if the original term refers to equilibrium models [17], while the merger remnant is not in equilibrium at this stage. Supramassive neutron star (SMNS) is used to refer to a uniformly rotating neutron star with rest mass larger than the upper limit for static L. Baiotti (B) International College and Graduate School of Science, Osaka University, 1-2 Machikaneyama-cho, Toyonaka, Japan E-mail: baiotti@ipc.phys.sci.osaka-u.ac.jp neutron stars in equilibrium. Strictly speaking, also the term SMNS refers only to equilibrium configurations [52,53].
The dominant overall deformation of the material object formed after the merger has m = 2, i.e., it is a bar-like deformation. 1 Gravitational waves carry away further energy and angular momentum from the bar. Non-axisymmetric deformations (m = 1) are also seen in simulations [9,30,54,63,72,73,106,123,145,148,149,156]. A secular increase of the central rest-mass density is also observed.
Remnants supported against collapse by differential rotation and thermal gradients may last for up to ∼ 1 s. The remnant is driven toward collapse or toward a uniformly rotating neutron star by several processes: angular-momentum transport related to magnetic-field effects (magneto-rotational instability and magnetic braking), cooling through neutrino emission, and the gravitational torque resulting from its non-axisymmetric structure (see, e.g., Refs. [99,168,173]). If a uniformly rotating star forms, it is likely that its mass is above the limit for nonrotating stars, and thus, it would then collapse to a black hole on timescales (∼ 10 − 10 4 s) related to angular-momentum loss due to secular mechanisms (mostly electromagnetic emission) [160]. If, instead, the total mass of the remnant is smaller than the upper limit for a static neutron star, the merged object may never collapse; namely, it may become a stable neutron star. This is probably a marginal case. See for example Fig. 1 of Ref. [14] for a schematic view of what said in this paragraph about the possible outcomes of the merger.
When a black hole is formed, a large amount of high-angular-momentum matter may remain outside of it in the form of an accretion disc. The projection of the disc on the equatorial plane is predominantly circular, but, because of its violent birth, it undergoes large oscillations, mostly in the radial direction. On average, such accretion discs have a peak rest-mass density around 10 12 − 10 13 g cm −3 , a horizontal extension of a few tens of km, and an aspect ratio of 1/3. The initial rest mass of the disc varies in the range 0.001 − 0.2M and influences the ejecta from the BNS and electromagnetic emission, including that of kilonovae [134] and short gamma-ray bursts, for which the existence of a massive disc around the newly formed Kerr black hole is thought to be necessary [29].
To help and interpret observations, we need solutions of the general-relativistic equations describing at least the dynamics of spacetime, matter, and magnetic fields, and possibly additional ingredients like neutrino radiation. Such equations need to be solved numerically. After overcoming big theoretical and technical difficulties in the discretization of the Einstein equations, 2 numerical relativity can nowadays provide some reliable results [14,70]. Codes can robustly compute the matter and spacetime dynamics for tens of orbits before the merger and hundreds of milliseconds after the merger, also thanks to the selection of appropriate, evolved gauges. The equations of state used in simulations are either piecewise polytropes [161] (often with the addition of a thermal part) or tabulated from nuclear-physics studies. The first general-relativistic simulations of merging neutron stars including quarks at finite temperatures have been performed recently [22,[24][25][26]35,90,139,140,199].
Some robust results have been obtained also in computations of ejecta from BNS mergers (especially for estimates of heavy-element production and kilonovae light curves) [176] and in the general-relativistic treatment of physical viscosity for BNSs [8,86,87,109,154,177,209]. Some other issues, instead, are still very open, in particular the effects of magnetic fields [14] and of radiation transport [193].
In the following, I will focus on how it is or may be possible to extract information on the neutron-star equation of state from analysis of gravitational radiation emitted from BNS mergers. The two main methods to link the observed gravitational waves to the neutron-star equation of state involve either tidal deformations during the last orbits before merger or the spectra of the gravitational radiation emitted from the post-merger object (if collapse to a black hole does not occur too early). Tidal deformations have been measured from GW170817 [1,183,184] and GW190425 [186], but post-merger gravitational waves have not been observed yet. 1 Stellar deformations are commonly described by decomposing the linear perturbations of the energy or rest-mass density as a sum of quasinormal modes that are characterized by the indices (l,m) of the spherical harmonic functions. The m mentioned in the text is the dominant term of such expansion. m = 0 is a spherical perturbation, m = 1 is a one-lobed (or one-armed) perturbation, and m = 2 is a bar-shaped perturbation. See Ref. [148] for a review. 2 In addition to the standard problems of any numerical simulation, there are multiple reasons for the increased difficulty inherent to general-relativistic simulations: (i) the formulation of the equations is not self-evident (time is not "simply" defined and very careful definitions of variables are needed to obtain a system that is strongly hyperbolic); (ii) physical singularities may be present and need special treatment; (iii) while not carrying physical information, gauges play an important role in numerical stability, for example in countering grid stretching. For explanations on all these issues, well-written textbooks are available [5,16,91,167,175].

Gravitational waves from the inspiral phase
Stars in a binary system are tidally deformed, and these deformations affect the orbital trajectory of the binary and thus the emitted gravitational waves [83,111,119]. Such tidal deformations depend also on the equation of state and in particular on the pressure gradients inside the star, which are related to the stiffness of the equation of state. It has also been shown that differences between equations of state of neutron stars that do not appear as differences in the tidal deformability are difficult to measure through inspiral gravitational-wave signals [162].
The lowest-order description of tidal deformations can be achieved by introducing a single tidal deformability coefficient, defined as the proportionality constant λ between the quadrupole moment of the star Q i j and the external tidal field E i j (the field generated by the companion star) (1) or, better, the dimensionless tidal deformability 3 where M is the stellar mass. Equivalently, this can be written as where κ 2 is the quadrupole Love number (whose values for different hadronic equations of state and masses are in the range 0.02 − 0.15 [151]) and R is the stellar radius. can be calculated as [58,96,201] where C ≡ M/R is the stellar compactness and y(r ) satisfies the differential equation where p and ρ represent pressure and mass density, respectively, and where γ (r ) is in the definition of the metric coefficients Equation (5) can be solved for a given equation of state together with the Tolman-Oppenheimer-Volkoff (TOV) equations [144], which describe spherically symmetric stars in static equilibrium in general relativity. Appropriate boundary conditions need to be imposed. The dimensionless tidal deformability for neutron stars may take values between a few and ten thousands, depending on the stellar mass and equation of state. A larger tidal deformability corresponds to a larger, less compact star that is easily deformable and vice versa. For a given neutron-star mass, the values of the dimensionless tidal deformability vary up to a factor of ten according to the equation of state assumed.
In practice, the dimensionless tidal deformability cannot be measured with the sensitivity of current and near-future detectors. 4 What has been and in the near future will be directly measured from BNS waveforms is the dominant tidal parameter in the post-Newtonian expansion of the waveform, where the gravitational-wave signal can be calculated approximately by imposing that the power radiated by a binary system in gravitational waves is equal to the change in the energy of the binary. Then, the lowest-order tidal effects appear as a term of the fifth post-Newtonian order, whose coefficient is given the name of mass-weighted average tidal deformability (also called effective tidal deformability) and has the expression [82,83,195] = 16 13 or, equivalently˜ where q ≡ m B /m A (≤ 1) is the ratio of the masses of the two stars in the binary.
A different method for obtaining information on the tidal deformability from detector data consists in extracting the (mass-independent) coefficients of a Taylor expansion of the tidal deformability about some fiducial mass [62,133,197,[202][203][204]; usually, 1.4M or 1.35M . Current detectors, however, do not have sufficient sensitivity for measuring any of these coefficients accurately enough.
To estimate the parameters of inspiralling BNS systems, it is necessary to perform matched filtering, which require theoretically predicted template waveforms (also called approximants), one for each value of the physical parameters considered in the analysis. Currently, the standard approximants in gravitational-wave analysis are post-Newtonian expansions, effective-one-body (EOB) [42,43], and the Phenom models [3,4]. The lowestorder post-Newtonian tidal effects appear as terms of the fifth order, but currently relevant post-Newtonian terms are known only up to fourth order in the dynamics [59,129] and only to order 3.5 in the gravitationalwave phasing [36]. This may lead to systematic errors [82,195,205], which other phenomenological methods have tried to bypass.
The EOB model [42,43] maps the two-body dynamics to those of a single effective particle in an effective potential. The dynamics and gravitational emission from a binary are computed by solving the coupled system of ordinary differential equations for the orbital motion, gravitational-wave generation, and radiation backreaction in the time-domain. By comparing with numerical-relativity simulations, refinements and calibrations of the models can be carried out [12,37,55,142]. It is thought that further improvements of the EOB models are currently necessary for a satisfactory description of the signal [64,67,100,169], at least for some types of BNS systems.
As the name states, Phenom models [3,4,108,172] are other phenomenological approaches, which approximate a set of hybrid waveforms constructed by matching numerical-relativity waveforms with analytical post-Newtonian waveforms [3,4,65,108,172]. Finally, reduced-order-modeling techniques [117] and others [15,115,116,147] have been proposed for faster computations.
Estimates of the tidal deformability have indeed been possible by analyzing GW170817 and GW190425. In these two first detections, Bayesian analyses indicated that the systematic uncertainties due to the modeling of matter effects (choice of representation of the equation of state and choice of approximant) are smaller than the statistical errors in the measurement [143,184,186]. However, further improved waveform approximants will be necessary to avoid related systematic errors in future higher signal-to-noise-ratio events [57,66,69,92,95,107,169].
More comprehensive Bayesian analyses, including not only gravitational-wave data, but also data from the related electromagnetic emission and other observational and experimental data are also available [23,41, 56,93,122,130,155,158,196].

Gravitational waves emitted after the merger
After a merger, the maximum density in the system increases of a factor of several and temperature increases up to ∼ 50MeV. The post-merger phase is particularly interesting also because such high-density hightemperature conditions cannot be easily investigated anywhere else. The energy emitted in gravitational waves after the merger may be higher than that emitted during the inspiral (if there is not a prompt collapse), but post-merger gravitational radiation is more difficult to detect, because its frequencies are higher (from 1 to several kHz [21,127,148,179]), and thus their signal-to-noise ratio in current detectors is smaller than that of the inspiral. They are unlikely to be measured by detectors like Advanced LIGO or Advanced Virgo, namely only for extremely rare sources within ∼ 30 Mpc [47,50,51,75,188,190,206].
Furthermore, turbulence, large magnetic fields, various physical instabilities, neutrino emission, phase transitions, viscosity, and other microphysical effects make numerical simulations of the merger and postmerger dynamics more difficult than those of the inspiral part. As a consequence, the accuracy of simulations of the post-merger is not so high as that achieved for the inspiral. In particular, while the spectrum of post-merger gravitational radiation is known somewhat reliably, currently, it is not so for the phase of post-merger signals (see, e.g., Ref. [40]). For these reasons, it is challenging to construct accurate templates of post-merger signals.
The basic idea for obtaining information on the supranuclear equation of state from gravitational waves emitted after the merger is that the main peak frequencies of the post-merger power spectrum correlate, in a rather equation-of-state-insensitive way, with properties (radius at a fiducial mass, compactness, etc.) of a zero-temperature spherical equilibrium star. These relations may be used to estimate constraints that future observations of the post-merger signal may set on the neutron-star radius [39,40,47,188,190].
The general structure of gravitational-wave spectra of BNS mergers consists of a lower-frequency part related to the inspiral, and a part at higher frequencies, with some peaks. Through numerical simulations of BNS mergers, the frequency of the main peak, f main peak , which corresponds to the fundamental mode of oscillation of the star, has been found to correlate with the tidal deformability [166,180], or with the radius of a 1.6M neutron star [21] with the given equation of state. These relations are affected by the spin of the neutron stars in the inspiral (see, e.g., Refs. [20,30,31,64,66,74]), but otherwise hold at the ≈ 10% level for equations of state without phase transitions and for both equal-mass and unequal-mass binaries [110]. The likely physical explanation of these empirical relations is that at a given separation, the tidal interaction is more attractive for larger values of the tidal deformability and thus such binaries merge at lower frequencies (larger separations). As a consequence, their remnants are less bound, and have larger angular-momentum support and extension at formation and thus lower frequency of the fundamental mode [31].
The location of the main peak is estimated to have an error of about 5% in numerical simulations at the available resolutions [110]. One of the secondary peaks at lower frequency appears only in (nearly) equal-mass BNS mergers and for softer equations of state [110], but if present and observed, its frequency might be used to gain information on the equation of state, according to a relation with the average compactness of the original stars in the binary [179,180] or with the tidal deformability [166].
Another smaller peak in the spectrum may be observed in relation to an m = 1 deformation (also called one-armed spiral instability; see note on page page 2) appearing in the merged object [148], particularly in the case of eccentric binaries [73]. This deformation, showing up in the spectrum at about half the frequency of f main peak , is found to be present generically in BNS simulations and to carry information about the equation of state, but its amplitude is much smaller than that of other peaks and it is unlikely that current gravitational-wave detectors measure it [156]. Future detectors may have the sensitivity to observe these signals.
The location of most other peaks 5 depends more strongly on numerical parameters and thus cannot be easily used for estimating physical properties.
Note that the merger remnant rapidly evolves toward a more stable configuration and, thus, especially in the first milliseconds after the merger, post-merger frequencies slightly change in time (see spectrograms in, e.g., Refs. [51, 61,64,85,127,166,180]). Prominently, f main peak increases up to the collapse to black hole, while its amplitude decreases. This is because the merged object loses angular momentum to gravitational radiation, and otherwise, the angular-momentum distribution is rearranged, so that axisymmetry and compactness increase [127,166,180]. To track the time evolution of these frequencies and thus assert reliably the spectral properties of the gravitational-wave signal in observations, a high enough signal-to-noise ratio is required.
A different way to get information on the mass-radius relation of neutron stars is based on some other empirical relations, obtained from numerical simulations of BNS mergers, between the threshold binary mass for prompt collapse, the maximum mass for a nonrotating neutron star, and its radius [18,19,23,25,26,113]. This method is said to be more robust, because whether a merger ended in prompt collapse or not may be assessed not only through gravitational radiation, but also through associated electromagnetic emissions, which may be observed more easily. The only other quantity this method relies on is the chirp mass, which can be measured precisely from inspiral gravitational waves.
The full LIGO-Virgo network operating at design sensitivity may provide in the future reasonable estimates of the dominant post-merger oscillation frequency [184]. In about 1 year of operation of third-generation detectors [153,165], the dominant-peak frequency may be measured to a statistical error of the order of 10Hz for certain equations of state, corresponding, through the above-mentioned empirical relations derived from numerical simulations, to a radius measurement to within tens of meters [206]. With further improvements, it will probably also be possible to extract subdominant frequencies.
It is to be noted that in this kind of analyses the error on the radius is currently dominated by the scatter in the relations coming from numerical relativity, rather than the statistical error of the signal reconstruction itself for signal-to-noise ratios larger than ∼ 10 (see, e.g., Refs. [40,47,75,188,190,206]).

Constraints on the high-density equation of state from GW170817 and other observations and experiments
Tens of works have analyzed the data from GW170817 (some of them together with other observational and experimental data) to set constraints in the space of possible high-density equations of state. Data from GW190425 have been found not to further affect constraints much (see, e.g., Ref. [68]). One of the ways to set such constraints consists in estimating from data the relation between stellar mass and radius. The stellar masses in a binary are known more directly from gravitational-wave data, but computation of radii needs further assumptions, including some on the type of equations of state considered in an analysis.
Common findings are that the upper bounds on radius and tidal deformability (ruling out the stiffer equations of state) are driven by gravitational-wave observations, while the lower bounds on radius and tidal deformability (ruling out the softer equations of state) are driven by the observation of high-mass pulsars and neutron-star radii by NICER. In most cases, estimates from GW170817, nuclear theory and experiments, X-ray observations, and combinations of them are found to be compatible. One exception are studies that analyze together data from GW170817, NICER, and the PREX-II experiment [2].
PREX-II is the second in the series of the Lead Radius Experiment, based on elastic scattering of polarized electrons. It has measured the neutron-skin thickness R skin of 208 Pb, namely, the difference between the distribution of the neutron root-mean-square radius and the proton root-mean-square radius in 208 Pb. It is known in nuclear physics 6 that the slope of the symmetry energy L displays a strong correlation with the neutron 6 To connect more directly with nuclear theory and experiments, nuclear physicists often use an empirical parameterization of the nuclear equation of state based on the Taylor expansion of the energy per nucleon E in terms of the isospin asymmetry parameter δ ≡ (ρ n − ρ p )/ρ (where ρ p and ρ n are the proton and neutron number densities and ρ ≡ ρ p + ρ n ) around the symmetric nuclear matter case δ = 0 and near the saturation density ρ 0 (see, e.g., Refs. [44,80,128,131,132,210,212]), which is ρ 0 2.7 × 10 14 g cm −3 (or, equivalently, 0.16 nucleons per fm −3 ). The expression is [192] where E(ρ, 0) corresponds to the energy of symmetric nuclear matter. E(ρ, 0) and E sym (ρ) are then expanded around the nuclear saturation density where y ≡ (ρ −ρ 0 )/3ρ 0 . These lowest-order parameters are known as the energy per particle E 0 , the incompressibility coefficient K 0 , the symmetry energy E sym,0 , its slope L, its curvature K sym , and its skewness Q sym at saturation density. Some of these parameters (E 0 , K 0 , and E sym,0 ) are either strongly constrained by terrestrial experiments or have no noticeable impact on the bulk properties of neutron stars (see, e.g., Ref. [124] for a review.).
skin thickness and there are hints that the symmetry energy L is also correlated to the radius of a neutron star with a given mass. Therefore, there have been attempts to carry out Bayesian studies comparing the tidal deformability observation from GW170817, neutron-star radii from NICER observations, and neutron-star radii evaluated through PREX [34, 79,163]. A mild tension has been found between astrophysical observations and PREX measurements. PREX finds a neutron-skin thickness around R skin = 0.283 ± 0.071 fm [2], while astrophysical observations suggest R skin 0.2 fm. Namely, experiments with atomic nuclei suggest a stiffer equation of state than the one implied by astrophysical observations for the higher densities encountered in neutron stars. The evolution of the equation of state from being stiffer at lower densities to being softer at higher densities may be indicative of a phase transition in the interior of neutron stars [34, 79,163].

Investigating phase transitions
Up to densities as high as twice the nuclear saturation density, the neutron-matter equation of state obtained by theoretical arguments and corroborated by laboratory experiments is expected to be reliable (see Ref. [28] for a review). Beyond that not much is known for sure and, as mentioned previously, matter could undergo a phase transition [27, 101,102] or a smooth crossover [28] to hyperonic or quark matter, possibly giving rise to a hybrid hadron-quark star, in which the equation of state is purely nucleonic for lower densities and contains deconfined quark matter at higher densities. The possibility of more than one phase transition has also been taken into account [6,7,28,94].
Such phase transitions would change drastically the properties of the star, in particular because they may produce drastic softening (or perhaps stiffening) of the equation of state. Therefore, neutron stars with similar masses but rather different tidal deformabilities may be observed and yield information about transitions to quark matter.
Tens of works have analyzed the data of GW170817 without assuming the absence of phase transitions as instead done in the original analyses and found that the data are also consistent with the coalescence of a BNS system containing at least one hybrid star (for a review, see, e.g., Refs. [13,46]).
Since post-merger density is higher than that of the parent stars, it is possible that phase transitions occur only after the merger. In this case, a measurement of the tidal deformabilities, of course, cannot contain information on phase transitions.
The easiest way to identify a phase transition that has a sufficiently strong discontinuity in the density, like one from nucleonic to deconfined quark matter, consists in inspecting the main frequency of the power spectral density of the post-merger phase, f main peak , which may be different from the one inferred from the measurement of the tidal deformability in the inspiral (namely, before the phase transition occurs) [35,139] or may change rapidly during the post-merger phase because of the abrupt speed-up of the rotation of the differentially rotating core of the remnant [24], related to the formation of a quark-matter core with a softer equation of state that makes the remnant more compact. The detection of such an abrupt change in the dominant post-merger gravitational-wave frequency would also allow to constrain the density at which the phase transition occurs and might be possible in future gravitational-wave detectors [24].
The lifetime of the merged object before collapse and the black-hole ringdown waveforms may also be rather different from the ones expected from purely nucleonic equations of state [139,140] and, in addition to changes in gravitational waves, the rearrangement of the angular momentum in the remnant resulting from the formation of a quark core could be accompanied by a significant burst of neutrinos and then by a gamma-ray burst [38,49,136], which, however, might not be qualitatively different from the electromagnetic counterparts of neutron-star mergers where phase transitions do not occur [22].
Another way to identify the presence of strong phase transitions is based on the relation between the total mass M threshold of the BNS that is at the threshold between prompt and delayed collapse to black hole and the dimensionless tidal deformability threshold of a star that has a mass M = M threshold /2. It has been found that only BNS that have undergone a phase transition can occupy a part of the plane (M threshold , threshold ) [25,26].
An attempt to classify different dynamics of post-merger evolutions in the presence of phase transition has been offered in Ref. [199].
Up to this point, I have focused on equations of state with phase transitions to deconfined quark matter. Another type of transitions that has received a lot of attention is the one to hyperonic matter. However, numerical-relativity simulations indicate that for BNS systems described with hyperonic equations of state, both the inspiral and the post-merger gravitational waveforms are, in the foreseeable future, indistinguishable from those obtained with the corresponding purely nucleonic equation of state [24,157]. Amplitude, phase modulation, and total luminosity (but not frequency) of the post-merger gravitational-wave signal are found to be different for hyperonic equations of state, but it is very difficult to extract information from these quantities, because they are more difficult to measure in detectors and more difficult to simulate accurately in numerical theoretical studies [157].  203 (1992)