Interpreting Binary Neutron Star Mergers: Describing the Binary Neutron Star Dynamics, Modelling Gravitational Waveforms, and Analyzing Detections

Gravitational waves emitted from the coalescence of neutron star binaries open a new window to probe matter and fundamental physics in unexplored, extreme regimes. To extract information about the supranuclear matter inside neutron stars and the properties of the compact binary systems, robust theoretical prescriptions are required. We give an overview about general features of the dynamics and the gravitational wave signal during the binary neutron star coalescence. We briefly describe existing analytical and numerical approaches to investigate the highly dynamical, strong-field region during the merger. We review existing waveform approximants and discuss properties and possible advantages and shortcomings of individual waveform models, and their application for real gravitational-wave data analysis.


Introduction
The collision of two neutron stars (NSs) belongs to the most violent events in the Universe and is connected to a variety of observables from gravitational wave (GW) and neutrino signatures up to electromagnetic (EM) signals ranging through the full frequency band. This variety of possible 'information channels' makes binary neutron star (BNS) mergers a perfect 'tool' to study supranuclear dense material and strong gravity regions.
The discovery of gravitational waves (GWs) [7] and electromagnetic (EM) signals [6,286,8] arising from a BNS coalescence detected in August 2017 marked a breakthrough in the field of multi-messenger astronomy. This multimessenger observation represented a major leap forward in a number of research areas. It enabled a new and independent measurement of the Hubble constant [5,278,11,171,136,177]; it proved NS mergers to be a major cosmic source of heavy r-process elements, e.g., [138,494,286,285]; it revealed that BNS mergers can act as the central engine for gamma-ray bursts (GRB) [8]; it enables an accurate measurement of the propagation speed of GWs [8]; it ruled out a number of alternative theories of gravity [200,39,139]; and it allowed to place limits on the equation of state (EOS) of cold matter at supranuclear densities, e.g., [7,55,30,386,467,361,462,486,10,12,167,447,135,137,141,445,106,177].
In addition, the increasing number of potential BNS candidates and the second confirmed detection of a BNS merger, GW190425 [13], suggest that many more BNS systems will be detected in the near future.
The extraction of information from the GW (and potentially EM) data relies on existing theoretical waveform models to perform a matched-filtering analysis, in which the data are cross-correlated with templates covering the possible physical parameter space, e.g. [530]. Cross-correlations have to be performed for a sufficient number of individual target waveforms (about a hundred millions of individual waveforms), which leads to large computational costs. Therefore, the computation of each individual waveform needs to be efficient and fast to ensure that the Bayesian parameter estimation of signals, containing several thousand GW cycles, is at all manageable. On the other hand, the construction of accurate template models requires a detailed understanding of the late stages of the BNS coalescence, where matter effects on the GW signal become increasingly prominent. In this regime, one has to solve Einstein's equations together with the equations of general relativistic hydrodynamics (GRHD). This challenging task is typically accomplished through numerical relativity (NR) simulations, see e.g. [221,201,38,180,35,295] and references therein.
Unfortunately, NR simulations are computationally expensive and cannot be used to create template banks for BNS GW signals. Therefore, other approaches are needed to search for and interpret measured GW signals. In the last years, there has been an increasing effort and success to construct analytical and semi-analytical models to describe the GW signals emitted during the BNS coalescence. Despite the improvement, including the computation of higher-order tidal corrections, dynamical tidal effects, and spin-tidal couplings, e.g., Refs. [163,414,41,15,321,283], the performance of Post-Newtonian (PN) waveform approximants becomes increasingly inaccurate towards the merger, e.g. [67,208,536,276,183,189,470], which shows the need for possible alternatives. Therefore, most state-of-the-art time-domain tidal waveform models [36,37,66,274,268,497,395,392,22,390] are based on the EOB description of the general relativistic two-body problem [100,155]. This approach has proven to be able to predict the BNS merger dynamics in large regions of the BNS parameter space, but recent NR data revealed configurations for which further improvements of the tidal EOB models might be required [274,178,22]. While one can expect that over the next years, these issues will be overcome due to further progress in the fields of NR, gravitational self-force (GSF), PN theory, and applications of methods from high-energy physics, the high computational cost for a single EOB waveform is yet another disadvantage. One possibility to speed up the EOB computation is the use of high-order post-adiabatic approximations of the EOB description to allow an accurate and efficient evaluation of the waveform up to a few orbits before merger [392]. The other possibility is the construction of a surrogate model [312,314]. Those models allow the fast computation of waveforms in the frequency domain and are well suited for direct use in parameter estimation pipelines.
In addition to PN and EOB approximants, alternative ways to describe tidal GW signals have also been proposed. Refs. [313,418,509] develop phenomenological black hole-neutron star (BHNS) approximants based on NR data. Ref. [48] transforms NR simulations of binary black hole (BBH) systems by adding PN tidal effects, and Refs. [326,325] develop a method to employ NR waveforms or computationally expensive waveform approximants, e.g. tidal EOB approximants, directly for parameter estimation.
Yet another approach to describe BNS systems was presented in Ref. [174], in which existing BBH models have been augmented by an analytical closedform expression correcting the GW phase to include tidal effects [174,183,290,181].
In this review article, we discuss the BNS and merger dynamics in Sec. 2. In Sec. 3, we briefly describe the use of NR simulations focusing on the imprint of the binary properties on the merger outcome and we review the current status in producing high-quality NR waveforms. In Sec. 4, we give a brief overview about the existing PN knowledge for the description of the BNS coalescence. In Sec. 5 and Sec. 6, we review EOB and phenomenological BNS models, respectively. Section 7 gives a short overview about the usage of BNS approximants for the parameter estimation of BNS signals. Throughout

The BNS coalescence
In contrast to BBH signals, BNS and BHNS mergers can produce a variety of EM counterparts such as sGRBs, a synchrotron afterglow, e.g., [194,404,398,331,396,473], as well as a kilonova, see [503,377] for reviews. The main difference between the GW signals emitted from the coalescence of BNS or BHNS systems and BBH systems are finite size effects that change the coalescence before the compact objects merge. Additionally, for BNS systems there is the possibility to have a postmerger GW signature. An illustration of the density evolution and the emitted GWs of a BNS merger are presented in Naturally, the full GW signal is composed of an inspiral and a postmerger evolution. These two regimes are separated by the moment of merger, i.e., the time when the amplitude of the GW signal reaches its maximum; cf. Fig. 1.

Inspiral
The inspiral of a BNS coalescence reaches frequencies up to about 1-2kHz, where the exact frequency depends on the binary properties. Current GW spin-orbit spin-spin signal-to-noise symm. mass ratio tidal chirp mass Fig. 2 Accumulation of information about mass, spin, and tidal parameters per logarithmic frequency interval. Left panel: aLIGO design sensitivity (zero detuned high power configuration) [4,3], right panel: Einstein Telescope [265,1]. Here, 'normalized measure' is the leading-order contribution to the integrands appearing in the Fisher information matrix, with each curve normalized to its maximum value. To generate these plots, we used the PN TaylorF2 model with the nonspinning point-mass terms up to the highest order available (3.5PN) [474,378] combined with the leading-order tidal [214] and spin effects [292,146,88]. Adapted from Ref. [256].
detectors are most sensitive at frequencies around ∼ 100Hz and get less sensitive with increasing frequencies. Because of this frequency-dependent sensitivity only the inspiral phase of BNS coalescences have been detected for GW170817 [9,12] and GW190425 [13]. While the detection of the postmerger signature is anticipated in the next years or decades, most information about the binary properties will come from the measurement of the inspiral. Therefore, the GW community has made significant efforts for a proper modeling of this regime. Here, we will focus mainly on the dominant tidal effects arising from the deformation of the NSs due to the external gravitational field of the companion, and matter effects introduced by the intrinsic rotation of the stars leading to a deformation related to the rotation frequency of the individual stars. These individual components affect the GW signal at different frequencies, cf. Fig. 2. While the pure tidal effects dominate the late inspiral, the spin-induced quadrupole moment is mostly detectable from the GW signal around 20−30Hz or at even lower frequencies for 3rd generation GW detectors.

Coupling of NS matter to GWs
During the early inspiral, the separation of scales in the binary system can be exploited for analytical approximations that enable tracing the information flow from matter properties to GW signals. In these approaches, the challenge of computing the dynamical spacetime is split into a coupled system of the binary dynamics, the asymptotic GWs, and the backreaction of GW losses on the dynamics. This requires using different approximation schemes in different regions of the spacetime and over different timescales. They are connected into a composite description through matched asymptotic expansions. For instance, the binary dynamics involve meshing a description adapted to the body zones close to each NS, where gravity is strong but the effects from Fig. 3 Cartoon depicting the definition of tidal deformability. The tidal field E due to the spacetime curvature of the companion causes the NS to deform as the matter adjusts to a new equilibrium configuration. The relevant quantity influencing the GWs is the induced change in the multipole structure of the NS's exterior spacetime Q. The multipoles are also impacted by spin effects, and dynamical tidal effects.
the presence of a companion is small, with a description of the interaction zone where the NSs behave almost as point masses with small corrections due to their finite size [213,442], see also [184,513,281,303,302,535]. For weakly self-gravitating bodies described by PN gravity see also the seminal series of papers by Damour, Soffel, Xu [164]. As will be discussed in detail in Sec. 4, the multipole moments defined for the spacetime in the vicinity of the NSs play a key role for communicating information about NS matter between these descriptions. The multipole structure is affected by a variety of tidal effects, spins, and more complicated spin-tidal interactions. In addition to affecting the dynamics, the NS' multipole moments also give rise to additional imprints on the asymptotic gravitational radiation. The radiation can be described by double perturbation expansion around flat spacetime and an infinite series of radiative multipole moments, as explained in detail in the review article [86]. The radiative moments are related in a complicated way, i.e., nonlinearly and non-locally in retarded time, to the total multipole moments of the binary system, which comprise contributions from the orbital motion and the NSs' multipoles. Problems such as the relativistic two-body problem that involve different scales can also efficiently be treated with effective-field-theory methods, see [335,436,466,219] for comprehensive reviews and references.

Dominant tidal effects
In Newtonian gravity, tidal effects arise from the response of the NS to the gradient of the companion's gravitational field across its matter distribution. From the perspective of the NS, the companion is orbiting and produces a timevarying tidal field that slowly sweeps up in frequency. This quasi-periodic tidal forcing can excite characteristic oscillation modes in the NS that depend on the properties of matter in its interior. These concepts carry over to a General Relativistic description, where the modes are the NS's quasi-normal modes. A NS has a broad spectrum of modes [300], several of which have sufficiently low frequencies to be relevant for the inspiral. The tidal excitation can either be a  full resonant excitation of a mode or can be approximately adiabatic when the mode is driven well below its resonance frequency. The net effect of a mode on the GWs depends on (i) its frequency, which determines at what stage of the inspiral the mode gets excited and how long the effect accumulates, and (ii) on the tidal excitation factor, which characterizes how strongly the mode couples to the tidal field (analogous to the tidal deformability coefficients for the fundamental modes, as will be explained below). For circular binary inspirals, the fundamental modes have by far the largest tidal coupling [371,484,299,316], with the effect of the quadrupole = 2 being the largest. During most of the inspiral the fundamental modes generally have a higher frequency than the circular orbit tidal driving frequency and the effect is dominated by the adiabatic contributions characterized by the tidal deformability parameters [214]. The tidal deformability λ for each -th multipole is defined in the limit that the timescale of variations in the tidal field are much faster than any internal timescales of the NS. It is the ratio of the induced multipole deformation in the NS's exterior spacetime Q to the strength of the multipolar tidal field E associated with the curvature due to the companion as will be discussed in Sec. 4: The associated effect on GWs steeply increases with smaller separation or higher GW frequency, however, the information accumulates over all frequencies.

Spin-induced multipole effects
Multipole moments that couple to the companion's tidal field and lead to EOS signatures in GWs also arise from the NS's rotation. For rotating black holes, which are vacuum gravity in General Relativity, the multipole structure of the exterior spacetime is fixed by the no-hair property in terms of only their mass and spin [113,261,247]: where M BH and S BH are the -th mass and current multipole moments of the black hole, and m and χ = cS/Gm 2 are its mass and dimensionless spin. For material objects like NSs, the structure (2) no longer applies and instead the quadrupole and higher multipole moments also depend on the EOS of the matter. For instance, the spin-induced quadrupole is given by The key characteristic matter-dependent coefficients κ are the rotational Love numbers [257,311,432,381,72]. Studies in Refs. [544,423] have shown that the complicated multipole structure of a spinning NS can be approximated by three moments. The effects of spin-induced multipole moments in GWs are quadratic and higher order in the NS's spin (see [306] for the latest update on existing knowledge) but scale as a lower power of the frequency or orbital separation that fundamental-mode tidal effects.

Post-merger dynamics
In contrast to BBH systems where the final state is always a BH, the postmerger of BNS allows for a variety of possible merger outcomes. In general the dynamics and the corresponding GW spectrum is complicated and influenced by several physical processes as thermal effects, turbulences, the magnetohydrodynamical instabilities, neutrino emissions, and possible phase transitions, e.g., [492,23,443,487,51,382,168,130,540]. Therefore, limited knowledge about this part of the BNS coalescence exists and NR studies are the only way of investigating this regime of the BNS coalescence. However, we note that even state-of-the-art simulations contain just a subset of the important physical properties influencing the postmerger and lack a clean error budget during this part of the simulation due to shocks, turbulences, or discontinuities present at or after the merger; see e.g. [383] for recent progress improving the accuracy of postmerger simulations. Independent of these caveats, the qualitative picture, as depicted in Fig. 5, is robust and the possible merger outcomes are generally classified in four different categories [50]: prompt collapse, hypermassive NS (HMNS), supramassive NS (SMNS), and massive NS (MNS).
If the total mass of the binary is large enough, then the system will undergo a prompt collapse to a BH, i.e., the collapse happens within a few milliseconds after the merger. In this scenario, only a small amount of matter is ejected and the final BH is surrounded by a low mass accretion disk. Therefore, one does not expect bright EM counterparts, e.g, [55,447,362,135,177], unless the system has a large mass ratio [295]. The distinction between prompt collapse configurations and non-prompt collapse configuration is mainly determined by the ratio of the total mass M and a given threshold mass M thr as computed in [52, 304,18].
If the mass is smaller than the threshold mass and the baryonic mass is below the baryonic mass supported by a spherical, non-rotating star, the remnant is called a MNS. If the baryonic mass is below the maximum baryonic mass of a rigidly rotating star, it is called a SMNS, otherwise, if the remnant is supported by differential rotation, it is classified as a HMNS. It is important to emphasize that the latter three classifications were derived for cold-equilibrium configurations with zero temperature EOSs [50], but has been applied to make qualitatively predictions, e.g. [272,558]. A proper classification incorporating the temperature dependent effects are currently missing, but see [284,105] for first attempts.
While this review focuses on the inspiral part of the BNS coalescence, we will briefly present some characteristic features of the postmerger evolution, which for future generations of GW detectors might by important to place constraints on the EOS at densities and temperatures exceeding those which are probed during the inspiral [444].
The postmerger waveform shows a non-monotonic amplitude and frequency evolution; cf. Fig. 1. The dominant and most robust feature of the postmerger GW signal is the main emission frequency, called f 2 or f peak [54,131,501,63,56,132,463,521]. In addition also a number of other, sub-dominant frequencies are present, e.g., Refs. [132, 56,463]. In general, the frequency of the GW signal emitted during the postmerger phase is determined by the rotational state of the remnant and the EOS. Based on a large set of NR simulations, a number of groups have proposed different quasi-universal relations, which, in principle, can be used to determine the part of the EOS governing the remnant's density evolution [444]. Recently, Refs. [121, 521, 96] employed a full Bayesian framework to investigate which SNR is required to extract information about the EOS from the postmerger GW signal and found that a minimum postmerger SNR of 3 to 5 is required. This corresponds to a total SNR of 100 to 200, which would allow very precise measurements about tidal effects from the inspiral GW signal.

Numerical Relativity Simulations
NR simulations are a fundamental tool to study GWs from compact binaries in the late-inspiral and postmerger phase of the coalescence. The main advantage of NR simulations is that the Einstein Equations can be solved directly by numerical discretization of the equations. Over the last years, the NR community has made tremendous progress on several fronts: (i) the exploration of the BNS parameter space by systematically varying the EOSs, total mass, mass ratio, spin, and eccentricity, e.g., [275,65,225,70,182,333,480,64,522,288,Fig. 5 Overview about the possible postmerger dynamics and the possible postmerger outcome. Mostly depending on the mass (but also the mass ratio and spin), the merger remnant can promptly collapse to a BH, otherwise the remnant can be a HMNS or SMNS collapsing on a longer timescale to a BH or, if the mass is low enough, it can form a stable NS. 176,191,287,179,499,173,241,193,190,446,420,192,524,384,468,295,124]; (ii) the development of numerical schemes allowing many-orbits simulations with high-accuracy GWs, e.g., [71,277,448,62,274,174,294]; (iii) the increasing realism of microphysical schemes that incorporate magnetic effects, finite-temperature EOSs, composition effects, and neutrino transport, e.g., [26,239,461,405,297,406,467,129,296,478,234,400,479,406,228,223,489,227,223,443,488,487,486,230,383]; (iv) the study of mass ejecta and EM counterparts, e.g., [273,53,538,333,446,182,176,228,231,427,271,95].
In the following, we will focus on the computation of the GW signal and therefore, purely on GRHD simulations. We refer the interested reader to, e.g., [38,35,460,464], for discussions about effects of magnetic fields or neutrino radiation. GW signals for a variety of different binary parameters: a high total mass setup, a high mass-ratio configuration, a precessing simulation, and a highly eccentricity configuration. The simulations are taken from the CoRe database and labeled accordingly. Adapted from [180].

Imprint of the binary properties and waveform catalogs
Waveform Data Bank 2 [295]. Both catalogs combine several hundred simulations and are freely available. In addition large simulation campaigns have been performed by other groups, e.g., Refs. [272,502,53]. All these simulations allow the NR community to obtain a better understanding of the imprint of the individual binary properties on the merger process. In the following, we want to discuss shortly the imprint of the total mass, mass ratio, spin, and eccentricity.

Total mass
In contrast to BBHs, BNS systems are not scale invariant and therefore can not be rescaled by the total mass since the mass enters in the computation of tidal interactions during the inspiral and determines the merger outcome; cf. Fig. 5. Based on the constraints derived from GW170817 one expects that the mass of an individual NS lies within M ∈ [1.0, 2.3M ] [327,403,462,486,467,361,491].
As discussed in Sec. 2, systems with large total masses are typically characterized by a possible prompt collapse directly after the moment of merger. Such a system will be very similar to a BBH mergers, see e.g. [251].
In addition to this prominent postmerger evolution, which is more similar to the BBH dynamics (top panel of Fig. 6) than a typical BNS merger, massive NSs will have small tidal deformabilities and consequently tidal effects are for massive BNS mergers less important than for lower-mass systems; cf. Fig. 4. Consequently a distinction between high mass BNS and low mass BBH systems will be problematic. This problem becomes even more pronounced since massive BNS mergers eject only a small amount of material and, thus, create no bright EM counterpart. Therefore, the correct astrophysical classification is complicated and is an active field of research, see e.g. [125].

Mass ratio
Based on population synthesis models a maximum mass ratio [q = M A M B ≥ 1] of q max 1.8−1.9 [186,179] might be possible for a BNS system. While observation of isolated NSs distribution would theoretically allow maximum mass ratios of q max 2, the largest observed mass ratio in BNSs is q ∼ 1.3 [364,330]. Over the last years, the NR community managed to simulate systems with large q, at the boundary of physically realistic scenarios, see [179,517] for simulations with q 2.
With increasing q, one finds that the emitted GW energy decreases during both the inspiral and the postmerger phase; cf. Fig. 14 of [182] and that the merger frequency decreases roughly as a function of ∝ q −1/2 , cf. [183] 3 .
Considering the post-merger evolution, one finds that, while the EOS determines most of the evolution, systems with larger mass ratios tend to have merger remnants which seem to survive slightly longer before gravitational collapse and BH formation. In combination with the longer lifetime, the mass of the formed accretion disk increases with increasing q [182,295,177]. Due to the large disk and ejecta mass of unequal mass systems one expects them to have brighter EM counterparts. In addition, one finds that systems in which the mass ratio is large enough to tidally disrupt the secondary star before the merger, the postmerger GW amplitude is significantly reduced, cf. Fig. 6 (second panel) and [521].

Spin
The maximum NS spin is not precisely known, since it depends on the unknown EOS, but EOS models predict breakup spins of up to χ ∼ 0.7, which corresponds to spin periods of less than 1 ms [348].
Observationally, a large number of highly-spinning NSs are known. The pulsar PSR J1748−2446ad [264] is, to date, the fastest spinning NS with a frequency of 716 Hz. PSR J1807-2500B [354] is the fastest spinning NS in a binary system with a rotation frequency of 239 Hz and PSR J1946+2052 [498] is the fastest spinning NS in a BNS system (59 Hz). Due to magnetic dipole radiation, J1946+2052 will spin down and will presumably have a dimensionless spin of about χ ∼ 0.05 at merger.
In the past, consistently evolved spinning BNS simulations have been presented in Refs. [64,179,176,172,384,524,192] and precessing NR simulations are presented only in [179,499,173,124].
With respect of the dynamical evolution and energetics of spinning NS binaries during the inspiral one finds that even for astrophysical realistic spin magnitudes (χ 0.05) spin-orbit interactions are larger than tidal effects up to or shortly before the contact of the two stars [64]. Therefore, neglecting spin effects might lead to a wrong description of the binary evolution if one of the two binary constituents is a pulsar, which would introduce biases during the determination of tidal parameters from GW observations, e.g. [208,17,471].
Studying the energetics during the postmerger, one finds that spin aligned systems can create merger remnants with larger angular momentum than nonspinning systems. Consequently, collapse to a BH happens at later times [124]. We note that simulations also indicate that characteristic frequencies in the postmerger, which generally also follow quasi-universal relations [57,56,502,501,63,132,463], might shift under the presence of spin.

Large Eccentricities
While most NS binaries are expected to be formed out of an already existing binary system, there is a chance of producing BNS systems by dynamical capture in dense stellar environments such as globular clusters [332]. Those binaries can have a non-negligible eccentricity when they enter the LIGO frequency band. Unfortunately, estimates for the rates of eccentric BNS mergers are highly uncertain. Ref. [519] gives a conservative estimate of a BNS merger rate of 0.5yr −1 Gpc −3 .
While eccentric BNS mergers will occur infrequently, the detection of at least one of these events would allow for an intensive study of the strongfield regime and the properties of NS matter. Of particular interest for highly eccentric systems is the possibility to constrain the EOS during the inspiral from density oscillations induced by close encounters of the two stars. In fact, when the NSs approach each other, the induced tidal deformation of the star causes a mode excitation. Once excited, most energy is released through the f -mode. Consequently a measurement of this excitation frequency would allow to place additional constraints on the unknown supranuclear EOS, cf. e.g. [525,299,241,379,126,424,552,123]. Furthermore, since information of the f -mode allows to place constraints on the binary properties, an accurate description also allows to understand possible resonance effects when the orbital motion approaches the f -mode frequency and due to a resonances [497,476,28].

High-quality NR simulation
As pointed out, there have been more than a thousand individual BNS simulations performed over the last decades, however, up to now, there is only a limited number of simulations which meet the accuracy standards necessary

D/M
Iter0-e = +0.0000,ṙ = +0.0000 to develop and test existing GW models; see e.g. [449,62,274,523,174,294,181,383,295] for further discussion. The main difficulty for the production of high-quality waveforms are: (i) the reduction of artificial eccentricity; (ii) the usage of high enough resolutions to reduce discretization errors; (iii) the development of numerical schemes showing clear convergence across multiple resolutions and configurations.

Eccentricity reduction
As an exemplary case addressing all listed items, we present a high-resolution simulation of [181]. The setup that we pick has been simulated for six different resolutions and describes an equal mass, non-spinning configuration with NS masses of M A = M B = 1.35 employing a piecewise-polytropic fit of the SLy EOS [454,179]. The eccentricity reduced initial configurations are obtained with the SGRID code [514,515,516,179,517]. SGRID is one of the few codes able to produce data with arbitrary eccentricities. Up to our knowledge, other codes capable of computing low-eccentric data are SpEC's spectral solver SPELLs [226,499] and the non-public extension of LORENE [2] discussed in Ref. [310]. The computation of low-eccentric data relies on an iterative procedure in which the binary's initial radial velocity and the eccentricity parameter (or the orbital frequency) are varied such that during an evolution the eccentricity reaches a certain threshold, see [179,517]. Figure 7 shows the eccentricity reduction procedure for this configuration. In most cases it is sufficient to perform 2 to 4 iteration steps to achieve final eccentricities 10 −3 .

Errors through numerical discretization
The configurations which we want to discuss in the following have been simulated with the BAM code [99,508,175,62] for six different resolutions covering the individual NSs with 64,96,128,192,256, and 320 points respectively.
The highest resolution (320 points in the finest refinement level) has a spatial resolution of 0.047M ≈ 70 m in the finest refinement boxes. The computational cost for all six resolutions is about ∼ 8 million CPU-hours, where as for all 3+1-NR codes the computational expenses increase with the fourth power of the number of grid points, i.e., the 320 point simulation requires more than 800 times the resources compared to the lowest resolution setup.
The real part of the GW signal for all resolutions is shown in the top panel of Fig. 8. It is clearly visible that the merger time increases with an increasing resolution. This behavior is understandable since numerical dissipation, which decreases with an increasing resolution, adds dissipation and accelerates the inspiral. Therefore, the merger for an 'infinite resolution' would be later than for the highest resolution shown in Fig. 8. Because of the presence of a clean convergence order, as visible in the middle panel, one can use a Richardson extrapolation to obtain a better guess for such an 'infinite resolution' scenario.
In the bottom panel we show among others the phase difference between the Richardson extrapolation using the highest resolutions, R(n 320 , n 256 ), and the phase of the highest resolution n 320 .
In addition, one can check the robustness of the extrapolation by computing the phase differences R(n 320 , n 256 )−n 320 and R(n 256 , n 192 )−n 256 in the bottom panel of Fig. 8. Rescaling the phase difference of R(n 320 , n 256 ) − n 320 assuming second order convergence shows again excellent agreement with R(n 256 , n 192 )− n 256 . This proves that the leading error term scales quadratically with respect to the resolution.
At merger, the phase difference R(n 320 , n 256 ) − n 320 is only 0.37 rad and provides a good error estimate of our simulation. An alternative, but non-conservative error can be approximated by the difference between two Richardson extrapolated waveforms (green line in the bottom panel). Throughout the whole inspiral, the difference between R(n 320 , n 256 ) and R(n 256 , n 192 ) is below 0.1 rad. However, due to the non-monotonicity of R(n 320 , n 256 ) − R(n 256 , n 192 ) it is unclear to what extent higher-order error terms are properly modeled.

Finite radius extraction
In addition to the discretization errors, one has to consider for the full error budget of a simulation that the GWs are extracted at a finite extraction radius, while in principle GWs should be extracted at future null infinity. We refer the interested reader to [84] and references therein.
The best solution to obtain properly extracted results is the Cauchy Characteristic Extraction (CCE) [83,459], see also [34,254].
Another option and computationally cheaper approximation is the use of a polynomial extrapolation. For this purpose, the GWs are extracted at different radii r i with i = 0...N and then the phase and amplitude are extrapolated using a polynomial of the order K with K < N : Using this extrapolation method and estimating the error as the difference between the extrapolated and non-extrapolated waveform, the finite radius extraction error is below 0.044 rad for our example shown in Fig. 8. Another possibility is it to consider the next-to-leading-order (NLO) falloff behavior of the Ψ 4 multipoles [352,62], and obtains an extrapolation formula for rḧ m which uses ψ m extracted at a given radius. It is interesting to note that, as shown explicitly in Fig. 12 of [62], the finite extraction radius decreases for an increasing GW frequency. Thus, it is smaller at the moment of merger than during the early-part of the inspiral.

Post-Newtonian Waveform modelling
The PN modeling of GWs from NS binary systems has significantly advanced in recent years. The description of spinning point-masses in PN theory, which applies for black holes up to the absorption of GWs by their horizons, have been iterated to high order. Progress has also been made on elucidating several subtleties that arise in higher-order calculations and on developing appropriate regularization schemes within different approaches for computing the dynamics. We refer the reader to the review article [86] for a comprehensive list of references on this decades-long effort; for other reviews see [232] and [475]. To date, the dynamics have been computed to 4PN order with three classes of PN methods: a Fokker-action approach, with the most recent new results obtained in [359,60,61], a Hamiltonian approach [150,151], and effective field theory methods [216,217,335,338,91]. Other approaches including post-Minkowski expansions and further applications of methods from high-energy theory have garnered much recent interest and seen rapid developments that already enabled computing higher order results, e.g. [78], but will not be discussed here. For high-order computations of the gravitational radiation, however, only results from a single method, the multipolar post-Minkowski expansion [87,86,510], are currently available. Other approaches have only been carried out up to second order, such as the direct integration of the relaxed Einstein equations [90,89,198] and effective field theory methods [334] based on tools developed in [237,356,236,235,82,81,335,438,465,218,243,242].
In this section, we will focus the discussion on modeling the additional contributions from the presence of matter on top of the point-mass descriptions. We will give details only for a few dominant effects that are included in current waveform models used in data analysis, specifically only the quadrupole and octupole tidal effects from the fundamental modes and their adiabatic limits. As outlined in the Introduction, analytical modeling of the binary inspiral is based on a tapestry of approximation schemes for the dynamical spacetime describing the binary system. These perturbative results can then be re-summed and combined with additional information from test-particle limits and NR into EOB models as will be described in the next section.
The starting point for approximations is to divide the problem into GW generation and propagation. In the distant wave zone, at distances from the binary much larger than the wavelength λ of the emitted GWs as well as the binary's size and gravitational radius, and outside the regions of significant spacetime curvature produced by the binary, the perturbation expansion is a post-Minkowski multipolar scheme that yields a solution for the radiation field [86]. The source of that radiation field considered here is a nonspinning binary system with comparable masses on quasi-circular orbits in the regime where the orbital separation between the bodies r is large compared to their characteristic size R. The orbital dynamics can be approximated by a PN expansion with dimensionless parameter ∼ GM/rc 2 ∼ v 2 /c 2 O(1) and is dominated by point-mass contributions. Small corrections from finite-size effects are included through a tidal expansion in the parameter α = R/r. This yields a double expansion where and α are treated as distinct parameters that correspond to different physical effects, although they scale in the same way with the orbital separation [213]. As we will discuss below, the finite size effects involve characteristic matter parameters such as tidal deformabilities that enter as coupling constants in the effective action for the dynamics and influence the radiative multipole moments of the binary and hence are imprinted in the GWs [214,266]. The link between these constants and the microphysics of matter are determined by considering the physics in the body zone, close to each NS, where the self-gravity of the NS is strong but the influence from the distant companion is a small perturbation. Solving the linearized Einstein field equations for perturbations to the NS's equilibrium configuration yields the perturbed NS matter distribution and spacetime geometry, and enables extracting the knowledge of the characteristic matter coefficients.
In this framework, Newtonian quadrupolar tidal effects in the dynamics first appear at O(α 5 ). Relativistic corrections have been computed up to O(α 5 2 ) for the quadrupole and O(α 7 2 ) for the octupole [263,77]. The knowledge of effects in the GW signals will be described in more detail below. We also note that the scalings discussed above apply only for adiabatic tidal effects, and in particular, the resonant excitation of a generic mode does not fit into that counting scheme.
For simplicity we first consider a NS with a point-mass companion. To linear order in the finite-size effects, the case of two extended objects can be recovered by adding the same contribution with the body labels interchanged. We will assume that the bodies are nonspinning and spherically symmetric in isolation. To motivate the relativistic calculations of multipole moments and their impact on the dynamics and GWs, we will first discuss the Newtonian description expressed in a form that carries over to General Relativity (GR) with appropriate modifications. We will specialize to the same assumptions made in models that are currently well-developed for data analysis, and limit the details to the dominant mass-quadrupolar tidal effects for simplicity; extensions to higher mass multipole moments are straightforward. At the end of this subsection we will discuss various other effects that have been studied and are expected to give rise to interesting subdominant contributions.

Multipole expansion: NS's mass moments and tidal moments
Below, we will first introduce the multipole moments generated by the NS itself, then consider the multipolar tidal moments due to the companion that are felt by the NS.

NS's mass multipole moments in Newtonian gravity
In Newtonian gravity the gravitational potential U generated by a mass density ρ at a field point x is a solution to Poisson's equation ∇ 2 U = −4πρ. The Green function solution is: For points x > x outside the mass distribution we can perform a Taylor series expansion around a moving reference point z(t) to write Eq. (6) as Here, L = a 1 a 2 · · · a denotes a string of indices, x L = x a1 · · · x a , and the summation over repeated indices is implied. Likewise, the notation for deriva- The expansion of the potential (7) can be written in a more suggestive form by defining the NS's Newtonian mass multipole moments by the following integrals where M Newt is the mass monopole and Q L the -th multipole moment. The angular brackets around indices denote the symmetric and trace-free projection of a tensor. For example, x <ij> = x i x j − δ ij |x| 2 , see Refs. [510,259] for more details. We have defined the moments as their trace-free parts because the contribution to Eq. (7) involves the derivative ∂ L |x − x | −1 which projects out the trace-free piece Q L . Because our reference point z i (t) was chosen as the NS's center of mass the dipole term ( = 1) is absent from the expansion (8).
With these definitions the exterior potential takes the form This expression in terms of Cartesian tensors is useful for computing the dynamics of a binary system. When calculating the deformation of the NS, it is more convenient to use spherical coordinates (x − z) i /|x − z| = (sin θ cos φ, sin θ sin φ, cos θ) and express (9) as an equivalent sphericalharmonic expansion where the relation between spherical and Cartesian multipole moments is [510]: The tensors Y m L are defined by Y m = Y m L n <L> and involve complex coefficients that convert between the unit vector n and spherical harmonic representations [510].

Multipole moments of the NS's spacetime in GR
The Newtonian expansion around the center-of-mass position can be generalized to an expansion around a reference center-of-mass worldline z µ (σ) where σ is an evolution parameter along the worldline. This is known as the worldlineskeleton approach [184], and also arises naturally in frameworks based on an effective action [335]. The NS's multipole moments are defined in a region of spacetime at large distance compared to the size of the NS. In this region, and when expressing the asymptotic metric in a local asymptotic frame [510,511], the time-time component of the metric involves a potential U eff that is analogous to the Newtonian gravitational potential U : For a spherical object U eff = M/r where r is the distance from the body, while for a nonspherical, nonspinning body it has the asymptotic form In this setting, the -th mass multipole moment is associated with the coefficient of the term that falls off as r −( +1) . This definition is equivalent to the Geroch-Hansen multipole moments for stationary spacetimes [249]. We note that the multipole moments can be defined by considering the energetics of test particle orbits [469], and that in a binary system a small ambiguity remains in the definition of the multipole moments [513,246].

Tidal moments in Newtonian gravity
Next, we consider the gravitational potential due to the binary companion that is felt by the NSs. We will denote it by U ext to indicate that the source of this potential is external to the NS. A Taylor expansion around the NS's center of mass leads to: The coefficients E L in (14) are the NS's tidal moments [513]. We are interested in specializing the general expansion (14) to the case that U ext is sourced by a binary companion whose gravitational potential we denote by U c . The NS's tidal moments are then given by

Tidal moments in GR
The relativistic tidal moments are projections of the curvature tensor due to the companion. The analogue to the Newtonian E L is the 'electric part' of the companion's Weyl tensor C µανβ as [513] where u µ = dz µ /dσ is the tangent to the NS's worldline and we have defined The tidal tensor E µν is symmetric and trace free. In the NS's rest frame u µ = (−1, 0, 0, 0) it is a purely spatial tensor E µν u ν = 0 and given by E ij = C titj . Higher multipole tidal moments are defined by where the semicolon denotes covariant derivatives.

Equations of motion and action principle
We next consider the effects in the interaction zone, where the dynamics are dominated by the point-mass motions with small corrections from (i) the coupling of the NS's multipole moments to the tidal field, and (ii) the energetics associated with the internal dynamics of the multipoles.

Tidal effects on Newtonian binary dynamics
The center-of-mass position of a NS in a binary system moves according to Newton's second law: Note that only the potential sourced by the companion contributes to the body's motion, as can be verified by direct calculation. The dynamics can be concisely summarized by an action principle. The action is constructed from the Lagrangian L = T − V , where T and V are the kinetic and potential energies for the system. These can be separated into contributions from the NS's center-of-mass motion and from its internal dynamics [535]: Expanding around the NS's center of mass, adding the contributions from the companion, and transforming to the barycentric frame of the binary system leads to Here, M and µ are the total and reduced mass respectively. The relative separation is r = z − z c with magnitude r = |r|, and the relative velocity is v 2 =ṙ ·ṙ. The action for the binary dynamics is then given by [214,452,316] where S pm = dt (µ/2)v 2 + µM/r describes the orbital motion of pointmasses and L int the internal dynamics of the multipole moments. We will assume that the multipoles arise only from the response to the companion's tidal field [74,458,317,316,458,557,556,301,298,255,381,299,214,210,164,485], which excites the oscillation modes of the NS. We will consider here only the fundamental (f -) modes. They behave as harmonic oscillators described by a Lagrangian [214,452,316,299,116] Here, the quantities ω 0 denote the f -mode frequencies and we have neglected contributions from other quadrupolar = 2 modes. The parameters λ are the tidal deformability coefficients. They are defined by considering the adiabatic limit, where the NS's internal time scales τ int ∼ ω −1 0 ∼ R 3 /m are fast compared to the time scale of variations in the tidal field τ orb ∼ r 3 /M as The tidal parameters λ are related to the NS's dimensionless tidal Love numbers k [353] by λ = 2 (2 −1)!! k R 2 +1 and depend on the NS's interior structure. For data analysis applications, it is useful to define the dimensionless tidal deformabilities Λ = λ m 2 +1 .
For adiabatically induced multipoles, dQ L /dt = 0, the internal Lagrangian is only the elastic potential energy associated with the deformation L int adiab = −Q L Q L /(2 !λ ). In this limit the mass quadrupole finite size effects on the dynamics depend only on the orbital variables and λ : with E L = −m c ∂ L r −1 .

Tidal effects on relativistic dynamics
The action describing tidal interactions in a relativistic binary system can be obtained by expressing the Newtonian result from Eq. (21) in a covariant form, inserting the appropriate redshift factors to ensure invariance under reparametrizations, and replacing all time derivatives by covariant derivatives along the center-of-mass worldline [497]: Here, total derivatives are denoted by D dσ = u β ∇ β , where ∇ α is the covariant derivative. Contributions from quadrupolar modes with higher radial nodes and other terms due to the incompleteness of quasi-normal modes have been omitted. The interaction terms in the action (26) also follow from an effectivefield-theoretical approach [244,197,497] when considering all possible terms consistent with the symmetries (general covariance, parity, and time reversal) and re-defining variables to eliminate accelerations. After decomposing Eq. (26) into the time and space components and imposing the constraints to eliminate all unphysical gauge degrees of freedom we obtain Here, the term L FD describes relativistic frame-dragging effects that arise kinematically when expressing (DQ µν /dσ) 2 in terms of (dQ ij /dt) 2 . The term L FD describes the coupling of the orbital angular momentum to the angular momentum (or spin) associated with the quadrupole To proceed with computing the dynamics from (27) and the GWs requires explicit expressions for the various quantities such as E ij and z and the framedragging terms appearing in the action. These have been computed in post-Newtonian (PN) theory [535,77,497], see also [335] for the effective field theory description, in the test-particle limit [79] or the gravitational self-force formalism [185,75,402,482]. An alternative approach based on an affine model has also been studied [210,367].

Computation of tidal deformability
We have seen that the action summarizing the dynamics involves the characteristic information about matter properties through parameters tidal deformability and, in the non-adiabatic case, the mode frequencies. The mode frequencies are routinely computed in calculations of stellar oscillations in Newtonian gravity and of quasi-normal modes in GR [300,209]. Thus we will focus here on the tidal deformability parameters defined in terms of spacetime multipole moments and their relation to the tidally perturbed interior structure of the NS.

Newtonian calculation of tidal deformability
We start from an equilibrium configuration for the NS in isolation. It is a solution to the Poisson equation for the gravitational potential together with conservation of mass and momentum: Here, a i ext is the acceleration due to external forces which vanishes for a star in isolation. Solving the system of equations, (29), for a given EOS model yields a background configuration of an isolated NS with pressure p 0 , density ρ 0 , and gravitational potential U 0 . Next, we consider the effect of external static tidal field E L . This causes the NS matter to adjust to a new static configuration with an exterior gravitational potential given by The linear response defining λ in (23) can equivalently be written as for each -th multipole and azimuthal quantum number m. To obtain λ , we need to consider only a single value of m, and it is easiest to choose m = 0.
We will need to compute the perturbed interior solution for the gravitational potential and matching it to the exterior description (30) at the NS surface. The equations (29) also describe the perturbed interior, with p = p 0 + δp, ρ = ρ 0 + δρ, U = U 0 + δU , and an external acceleration a ext = ∇U tidal . The fluid perturbation can be represented by a Lagrangian displacement ξ(x, t), which maps a fluid element at position x in the unperturbed NS to x + ξ(x, t) in the perturbed star. Expanding the Euler equation to linear order in the perturbations leads to the differential equation This expression is used in general for calculating stellar oscillations. Here we are interested in specializing to static perturbations with (ξ = 0). For a barotropic EOS relation p = p(ρ) we can eliminate δp from (32) in favor of δρ, combine all terms that involve δρ into a total derivative and integrate. We then decompose the perturbations as δρ = f (r)Y m (θ, φ) and δU tot = H(r)Y m (θ, φ), substitute into the linearized Poisson's equation, and after manipulations obtain the differential equation for H(r) in the region r ≤ R: Except for special choices of the EOS, solving (33) requires numerically integrating in the NS interior, with the boundary condition that ensures regularity at the center of the star, H ∝ r for r → 0. For r > R the exterior solution is To extract the Love numbers we eliminate E m by considering the logarithmic derivative We solve for k by using (34) in (35) and matching the result for y(R) obtained from the interior and exterior solutions at the star's surface: For an incompressible n = 0 polytrope the exact solution is k = 3/(4 − 4).

Calculation of tidal deformability in GR
The definition of λ from Eq. (23) is general and holds with the relativistic definitions of E L and Q L . The computation of λ in GR follows a similar method as the Newtonian case but requires replacing the gravitational potential by the spacetime metric, the Poisson equation by the Einstein field equations, the density by the stress-energy tensor, and the continuity and Euler equations by the relativistic energy-momentum conservation. This complicates the equations so we will only outline the process here.
The NS matter is described by the stress-energy tensor T µν = (ρ + p)u µ u ν + pg µν , where u µ is the fluid's four-velocity and ρ the energy density. In the NS's rest frame, the normalization condition u µ u µ = −1 implies that u µ = (e −ν/2 , 0, 0, 0). The Einstein field equations and energy-momentum conservation ∇ ν T µν = 0 determine the NS structure This system reduces to the Tolman-Oppenheimer-Volkoff equations that in general must be solved numerically except in special cases. When considering the tidal perturbations of interest here, the metric is the background plus the linear perturbations g µν = g (0) µν + h µν dx µ dx ν . For our purposes we can decompose the perturbation as [455,512,279,170] The initial condition at the center is that H ∝ r to ensure regularity of the solution.
Outside the star, the metric perturbation has the general form H = a Q Q 2 (x) + a P P 2 (x), where x ≡ r/M − 1 and P 2 (x) and Q 2 (x) are the normalized associated Legendre functions of the first and second kinds respectively. They are normalized such that for x → ∞, P 2 (x) ∼ x and Q 2 (x) ∼ x −( +1) . The constants a P and a Q are determined by matching the logarithmic derivative of the interior and exterior solutions, at the NS surface [266]. Comparing with the definition of Q and E in the asymptotic metric (13) and the definition of the Love numbers shows that For the dominant quadrupolar effect the explicit expression is where y is evaluated at the surface r = R. For an incompressible star with ρ = const or p = Kρ 1+1/n with n = 0, the density profile is a step function and the matching of the interior and exterior solutions must be modified as follows [154]. From Eq. (41) in the interior we obtain y in (R). The density discontinuity introduces a correction to y just outside the star y out given by  166,346,347,112]. Tidal parameters can also distinguish exotic objects, matter halos, and modified gravity [108,481,374,526,412,370,196,401,49,195,117,369,318,165,240,441,107,483,413,188,250,553,114,109,59]. Substantial recent interest has quasi-universal relations [546,547] that link dimensionless parameters characterizing various global properties of the NS in an approximately EOS-independent way [328,423,545,422,421,260,115,365,24,127,417,496,260,97,493,543,457,118,119,122,308,472].

Other matter effects
In the discussion above, we specialized to tidal effects for a nonspinning NS from its fundamental modes. A NS has a rich spectrum of modes. Those with sufficiently low frequencies could be tidally excited during an inspiral and yield spectroscopic information about matter in NS interiors; several examples have been studied to date [440,27,269,215,485,554,316,299,519,520,431]. Effects from nonlinear mode coupling have also been studied [542,199,322]. In GR, besides the tidal couplings to the gravito-electric fields E L that induce mass multipole moments, there are also gravitomagnetic tidal fields that induce current multipoles. These effects have no Newtonian analogue but their mathematical description is similar to the gravitoelectric case, although the interpretation and adiabatic limits are more subtle. In the local frame of the NS, the gravitomagnetic part of the curvature due to the companion is given by B L = 3 <a1jk C jk a20;a3···a > /2( + 1)( − 2)! where ijk is the completely antisymmetric permutation tensor. The induced current moments S L appear in the time-space part of the asymptotic metric in a local asymptotic frame, see Refs. [207,323,414,322,215] for more details. Similar to the tidal deformability coefficients λ , the relation between S L and B L is characterized by a set of gravitomagnetic Love numbers σ = −S L /B L .
When considering rotating NSs, the additional spin degrees of freedom introduce a yet richer phenomenology. Spins give rise to new spin-tidal couplings [283,15,415,416,233,320,324,181] and can significantly impact the tidal excitation of quasi-normal modes [187,269,229,355]. As discussed in the previous section, eccentricity also greatly enriches matter phenomena and mode excitations [241,126,551,531]. While this review focuses on double neutron star systems, we mention here for completeness that for NS-BH binaries there exist parameter regimes of mass ratios, spins, and NS compactness where the NS is tidally disrupted. This leads to a sudden shutoff of the GW signal, with the disruption frequency depending on the parameters including the EOS [527,490,419,211,366,224,291,313], see Ref. [490] for a review. An additional distinction between NSs and black holes is that the neutron star has a surface while a black hole possesses an event horizon that absorbs all incoming GWs; this effect is also imprinted in the GWs [434,507,500,368,433,202,435,388,69,504,280].

Matter signatures in GW signals
The asymptotic GWs are computed by combining an expansion using powers of G as the bookkeeping parameters around flat space with a multipolar decomposition [86]. The radiative post-Minkowski solution can be formally written in a radiative coordinate system (u, X) where u = T − |X|/c is the retarded time, as Here, P ijab (N ) is a transverse-traceless projection operator, N is the vector pointing between the source and the observer, D is the distance to the source, and U L and V L are radiative mass and current multipole moments. Their relation to the multipole moments of the binary system is nonlinear and nonlocal in time, and complicated in general. At the leading Newtonian order, there is a simple relation to the source quadrupole U Newt ij = d 2 Q T ij /du 2 , where the total quadrupole moment of the system is the sum of the orbital and NS contributions Q T ij = µr 2 n <ij> + Q ij . In Eq. (45) we have displayed explicit factors of c because they indicate the PN order, with each factors of c −2 being considered one PN order. Using (45), the energy and angular momentum fluxes in GWs can also be expressed in terms of the moments U L and V L [86].
An approximation for the GW phase evolution can be computed by imposing that the power radiated by a binary system is balanced by a change in the energy of the binary. We will briefly review this approach to obtain the lowestorder matter signatures in GWs. The energy of the binary is E = E pm + E tidal , where the subscript "pm" denotes point-masses. Specializing to circular orbits r =ṙ = 0 andφ = Ω,φ = 0, and to perturbatively solving the dynamics described by the action (25), and eliminating r in terms of the frequency variable x = (M Ω) 2/3 we find the following tidal contribution in the limit of adiabatic tides [535] Here, the subscripts A, B label the two objects, and the factor outside the brackets is the result for Newtonian point masses. The leading order adiabatic tidal corrections contributions to the energy flux from the mass quadrupole radiation are [533] By requiring that P GW be balanced by a change in the energy E of the binary one can derive the evolution equations There are several ways to solve for φ in a PN approximation. For example, one can numerically solve Eqs. (48) for φ(t) and x(t) after first expanding the ratio P GW /(dE/dx) about x = 0 to the consistent PN order. These waveforms are known as TaylorT4 approximants as reviewed in Ref. [103], where the point-mass terms are also given explicitly. The adiabatic quadrupolar tidal corrections that add linearly to the point-mass contributions: Other possibilities to perturbatively solve (48) and obtain the tidal contributions to different approximants for the gravitational waveform are detailed in the Appendix of Ref. [537]. Another widely utilized class of template waveforms for data analysis are TaylorF2 waveforms which provide a fully analytic frequency-domain model and are thus very fast to generate. The derivation is explained e.g. in Ref. [140] and leads to a signal of the form where f is the GW frequency, A ∝ M 5/6 /D, where M is the chirp mass M = η 3/5 M , and D is the distance between the GW detector and the binary. Extrinsic parameters of the source such as the location on the sky are also contained in A, where higher PN order corrections to the amplitude enter. The point-mass phase ψ pm to the current best knowledge for nonspinning binaries is given e.g. in Eq.(3.18) of [103]. The leading order tidal contributions to ψ tidal can be expressed as: x 5 − 4000 9 The parameterΛ plays an analogous role as the chirp mass. For equal-mass BNSsΛ reduces to Λ of the individual neutron stars. Other combinations of the two parameters Λ A and Λ B are also in use and have advantages in different contexts, e.g. to characterize the dominant effect in the conservative dynamics [155], or to improve the measurability [548,549]. This shows that for a double neutron-star system GW measurements are most sensitive to a weighted average of the deformability parameters of the two objects, which complicates the interpretation of measurements.
The effects of the fundamental mode excitations can be approximately included in this waveform model by adding to (51a) the following modefrequency-dependent and regularized contributions [476]: Finally, we give the full analytical knowledge to the highest complete PN order O(α 5 3/2 ), i.e., 1.5PN corrections to tidal effects for the quadrupolar adiabatic phasing [163]: with the dimensionless EOB tidal parameter κ A (defined below) and x(M ω) = (M ω/2) 2/3 . The individual coefficients c A i are and similarly with A ↔ B. Partial knowledge is available at higher orders. The O( 5/2 ) terms are known, however, the lower order O( 2 ) contributions have only been computed for the dynamics but not the fluxes. As all pieces of known contributions are used in the phenomenological models discussed below, we also give the two additional terms that contribute with a corresponding power of x to (53): The effects of spin-multipole moments also couple to the companion's tidal tensors and influence the GWs in an analogous way as tidally-induced multipole moments except the scalings with parameters are very different. For the most recent developments on modeling the effects of spin-induced multipole moments on the GW signals see Refs. [306,92,363,337,438,437,336,102]. The leading-order contribution to the GW phase from the spin-induced quadrupole is given by [306]: Figure 9 illustrates the various contributions to the GW phase discussed above.
From the discussion above, it is apparent that the fractional corrections to the Newtonian point-mass results due to tidal effects scale as a high power of the frequency, x 5 and higher, where x = (πM f ) 2/3 , which corresponds to O(α 5 ) in the expansion discussed in the beginning of this section. This means that in a loose PN counting that is based on assigning a PN order to each power of x, tidal effects first enter effectively as 5PN corrections. As explained above, physically and formally the PN corrections to the point-mass dynamics and GWs and finite size effects are distinct. Nevertheless, the net effect of unknown PN point-mass terms can impact measurements of matter effects, as will be discussed in more detail in Sec. 7. There are two classes of effective or phenomenological models for black hole binaries that effectively include all PN orders for the point-mass dynamics in an approximate way. These are the effective one body (EOB) model [100,101] and the so-called "Phenom" models [19,20], both of which aim to combine the available information on the relativistic two-body problem from different regimes into a single framework to generate waveforms for data analysis as we will discuss next.

Tidal Effective-One-Body Models for binary inspirals
The EOB approach [100,101] (see also the reviews [159,156,104,143]) combines results from the PN theory, valid for arbitrary mass ratios, with information about the strong-field effects in the test-particle limit. The EOB dynamics involves a nonlinear mapping of the phase space of the binary system to that of an effective particle on a nearly geodesic trajectory in an effective spacetime. The effective particle's mass is the reduced mass while the mass of the spacetime is the total mass of the binary system. This setting reduces to testmass motion in a BH spacetime in the test-mass limit. For finite mass ratio, corrections are added to the metric potential that depend on the symmetric mass ratio ν ∈ [0, 1/4] and thus deform the spacetime away from those of a BH. In addition, non-geodesic corrections that depend on higher-thanquadratic powers of momentum appear. These modification are dictated by the requirement to reproduce the known PN results in the weak-field limit. This theoretical structure achieves an extremely concise resummation of PN information into a small number of terms. There is still considerable remaining freedom in the model, for instance on choosing the functional form for the potentials which are generally expressed as non-analytic functions. These choices are constrained by considerations such as a pathology-free physical behavior of quantities such as energetics of the system under variations of parameters. In addition, extra terms that parameterize unknown higher-order information are added to the model that are fixed by comparing to results from NR and gravitational self-force [43,44,76,32,31,42].
The dissipative sector of current EOB models is based on a concise factorized representation of the gravitational waveforms [389,152,148,376,506,505,410,394,68,148,408] that accurately accounts for the modifications to the GW propagation from the spacetime curvature due to the total mass of the binary. These GW modes are explicit algebraic expressions that depend on the instantaneous EOB dynamics. The evolution of the orbits is computed by imposing energy and angular momentum balance. Calculating a full waveform therefore requires solving the time-domain coupled system of the Hamiltonian equations, the GWs, and the corresponding fluxes that give dissipative contributions to the equations of motion over an inspiral.
For spinning BBHs, there are two families of EOB models, one by A. Buonanno's group and collaborators, e.g. [93,506,505,410,409,45,47,33,46,407,133], and one by T. Damour, A. Nagar, and their collaborators, e.g. [158,157,160,153,162,148,161,376,393,394,142,149,387,40,391,393,157,69,388,387]. These BBH EOB models both reproduce the PN and test-mass results in the appropriate limits, use the same energy mapping to the effective phase spaces, and the effective description of a test particle in an effective spacetime. They differ in choices made for the remaining freedom, for instance, the non-analytic functions assumed for the metric potentials and similar aspects that are not fixed by any analytical knowledge, as well as how effects of spins are included, i.e., imposing the exact limit of a spinning particle in Kerr spacetime or not. Both of these models also incorporate information from NR  Fig. 10 Tidal effects in the EOB potential A(r) which determines the energetics of circular orbits. The grey curve illustrates the test-particle limit of a Schwarzschild BH. Comparing this to the black curves shows the change in the potential due to finite mass ratio, where two different choices of non-analytic functions for resumming the available PN information are shown as the solid and dotted curves. These point-mass results change to those illustrated by the red and blue curves when adiabatic tidal effects are included, where the quadrupole (blue) gives the main contribution and the octupole (red) a small additional imprint. The curves terminate at the light ring of the tidal EOB model. The bottom panels shows the difference in the potential with respect to the Schwarzschild BH.
by calibrating parameterized, unknown high-order PN coefficients to NR data to obtain an accurate description of GW signals.
Tidal effects can be included in EOB models in different ways [22,395,66,67,36,155,268,497,75,77]. For simplicity, we will give details for the nonspinning BBH baseline models and note that the inclusion of tidal effects carries over to the spinning baseline models. The effect of the spin-induced multipole moments is also included in these EOB models [390,314]. As noted above, several other matter effects from spins together with tidal effects have not been calculated in detail. The qualitative impact of tidal effects on the EOB dynamics is illustrated in Fig. 10, which shows the gravitational potential A(r).

Brief summary of EOB models for nonspinning binaries
The conservative dynamics within the EOB approach are described by the Hamiltonian where the effective Hamiltonian H eff is that of a particle of mass µ. This energy mapping also emerges in quantum-electrodynamics [98], scattering theory [144,532,534,145], and from explicit comparisons of the action variables [100]. For nonspinning binaries the the effective Hamiltonian is given by Here, p φ is the azimuthal angular momentum and p r * = p r / √ D is the rescaled radial momentum; all the momenta are per unit reduced mass. The EOB potentials A and D, describing the effective spacetime, can be written as A = A pm + A tidal and D = D pm . The different choices for the point-mass potentials A pp are either Eqs. (A1)-(A2h) of Ref. [497] with the calibration parameter determined most recently in [93] or the resummation from [22,395,66]. The potential D is either taken from Eq. (A4) of Ref. [497] or [395].
The EOB equations of motion are The factors of A/ √ D appear because p r * and r are not canonically conjugate variables. The gravitational radiation reaction forces F φ and F r are constructed from where the energy fluxĖ rad is given bẏ The sum is only over positive m since |h F −m | = |h F m |. The factorized EOB waveforms for point masses have the form [152,148,147] Each of these terms depends on the EOB canonical variables and the parameters of the binary, see e.g. [93,395] for the explicit up to date expressions.

Adiabatic tidal effects using PN information
Adiabatic tidal effects on the energetics of the binary system can be captured in the EOB approach through tidal contributions to the potentials of the form [77,155,535] Here A (A) are the Newtonian −th multipolar tidal potentials given by The 2PN corrections to the tidal potentialÂ (A) from Eqs. (6.6) and (6.18) of Ref. [77] are given bŷ Tidal effects also influence the dissipative sector since the moving multipole moments of the NS contribute to the gravitational radiation. In the EOB model this is accounted for by adding to the waveform modes of Eq. (63) a tidal contribution h tidal m to the GW modes so that The explicit results for the adiabatic tidal terms h tidal m are given in Eqs. (A14)-(A17) of Ref. [163] for ≤ 3.

Dynamical f-mode tidal effects
The effects of dynamic tides from the NS's f −modes within the EOB framework are discussed in detail in Refs. [497,268]. Here, we consider only the approximate model that is used in practical data analysis applications. In this model, the potential from Eqs. (64) and (66) is used but with k multiplied by a frequency-dependent enhancement factor such that k → k k dyn witĥ with Ω = 3/8, and Ω = M 1/2 r −3/2 . The quantity Q m is given by where the functions F S and F C are Fresnel sine and cosine integrals respectively using the conventions in Mathematica. The quantitiest and m are defined aŝ Similar to the treatment for the conservative dynamics, the effect of dynamic f −mode tides can be incorporated in the dissipative sector in an approximate way by multiplying the occurences of k in h tidal m in Eq. (67) by an effective function k → k kdiss dyn . For = 2 this function is given bŷ wherek (A) 2 dyn is the enhancement function for body A in the conservative dynamics from Eq. (68).

Adiabatic tidal effects with strong-field enhancement
Gravitational self-force calculations have recently computed tidal invariants that contain information about strong-field tidal effects in the limit of small mass ratios, to linear order in m A /M [185,75]. These results have been augmented in Refs. [75,66,22,395] by a term ∝ m 2 A /M 2 that would describe currently unknown second-order self-force effects. Specifically, in this model Eq. (64) is employed with the potentialÂ (A) given bŷ is obtained from Eqs. (7.24)-(7.27) of Ref. [75]. In the model of Refs. [66,22,395] the choices for the unknown parameters are a 2GSF 2 = 337r −2 /28 and p = 4. The radius of the light ring r LR is obtained from the conservative EOB dynamics by solving with the potentialÃ = A pp + A tidal adPN from the PN model in Eqs. (64) and (66). These tidal models are included in a slightly different EOB model than the dynamical f -mode effects. There are a few minor differences in the pointparticle sector such as the factor N m in the waveform modes of Eq. (63), a different resummation of the potentials A pm and D, and the arguments of factors in h F m which involve powers of v φ = (∂H EOB /∂p φ ) −2/3 Ω are evaluated for circular orbits instead of v = Ω 1/3 , Padé-resummation of higher-order hereditary terms in h F m , and a different calibration of the BBH model to NR.

Phenomenological models
Over the last years, there has been an active development of phenomenological BNS waveform models. The main idea behind the existing models is the possibility to augment BBH approximants with a tidal description incorporating finite size effects of the NSs. Thus, a fundamental assumption is that the total GW phase can be decomposed into where φ pm denotes the point-mass, φ SO the spin-orbit, φ SS the spin-spin, and φ T the tidal phase contribution. We will focus in the following on an accurate representation of the tidal phase φ T that is required for a proper modelling of BNS systems. However, in addition and as pointed out in Sec. 2, also the spin-spin interaction, φ SS , or higher order spin-contributions incorporate information about the EOS due to the presence of the quadrupole moment. Therefore, φ SS needs to be adjusted to describe the BNS coalescence in comparison to BBH systems. As in the time domain, the frequency domain GW phase is typically decomposed as In the following subsections, we will discuss different phenomenological approaches to describe φ T and ψ T .

NRTidal-approximants
We start our discussion with the newest NRTidal version (NRTidalv2) described in Ref. [181]. Earlier versions use similar principles and assumptions and therefore will not be discussed to avoid repetition, but we refer to [174,183,172] for more details.
The main idea of the NRTidal [174,183] approach is to provide a closedform approximation for the tidal phase contribution φ T and ψ T . In addition NRTidalv2 incorporates spin-spin and cubic-in-spin effects up to 3.5PN, these contributions depend on the quadrupole and octupole moment of the individual NSs. These effects are included following the PN framework [181]. Furthermore for NRTidalv2, an additional amplitude correction is added.
As discussed, tidal phase contributions enter at the effectively 5th PN order and analytic knowledge exists up to the 7.5th PN order (with incomplete knowledge at 7PN); cf. Sec. 4. Since the current NRTidal models use only leading and next-to-leading order mass-ratio effects, we restrict the parameters c A,B 1 , c A,B 3/2 , c A,B 2 , c A,B 5/2 in Eqs. (54c) to their equal-mass values. Thus, the superscripts A and B are discarded in the following.
As for the time domain, some of the Padé parameters are constrained to obtain the correct PN limit up to 7.5PN (setting unknown terms at 7PN to zero):ñ The remaining, fitting coefficients arẽ In addition to the phase correction the NRTidalv2 approach also includes tidal amplitude corrections. While the GW amplitude is in general less important than a correct phase model, one needs a realistic and accurate estimate of the GW amplitude to measure the source distance properly, thus, additional tidal amplitude corrections might play a role in the measurement of the Hubble constant [477,5] or in determining if a GW source is gravitationally lensed [411].
The amplitude correction is directly derived in the frequency domain and calibrated in such a way that the 6PN tidal corrections as well as the phenomenological corrections of [290] are recovered: with d = 13477.8.
To truncate the waveform after the moment of merger, a Planck taper [372] is added at the merger frequency with n 1 = 3.354×10 −2 , n 2 = 4.315×10 −5 , d 1 = 7.542×10 −2 , d 2 = 2.236×10 −4 . The parameter ω 0 = 0.3586 is chosen so that for the equal-mass cases the correct non-spinning BBH limit [262] is recovered. The tapering ends at 1.2 or 1.3 times the merger frequency, ω mrg , depending on the exact implementation and NRTidal model version. The final amplitude is given byÃ we note again thatÃ NRTidal T = 0 for the first/original NRTidal version. Because of the smooth frequency and amplitude evolution even after the moment of merger, the Planck taper only introduces negligible errors and does not lead to biases in the parameter estimation; cf. Ref. [189].

Phenomenological model of Kawaguchi et al.
The phenomenological waveform model extension of Kawaguchi et al. [290] follows a similar approach as the NRTidal framework, i.e., it uses the assumption that the GW phase can be split into individual components, Eq. (75). The main difference between [290] and the NRTidal model are: (1 + aΛ 2/3 x p ) 1 +c 1 x +c 3/2 x 3/2 +c 2 x 2 +c 5/2 x 5/2 (91) to include non-linear tidal effects instead of higher order PN effects as in Eq. (84); (iii) The model is valid up to a frequency of 1000Hz to ensure that no postmerger frequency contributions might effect the inspiral approximant.
In addition, the model includes a tidal amplitude correction of the form

Waveform-model Comparison
We compare different frequency-domain tidal models in Fig. 11. The plot shows the 6PN, 7PN, and 7.5PN tidal predictions. We find that the 7.5PN approximant is the least attractive and that all PN predictions are less attractive than the phenomenological GW models. The original NRTidal approximant reduces to the 6PN curve, while NRTidalv2 and the model of Kawaguchi reduce, by construction, to the 7.5PN approximation. It also becomes obvious that up to a frequency of 1000Hz or even beyond, the NRTidalv2 and Kawaguchi model predict almost identical tidal effects. This changes for very large tidal deformabilities, but for configurations similar to GW170817 as shown in Fig. 11 withΛ = 392.1, systematic waveform uncertainties seem to be negligible; see also [399]. Extrapolating the model of Kawaguchi et al. [290] beyond its validity region up to the merger leads to an overestimation of tidal effects. Considering the difference between NRTidal and NRTidalv2, we find that NRTidal predicts larger tidal effects than NRTidalv2, consequently, when the model is employed for the analysis of real GW signals, one can expect that it will predict smaller tidal deformabilities to compensate for the larger tidal phase contribution. Fig. 11 Comparison of the tidal phase of different PN and phenomenological waveform approximant as discussed in the main text; cf. [181] for further details.

Analyzing BNS signals
As emphasized in Sec. 2, the primary difference of a BNS coalescence from a BBH coalescence comes from the presence of finite-size effects. From the point of view of estimating source properties, this involves parameterizing the matter effects and performing data analysis on a larger parameter space.

Estimating source properties
The first analysis to establish the possibility of using a Bayesian approach to estimate BNS source properties was presented in [169] using simulations of BNS mergers. This initial work was extended in [17], where higher tidal PN orders and the dominant spin-induced quadrupole moment were included. These works found that constraints on the binary properties can be improved by combining multiple observations. Ref. [315] put forward a way to parametrize the EOS and combined the tidal deformability parameter from multiple events to constrain EOS information. All this work however focused on using the PN inspiral-only waveform model, TaylorF2. Hybrid waveforms were constructed using NR simulations and analytical inspiral waveforms and improvement of measurability of matter effects were investigated in [453]. The first parameter estimation simulations using different waveform models were performed towards a study of systematics connected to GW170817 in [12]. Since the observation of GW170817, more studies on parameter estimation of BNSs have been carried out, as discussed below.
Bayesian inference relies on using Bayes' theorem to estimate the posterior probability distribution functions. This is done by calculating the likelihood of the observed data and folding in prior assumptions about our system, in our case, the parameters characterizing the source [529,530]. The information from individual parameters characterizing the source is encoded in posterior probability distributions Here, θ represents the parameter set and H s is the hypothesis that a GW signal depending on the parameters θ is present in the data d. The parameter set of a BNS source consists of the usual parameters describing a BBH source {m A , m B , χ A , χ B , θ, φ, ι, ψ, D L , t c , ϕ c }, and in addition the tidal deformability components Λ A and Λ B . The likelihood of obtaining a signal h(t) in data stream d(t) = h(t) + n(t), is proportional to In addition to inferring parameters characterizing a BNS system, Bayesian analysis also allows us to perform hypothesis selection. This enables us to do model selection between different EOSs. Evidence for each hypothesis is computed by assuming an EOS and estimating source properties of a BNS source. By computing evidence for another EOS, the odds in favor of one EOS over another is calculated by taking the ratio of evidences. This is also called Bayes' factor between two EOSs. For any single EOS, the evidence is given by By computing p(d|H EOS2 , I) for another EOS, the ratio between their evidences Ref. [169] looked into combining multiple sources to distinguish between the EOSs they follow and concluded that ∼ 20 sources would be sufficient to tell apart a stiff, moderate, or soft EOS. With the additional refinements in their analysis and waveform models in [17], it was found that imposing a stricter prior on NS masses, Gaussian in this case, requires many more signals to make an accurate distinction between different EOSs. In the future era, the postmerger phase of a GW signal will have a higher SNR and therefore be visible in the detector's band, leading to a measurement of the frequencies of the post-merger spectrum; see e.g. [94,121,521,96]. Besides extracting the source parameters, future detections might also allow us to extract EOS-insensitive quasi-universal relations from GW data [472].

Analyses of real signals
Estimating the source properties from GW170817 provided the first constraints on NS-EOS from a GW observation. Within a Bayesian setting, conservative estimates on tidal deformability and hence the EOS, were placed using GW170817. This was done by using prior distributions independently on the two dimensionless component tidal deformabilities [7,12], i.e., even allowing the that the two stars have different EOSs. The analyses in [12] used different waveform models, combining different BBH baselines and different tidal prescriptions to arrive at comparable estimates. The waveform models used were TaylorF2, IMRPhenomD_NRTidal, IMRPhenomPv2_NRTidal, and SEOBNRv4_ROM_NRTidal. The models differed in their treatments of the pointmass baselines, the spin-effects as well as the tidal information. The consistency among results from these models confirmed that systematic biases from model inaccuracies were not dominant for GW170817 with the present detector sensitivities.
The analyses placed a uniform prior on each of Λ i between [0,5000] and estimatedΛ by effectively using a prior uniform onΛ. Two different ranges of uniform spin priors were used; following observations of Galactic BNS systems, an upper limit of 0.05 on the individual components' spin magnitudes were placed. In addition, to be non-conservative, upper limits of 0.89 were used. Estimates for intervals of all intrinsic parameters including component masses as well as the tidal deformabilityΛ were found to be stronger for the narrower spin priors. With wider spin priors, the component masses were estimated to lie between 1.00 M and 1.89 M . With narrower priors on spins, the component masses were found to lie between 1.16 M and 1.60 M . The estimate on the upper bound ofΛ was 630 for wider spin priors, whereas within 90% credible interval, an estimate of 300 +420 −230 was made. The posteriors onΛ from the high spin priors using the four different waveform models mentioned above is shown in Fig. 12. Constraints on NS radii and EOS were obtained in [10]. This analysis was restricted to spins observed in NSs in Galactic binaries and therefore to a narrower range. Two methods were used to constrain the NS radii and EOS. Firstly, EOS-insensitive relations were used to employ priors on tidal deformabilities to enforce that both NS components follow the same EOS. Secondly, a method to parameterize directly the EOS p(ρ) was used. Radii estimates of R A = 10.8 +2.0 −1.7 km and R B = 10.7 +2.1 −1.5 km at 90% credible intervals were obtained for the primary and secondary components by imposing the EOSinsensitive relations. Further constraining bounds of R A = 11.9 +1.4 −1.4 km and R B = 11.9 +1.4 −1.4 km were obtained from the method of parametrization of the EOS, in addition to imposing the maximum mass of a NS to be > 1.97M . Ref. [110] updated quasi-universal relations based on priors imposed from constraints on presure and density of the EOS from analyzing GW170817. A similar analysis was also carried out independently in [167] although there are differences in the waveform models as well as the universal relations used, the results have been overall consistent. Preliminary conclusions were drawn in [7,  12] about EOSs preferred by the GW data. A detailed study of model selection was carried out in [14]. The data in this study was analysed with two ranges of priors, the narrow prior included uniform priors on spin magnitudes up to 0.05, as well as a Gaussian prior on component masses. The wider prior included uniform spin magnitudes up to 0.7, and a uniform prior on component masses. Only the stiffest EOSs could be ruled out completely from GW data. Using only the GW data without information from it's EM counterpart, we cannot rule out the BBH hypothesis either.
Ref. [289] pointed out caveats on quoting lower bounds on quantities like the tidal deformability from use of uniform prior in Bayesian inference. However, we note that the upper bound onΛ is quoted to be constrained from GW data analysis. In [13], the second observation of a BNS, GW190425, in the GW spectrum and the heaviest chirp-mass of a BNS was reported. Being such a heavy mass event, the constraints on tidal information did not improve on the existing ones from GW170817. However the constraints extracted were found to be consistent with those obtained from GW170817. Due to the high masses, the possibility that this event was gravitationally lensed [411] or that one or two components have been BHs can not be ruled out, e.g., [252,309]. Incorporating knowledge about the EOS and the fact that no EM counterpart has been observed can be used to further constrain the properties of GW190425 [134, 220,385] 7.3 Uncertainties from waveform modeling Although tidal estimates between different waveform models were found to be in agreement for GW170817, with improved detector sensitivities, the systematic uncertainties will start becoming dominant [10]. Ref. [470] simulated non-spinning GW170817-like BNS sources in the era of Advanced LIGO and Advanced Virgo's design sensitivity to look for biases in the estimates of the tidal deformability parameter. This analysis was performed by using different BBH baselines combined with two different tidal prescriptions, the one derived from PN, and the NRTidal framework. The SNRs of these simulations were ∼ 90 and the analysis focused on BBH baselines derived from PN, EOB, and phenomenological models. The simulated sources were an equal-mass and an unequal-mass binary, both following the APR4 EOS. To estimate the effect of changing point-mass BBH baselines, the NRTidal framework was used by varying all four underlying models (IMRPhenomPv2, IMRPhenomD, TaylorF2, and SEOBNRv4_ROM). For the equal-mass source, the injected tidal deformability was recovered within 90% credible interval whereas in case of the unequalmass source, only the two phenomenological models recovered the injected value within a 90% credible interval; note that the phenomenological model IMRPhenomPv2_NRTidal was used for injection. Other non-tidal intrinsic parameters like total mass, mass ratio, and the effective spins (in this case, 0) were found to be recovered reliably and consistently across all models and for both the sources. In addition, systematics arising from using the two different tidal prescriptions by keeping the point-mass baseline the same, were also investigated. It was found that increasing the PN order did not make the estimates more accurate, instead, the half-integer PN orders result in making the binary's interaction potential repulsive and therefore lead to an overestimate of the tidal deformabilities; cf. Fig. 11. Even for this case, one finds that nontidal intrinsic parameter estimates of total mass, mass ratio, and the effective spin are reliably and consistently recovered across all tidal descriptions. Out of the two sources, it was shown that unequal-mass sources showed greater biases in tidal deformability estimates. The biases in tidal deformability for both the equal-mass and unequal-mass source for the Phenom, EOB, and TaylorF2 baselines when the tidal information is kept the same, NRTidal, in all cases, are shown in Fig. 13. Fig. 14 shows the bias in tidal deformability when the baseline is the same, IMRPhenomPv2, and the tidal descriptions are changed.
As an extension to this work, Ref. [471] investigated the biases on tidal parameter estimates from highly spinning BNSs, also in the era of design sensitivity, and the importance of spin-induced quadrupole moments. Particularly, it was found that neglecting EOS-dependent spin-induced quadrupole moments leads to a bias in parameter estimation once the dimensionless spin magnitude of a component NS is χ > 0.2. In addition, [471] studied different spin magnitudes and orientations as well as the advantage of having envisaged EM counterparts like a short GRB or a kilonova or both were studied. For aligned-spin configurations, in-plane spins, and spin vectors oriented at 45 • to the orbital angular momentum vector, the waveform models including the EOS-dependent spin-induced quadrupole moment performed better in estimating the tidal deformability parameter as spin magnitudes increased. In addition, having EM counterparts produced improved slightly the estimates of tidal deformability, however, no noticeable improvement in estimates of non-tidal parameters were observed.
While the above works focused on tidal estimates from the inspiral only, Ref. [189] used hybrid simulations to include the inspiral as well as post-merger regimes in a BNS coalescence. It was shown that non-tidal parameters were estimated reliably for SNRs as low as 25 and for soft EOSs, even when analyzing the signal with a point-particle waveform model without tidal effects. For stiffer EOSs or larger SNRs, significant biases in estimates of non-tidal parameters were observed when analyzing a BNS signal with a point-particle waveform model. Additionally, it was found that neglecting the post-merger effects in the recovery of a BNS signal would not lead to significant biases in parameter estimation in the second generation era of detectors. In [375], the authors calculated a higher order PN term at 5.5PN and included it for parameter estimation of tidal waveforms on a non-spinning injection. The estimates using this waveform model which includes the higher order PN term is more consistent with the signal and the estimates obtained from the full IMRPhenomD_NRTidal model than lower order PN approximants.

A GW data-analysis outlook
The envisaged era of third generation (3G) detectors will lead to increase in sensitivities at frequencies as low as f ∼ 5 Hz. This will result in GW signals spending several cycles in the detector's band and much higher SNRs, BNS signals may spend up to hours in the bandwidth of the detector. As the signal  enters the band much earlier, this will give the opportunity of early-warning alerts for possible EM counterparts [21]. Within the large timescales of the signal lasting in band, the detectors' power spectrum may not be stationary and the Earth may also involve further subtleties in GW data analysis with such long signals. Ref. [456] investigated overlapping signals in the third generation era, the work however focused on the problem of GW detections only, finding that overlap between two or more signals did not affect the performance of search algorithms unless the lower frequency cutoff was reduced. For long signals, the possibility of BNS signals being overlapped with other GW signals from BNSs or other GW sources is much higher. As no parameter estimation studies have been done on such overlapping signals, it remains to be seen whether retrieving source properties will be affected. The above work also identified critical issues for matched filtering for a signal lasting in band for up to a few days, making our current availability of computational resources inadequate. While techniques like faster sampling algorithms and waveform models will improve the ability to perform data analysis, we need better approaches to be able to extract information from the signals in the 3rd generation detector era.

Summary
In the last years, GW astronomy has been revolutionized by the direct detection of GWs from BH and NS binary systems. In particular the possibility to detect GW and EM signals in combination is a driving force in the field of multi-messenger astronomy and highlights the importance of the detection of BNS configurations. In this review, we have focused on the measurement of tidal information from GW signals in the current GW detector era, i.e., our discussion focused primarily on the inspiral and on matter effects that are incorporated in state-of-the-art models. We reviewed key aspects in the field of NR (Sec. 3), PN theory (Sec. 4), the EOB formalism (Sec. 5), and phenomenological waveform approximants (6) to model GW signals for BNSs. We finalized our discussion by a short review about the detected GW signals GW170817 and GW190425, and by briefly describing GW data analysis results focusing on the imprint and extraction of tidal effects from GW signals.
In conclusion, the planned improvements in experimental capabilities in the next years to decades will open unprecedented opportunities for fundamental physics and astrophysics with BNSs. However, capitalizing on the science potential will require significant further effort on addressing the challenges in numerical and analytical modeling as well as data analysis outlined above.