1 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, 8, 286] arising from a BNS coalescence detected in August 2017 marked a breakthrough in the field of multi-messenger astronomy. This multi-messenger observation represented a major leap forward in a number of research areas. It enabled a new and independent measurement of the Hubble constant [5, 11, 137, 172, 178, 278]; it proved NS mergers to be a major cosmic source of heavy r-process elements, e.g., [139, 285, 286, 501]; it revealed that BNS mergers can act as the central engine for gamma-ray bursts (GRB) [8]; it enabled an accurate measurement of the propagation speed of GWs [8]; it ruled out a number of alternative theories of gravity [39, 140, 201]; and it allowed to place limits on the equation of state (EOS) of cold matter at supranuclear densities, e.g.,  [7, 10, 12, 30, 56, 107, 136, 138, 142, 168, 178, 358, 366, 391, 452, 454, 469, 474, 493].

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. [538]. 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. [35, 38, 181, 202, 222, 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. [15, 41, 164, 283, 321, 419], the performance of Post-Newtonian (PN) waveform approximants becomes increasingly inaccurate towards the merger, e.g. [68, 184, 190, 209, 276, 477, 544], which shows the need for possible alternatives. Therefore, most state-of-the-art time-domain tidal waveform models [22, 36, 37, 67, 268, 274, 395, 397, 400, 505] are based on the EOB description of the general relativistic two-body problem [101, 156]. 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 [22, 179, 274]. 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 [397]. 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. References  [313, 423, 517] develop phenomenological black hole-neutron star (BHNS) approximants based on NR data. Reference [48] transforms NR simulations of binary black hole (BBH) systems by adding PN tidal effects, and Refs. [325, 326] 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. [175], in which existing BBH models have been augmented by an analytical closed-form expression correcting the GW phase to include tidal effects  [175, 182, 184, 290].

In this review article, we discuss the BNS and merger dynamics in Sect. 2. In Sect. 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 Sect. 4, we give a brief overview about the existing PN knowledge for the description of the BNS coalescence. In Sect. 5 and Sect. 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 this work geometric units are used by setting \(G=c=1\), and, additionally, \(M_\odot =1\). In some places units are given explicitly to allow a better interpretation. Further notations are \(M=M_A + M_B\) for the total mass of the system, \(\chi _A,\chi _B,\varLambda _A,\varLambda _B\) for the dimensionless spins and tidal deformabilities of the individual stars. The mass ratio of the system is \(q=M_A/M_B\) and the symmetric mass ratio is \(\nu =M_A M_B/(M_A+M_B)^2\). We define the labeling of the individual stars so their masses satisfy \(M_A\ge M_B\).

Fig. 1
figure 1

NR simulation of a BNS merger showing the GW signal and the matter evolution. Top panel: GW signal emitted during the last orbits before the merger (late-inspiral phase) and during the postmerger phase of the BNS coalescence. Bottom panel: Rest-mass density evolution for the inspiral (first panel), the merger (second panel) and the postmerger phase after the formation of the black hole (third panel)

2 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., [195, 331, 401, 403, 409, 480], as well as a kilonova, see [382, 511] 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 Fig. 1, where we show approximately \(60~{\mathrm{ms}}\) around the moment of merger. The figure includes the density profile obtained from a NR simulation during the late inspiral (bottom left panel), around the moment of merger (bottom middle panel), and the final BH surrounded by an accretion disk (bottom right panel).

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.

2.1 Inspiral

The inspiral of a BNS coalescence reaches frequencies up to about 1–2 kHz, where the exact frequency depends on the binary properties. Current GW detectors are most sensitive at frequencies around \(\sim 100\, {\mathrm{Hz}}\) 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–30 Hz or at even lower frequencies for 3rd generation GW detectors.

Fig. 2
figure 2

Accumulation of information about mass, spin, and tidal parameters per logarithmic frequency interval. Left panel: aLIGO design sensitivity (zero detuned high power configuration) [3, 4], right panel: Einstein Telescope [1, 265]. 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) [383, 481] combined with the leading-order tidal [215] and spin effects [89, 147, 292]. Adapted from Ref. [257]

2.2 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 the presence of a companion are small, with a description of the interaction zone where the NSs behave almost as point masses with small corrections due to their finite size [214, 449], see also  [185, 281, 302, 303, 521, 543]. For weakly self-gravitating bodies described by PN gravity see also the seminal series of papers by Damour, Soffel, Xu [165]. As will be discussed in detail in Sect. 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 [87]. 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 [220, 335, 442, 473] for comprehensive review and references.

Fig. 3
figure 3

Cartoon depicting the definition of tidal deformability. The tidal field \({{{\mathcal {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 \({{{\mathcal {Q}}}}\). The multipoles are also impacted by spin effects, and dynamical tidal effects

2.2.1 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 time-varying 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 [299, 316, 376, 491], with the effect of the quadrupole \(\ell =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 [215]. The tidal deformability \(\lambda _\ell \) for each \(\ell \)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_\ell \) to the strength of the multipolar tidal field \({{{\mathcal {E}}}}_\ell \) associated with the curvature due to the companion as will be discussed in Sect. 4:

$$\begin{aligned} Q_\ell =-\lambda _\ell \, {{{\mathcal {E}}}}_\ell \end{aligned}$$
(1)

The associated effect on GWs steeply increases with smaller separation or higher GW frequency, however, the information accumulates over all frequencies (Fig. 3).

2.2.2 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 solutions of General Relativity, the multipolar structure of the exterior spacetime is fixed by the no-hair property in terms of only their mass and spin [114, 248, 261]:

$$\begin{aligned} M^{{\mathrm{BH}}}_\ell +i S^{{\mathrm{BH}}}_\ell =m(i\chi m)^\ell , \end{aligned}$$
(2)

where \(M^{{\mathrm{BH}}}_\ell \) and \(S^{{\mathrm{BH}}}_\ell \) are the \(\ell \)th mass and current multipole moments of the black hole, and m and \(\chi =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

$$\begin{aligned} M_{\ell =2}^{{\mathrm{NS}}}\approx -\kappa \chi ^2 m^{3}, \end{aligned}$$
(3)

The key characteristic matter-dependent coefficients \(\kappa \) are the rotational Love numbers [73, 258, 311, 386, 437]. Studies in Refs. [428, 551] 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 than fundamental-mode tidal effects.

2.3 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.,  [23, 52, 131, 169, 387, 450, 494, 499, 547]. 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. [388] for recent progress improving the accuracy of postmerger simulations.

Independent of these caveats, the qualitative picture, as depicted in Fig. 4, is robust and the possible merger outcomes are generally classified in four different categories [51]: prompt collapse, hypermassive NS (HMNS), supramassive NS (SMNS), and massive NS (MNS).

Fig. 4
figure 4

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

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, [56, 136, 178, 367, 454], 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_{{\mathrm{thr}}}\) as computed in [18, 53, 304].

If the mass is smaller than the threshold mass and the baryonic mass is below the maximum 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 [51], but has been applied to make qualitative predictions, e.g. [272, 564]. A proper classification incorporating the temperature dependent effects are currently missing, but see [106, 284, 502] 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 [451].

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_{{\mathrm{peak}}}\)  [57, 64, 132, 133, 451, 470, 509, 529]. In addition also a number of other, sub-dominant frequencies are present, e.g., Refs. [57, 133, 470]. 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 [451]. Recently, Refs.  [97, 122, 529] 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–5 is required. This corresponds to a total SNR of 100–200, which would allow measurements about tidal effects from the inspiral GW signal.

3 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.,  [65, 66, 71, 125, 174, 177, 180, 183, 191,192,193,194, 226, 242, 275, 287, 288, 295, 333, 389, 425, 453, 475, 487, 507, 530, 532]; (ii) the development of numerical schemes allowing many-orbits simulations with high-accuracy GWs, e.g., [63, 72, 175, 274, 277, 294, 455]; (iii) the increasing realism of microphysical schemes that incorporate magnetic effects, finite-temperature EOSs, composition effects, and neutrino transport, e.g., [26, 130, 224, 228, 229, 231, 235, 240, 296, 297, 388, 405, 410, 411, 450, 468, 474, 485, 486, 493,494,495,496]; (iv) the study of mass ejecta and EM counterparts, e.g.,  [54, 96, 177, 183, 229, 232, 271, 273, 333, 432, 453, 545].

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., [35, 38, 468, 472], for discussions about effects of magnetic fields or neutrino radiation.

3.1 Imprint of the binary properties and waveform catalogs

A number of different groups have performed a large set of BNS simulations over the last years. The most ambitious simulation campaigns have led to the construction of the CoRe databaseFootnote 1 [181] and the SACRA Gravitational Waveform Data BankFootnote 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. [54, 272, 510]. 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.

3.1.1 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. 4.

Based on the constraints derived from GW170817 one expects that the mass of an individual NS lies within \(M \in [1.0,2.3 M_\odot ]\)  [327, 366, 408, 469, 474, 493, 498].

As discussed in Sect. 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. [252].

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. 5. 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. [126].

Fig. 5
figure 5

Examples of the link between the microphysics of cold, dense matter as encoded in the equation of state (EOS) and global NS properties such as tidal deformability as a function of mass. The connection occurs through the Einstein Field equations and energy–momentum conservation for linearized perturbations to an equilibrium NS configuration. The quantity \(\varLambda =G\lambda _{\ell =2}(c^2/Gm)^5\) is a dimensionless tidal deformability coefficient characterizing the ratio of the induced quadrupole to the strength of the tidal field in the adiabatic limit. The colored curves correspond to a few illustrative EoS examples. The parameter that varies along the tidal deformability curves is the central density. For each EoS, there is a maximum central density and mass for which the material can resist gravitational collapse to a black hole. Adapted from Refs. [249, 267]

Fig. 6
figure 6

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 [181]

3.1.2 Mass ratio

Based on population synthesis models a maximum mass ratio [\( q = \frac{M_A}{M_B} \ge 1\)] of \(q_{\max } \simeq 1.8 {-} 1.9\) [180, 187] might be possible for a BNS system. While observation of isolated NSs distribution would theoretically allow maximum mass ratios of \(q_{\max } \simeq 2\), the largest observed mass ratio in BNSs is \(q\sim 1.3\) [330, 369]. Over the last years, the NR community managed to simulate systems with large q, at the boundary of physically realistic scenarios, see [180, 525] for simulations with \(q > rsim 2\).

With increasing q, one finds that the emitted GW energy decreases during both the inspiral and the postmerger phase; cf. Fig. 14 of [183] and that the merger frequency decreases roughly as a function of \(\propto q^{-1/2}\), cf. [184].Footnote 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 [183, 295, 295]. 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 [529].

3.1.3 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 \(\chi \sim 0.7\), which corresponds to spin periods of less than 1 ms [352].

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 [359] is the fastest spinning NS in a binary system with a rotation frequency of 239 Hz and PSR J1946+2052 [506] 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 \(\chi \sim 0.05\) at merger.

In the past, consistently evolved spinning BNS simulations have been presented in Refs. [65, 173, 177, 180, 193, 390, 533] and precessing NR simulations are presented only in [125, 174, 180, 507].

With respect to the dynamical evolution and energetics of spinning NS binaries during the inspiral one finds that even for astrophysical realistic spin magnitudes (\(\chi \lesssim 0.05\)) spin–orbit interactions are larger than tidal effects up to or shortly before the contact of the two stars [65]. 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. [17, 209, 478].

Studying the energetics during the postmerger, one finds that spin aligned systems can create merger remnants with larger angular momentum than non-spinning systems. Consequently, collapse to a BH happens at later times [125]. We note that simulations also indicate that characteristic frequencies in the postmerger, which generally also follow quasi-universal relations  [57, 58, 64, 133, 470, 509, 510], might shift under the presence of spin.

3.1.4 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. Reference [527] gives a conservative estimate of a BNS merger rate of \(0.5\, {\mathrm{yr}}^{-1} \,{\mathrm{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 strong-field 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. [124, 127, 242, 299, 384, 429, 533, 558]. In addition, it also allows to understand possible resonance effects when the orbital frequency approaches the f-mode frequency and resonant effects influencing the evolution of the system [28, 483, 505].

3.2 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 to develop and test existing GW models; see e.g. [63, 175, 182, 274, 294, 295, 388, 456, 531] for further discussion. The main difficulty for the production of high-quality waveforms are:

  1. (i)

    the reduction of artificial eccentricity;

  2. (ii)

    the usage of high enough resolutions to reduce discretization errors;

  3. (iii)

    the development of numerical schemes showing clear convergence across multiple resolutions and configurations.

3.2.1 Eccentricity reduction

As an exemplary case addressing all listed items, we present a high-resolution simulation of [182]. 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 [180, 461]. The eccentricity reduced initial configurations are obtained with the SGRID code [180, 522,523,524,525]. 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 [227, 507] 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 [180, 525].

Figure 7 shows the eccentricity reduction procedure for this configuration. In most cases it is sufficient to perform 2–4 iteration steps to achieve final eccentricities \(\lesssim 10^{-3}\).

3.2.2 Errors through numerical discretization

The configurations which we want to discuss in the following have been simulated with the BAM code [63, 100, 176, 516] for six different resolutions covering the individual NSs with 64, 96, 128, 192, 256, and 320 points respectively.

Fig. 7
figure 7

Eccentricity reduction procedure. Top panel: NS track for one star for the different iteration steps during the eccentricity reduction. The configurations need to be evolved for a few orbits to estimate the parameters \((e,{\dot{r}})\) that are required for the next iteration step. Bottom panel: Proper distance of the system as a function of time. The artificial eccentricity is clearly visible as oscillations in the proper distance

The highest resolution (320 points in the finest refinement level) has a spatial resolution of \(0.047 M_\odot \approx 70 \ {\mathrm{m}}\) in the finest refinement boxes. The computational cost for all six resolutions is about \(\sim 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.

Fig. 8
figure 8

Top panel: Real part of the GW signal for the six different resolutions employing 64, 96, 128, 192, 256, and 320 points in the refinement levels covering the individual NSs as a function of the retarded time u. Middle panel: Phase difference between different resolutions. Bottom panel: Phase difference between different Richardson extrapolated waveforms or between a Richardson extrapolated waveforms and the waveform from an individual resolution. The vertical lines in each panel refer to the time of merger, i.e., the peak time of the GW amplitude for the individual resolutions. The dashed lines in the bottom two panels show the phase difference scaled to the next lowest pair of resolutions assuming second order convergence. Adapted from [182]

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\ {\mathrm{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\ {\mathrm{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.

3.2.3 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 [85] and references therein.

The best solution to obtain properly extracted results is the Cauchy Characteristic Extraction (CCE) [84, 467], see also [34, 255].

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\ldots N\) and then the phase and amplitude are extrapolated using a polynomial of the order K with \(K<N\):

$$\begin{aligned} f(u; r_j) = f_0(u) + \sum ^K_{k=1} f_k(u) r^{-k}_j \ . \end{aligned}$$
(4)

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\ {\mathrm{rad}}\) for our example shown in Fig. 8.

Another possibility is it to consider the next-to-leading-order (NLO) falloff behavior of the \(\varPsi _4\) multipoles [63, 356],

$$\begin{aligned} r\,\psi _{\ell m} = r\,{\ddot{h}}_{\ell m} + \frac{\ell (\ell +1)}{2r}\; r\,{\dot{h}}_{\ell m} + (r^{-2}) \ , \end{aligned}$$
(5)

and obtains an extrapolation formula for \(r\,\ddot{h}_{\ell m}\) which uses \(\psi _{\ell m}\) extracted at a given radius.

It is interesting to note that, as shown explicitly in Fig. 12 of [63], 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.

4 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 [87] for a comprehensive list of references on this decades-long effort; for other review see [233] and [482]. 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 [61, 62, 364], a Hamiltonian approach [151, 152], and effective field theory methods [92, 217, 218, 335, 342]. 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. [79], 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, 88, 519], 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, 91, 199] and effective field theory methods [334] based on tools developed in [82, 83, 219, 236,237,238, 243, 244, 335, 361, 445, 472].

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 \(\lambda \) 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 [87]. 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 \(\epsilon \sim G M/rc^2\sim v^2/c^2 \ll O(1)\) and is dominated by point-mass contributions. Small corrections from finite-size effects are included through a tidal expansion in the parameter \(\alpha =R/r\). This yields a double expansion where \(\epsilon \) and \(\alpha \) are treated as distinct parameters that correspond to different physical effects, although they scale in the same way with the orbital separation [214]. 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 [215, 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(\alpha ^5)\). Relativistic corrections have been computed up to \(O(\alpha ^5\epsilon ^2)\) for the quadrupole and \(O(\alpha ^7\epsilon ^2)\) for the octupole [78, 263]. 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.

4.1 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.

4.1.1 NS’s mass multipole moments in Newtonian gravity

In Newtonian gravity the gravitational potential U generated by a mass density \(\rho \) at a field point \({\varvec{x}}\) is a solution to Poisson’s equation \(\nabla ^2 U=-4\pi \rho \). The Green function solution is:

$$\begin{aligned} U(t,{\varvec{x}})=\int d^3 x^\prime \frac{\rho (t, {{\varvec{x}}}^\prime )}{|{{\varvec{x}}}-{{\varvec{x}}}^\prime |}. \end{aligned}$$
(6)

For points \({\varvec{x}}>{\varvec{x}}^\prime \) outside the mass distribution we can perform a Taylor series expansion around a moving reference point \({\varvec{z}}(t)\) to write Eq. (6) as

$$\begin{aligned} U=\int d^3 x^\prime \rho (t, {{\varvec{x}}}^\prime ) \sum _{\ell =0}^\infty \frac{(-1)^\ell }{\ell !} (x^\prime -z)^L\partial _L\frac{1}{|{{\varvec{x}}}-{{\varvec{z}}}|}. \end{aligned}$$
(7)

Here, \(L = a_1 a_2 \cdots a_\ell \) denotes a string of \(\ell \) indices, \(x^L=x^{a_1} \cdots x^{a_\ell }\), and the summation over repeated indices is implied. Likewise, the notation for \(\ell \) derivatives is \(\partial _L=\partial _{ x^{a_1}} \ldots \partial _{ x^{ a_\ell }}\) with \(\partial _x=\partial /\partial x\).

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

$$\begin{aligned} M_{{\mathrm{Newt}}}=\int d^3x\;\rho (t, {\varvec{x}}),\quad Q^{L}_{{\mathrm{Newt}}}=\int d^3x\;\rho (t, {\varvec{x}}) (x-z)^{<L>}, \end{aligned}$$
(8)

where \(M_{{\mathrm{Newt}}}\) is the mass monopole and \(Q^L\) the \(\ell \)th multipole moment. The angular brackets around indices denote the symmetric and trace-free projection of a tensor. For example, \(x^{<ij>}=x^ix^j-\delta ^{ij}|{\varvec{x}}|^2\), see Refs. [259, 518] for more details. We have defined the moments as their trace-free parts because the contribution to Eq. (7) involves the derivative \(\partial _L |{{\varvec{x}}}-{{\varvec{x}}}^\prime |^{-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 (\(\ell =1\)) is absent from the expansion (8). With these definitions the exterior potential takes the form

$$\begin{aligned} U(t,{\varvec{x}})=\frac{m}{|{{\varvec{x}}}-{{\varvec{z}}}|}+\sum _{\ell =2}^\infty \frac{(-1)^\ell }{\ell !}Q^L\partial _L\frac{1}{|{{\varvec{x}}}-{{\varvec{z}}}|}. \end{aligned}$$
(9)

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/|{{\varvec{x}}}-{{\varvec{z}}}|=( \sin \theta \cos \phi ,\sin \theta \sin \phi ,\cos \theta )\) and express (9) as an equivalent spherical-harmonic expansion

$$\begin{aligned} U(t,{\varvec{x}})=\frac{m}{|{{\varvec{x}}}-{{\varvec{z}}}|}+\sum _{\ell =2}^{\infty }\sum _{m=-\ell }^\ell Q_{\ell m}\frac{Y_{\ell m}(\theta , \phi )}{|{{\varvec{x}}}-{{\varvec{z}}}|^{\ell +1}}, \end{aligned}$$
(10)

where the relation between spherical and Cartesian multipole moments is [518]:

$$\begin{aligned} Q_{\ell m}=\frac{4\pi }{2\ell +1}{{{\mathcal {Y}}}}_{L}^{*\, \ell m}Q_L. \end{aligned}$$
(11)

The tensors \({{{\mathcal {Y}}}}_{L}^{\ell m}\) are defined by \(Y_{\ell m}={{{\mathcal {Y}}}}_{L}^{\ell m}n^{<L>}\) and involve complex coefficients that convert between the unit vector \({\varvec{n}}\) and spherical harmonic representations [518].

4.1.2 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^\mu (\sigma )\) where \(\sigma \) is an evolution parameter along the worldline. This is known as the worldline-skeleton approach [185], 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 [518, 519], the time–time component of the metric involves a potential \(U_{{\mathrm{eff}}}\) that is analogous to the Newtonian gravitational potential U:

$$\begin{aligned} g_{tt}=-\left( 1-2 U_{{\mathrm{eff}}}\right) , \end{aligned}$$
(12)

For a spherical object \(U_{{\mathrm{eff}}}=M/r\) where r is the distance from the body, while for a nonspherical, nonspinning body it has the asymptotic form

$$\begin{aligned} \lim _{r\rightarrow \infty } U_{{\mathrm{eff}}}=\frac{M}{r}+ \frac{3n^{<ij>}Q_{ij}}{2r^{3}}+{{{\mathcal {O}}}}(r^{-4}). \end{aligned}$$
(13)

In this setting, the \(\ell \)th mass multipole moment is associated with the coefficient of the term that falls off as \(r^{-(\ell +1)}\). This definition is equivalent to the Geroch–Hansen multipole moments for stationary spacetimes [250]. We note that the multipole moments can be defined by considering the energetics of test particle orbits [476], and that in a binary system a small ambiguity remains in the definition of the multipole moments [247, 521].

4.1.3 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^{{\mathrm{{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:

$$\begin{aligned} U^{{{\mathrm{ext}}}}(t,{\varvec{x}})= & {} U^{{{\mathrm{ext}}}}(t,{\varvec{z}})+(x-z)^j \left[ \partial _jU^{{{\mathrm{ext}}}}(t,{\varvec{x}})\right] _{{\varvec{x}}={\varvec{z}}}\nonumber \\&-\sum _{l=2}^\infty \frac{1}{\ell !}(x-z)^{L}{{{\mathcal {E}}}}^L_{{\mathrm{Newt}}}. \end{aligned}$$
(14a)

The coefficients \({{{\mathcal {E}}}}^L\) in (14) are the NS’s tidal moments [521]. We are interested in specializing the general expansion (14) to the case that \(U_{{\mathrm{ext}}}\) is sourced by a binary companion whose gravitational potential we denote by \(U_{{\mathrm{c}}}\). The NS’s tidal moments are then given by

$$\begin{aligned} {{{\mathcal {E}}}}^L_{{\mathrm{Newt}}}=- {\partial _{L}}U_{{\mathrm{c}}}(t,{\varvec{x}})\mid _{{\varvec{x}}={\varvec{z}}(t)}. \end{aligned}$$
(14b)

4.1.4 Tidal moments in GR

The relativistic tidal moments are projections of the curvature tensor due to the companion. The analogue to the Newtonian \({{{\mathcal {E}}}}_L\) is the ‘electric part’ of the companion’s Weyl tensor \(C_{\mu \alpha \nu \beta } \) as [521]

$$\begin{aligned} {{{\mathcal {E}}}}_{\mu \nu }=z^{-2}C_{\mu \alpha \nu \beta } u^\alpha u^\beta \end{aligned}$$
(15)

where \(u^\mu =dz^\mu /d\sigma \) is the tangent to the NS’s worldline and we have defined

$$\begin{aligned} z=\sqrt{-u_\alpha u^\alpha } \end{aligned}$$
(16)

The tidal tensor \({{{\mathcal {E}}}}_{\mu \nu }\) is symmetric and trace free. In the NS’s rest frame \(u^\mu =(-1,0,0,0)\) it is a purely spatial tensor \({{{\mathcal {E}}}}_{\mu \nu } u^\nu = 0\) and given by \({{{\mathcal {E}}}}_{ij}=C_{titj}\). Higher multipole tidal moments are defined by

$$\begin{aligned} {{{\mathcal {E}}}}_L=\frac{1}{(\ell -2)!} C_{\langle 0 a_10 a_2;a_3 \cdots a_\ell \rangle }, \end{aligned}$$
(17)

where the semicolon denotes covariant derivatives.

4.2 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.

4.2.1 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:

$$\begin{aligned} m\ddot{z}^j=m\frac{\partial }{\partial x^j }U^{{{\mathrm{ext}}}}(t,{\varvec{x}})\mid _{{\varvec{x}}={\varvec{z}}}-\sum _{\ell =2}^\infty \frac{1}{\ell !}Q^{L}{{{\mathcal {E}}}}_{jL}, \end{aligned}$$
(18)

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 \({{{\mathcal {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 [543]:

$$\begin{aligned} T=\frac{1}{2}\int d^3 x \rho \dot{{{\varvec{z}}}}^2+T^{{\mathrm{int}}}, \qquad V= \frac{1}{2}\int d^3 x \rho U_{{\mathrm{c}}}+V_{{\mathrm{int}}}. \end{aligned}$$
(19)

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

$$\begin{aligned} T=\frac{1}{2}\mu v^2+T_{{\mathrm{int}}}, \qquad V=-\frac{ \mu M}{r}+\sum _{\ell \ge 2}\frac{1}{\ell !}Q_{L}{{{\mathcal {E}}}}_{L}+V_{{\mathrm{int}}}. \end{aligned}$$
(20)

Here, M and \(\mu \) are the total and reduced mass respectively. The relative separation is \({\varvec{r}}={\varvec{z}}-{\varvec{z}}_{{\mathrm{c}}}\) with magnitude \(r=|{\varvec{r}}|\), and the relative velocity is \(v^2=\dot{{\varvec{r}}}\cdot \dot{{\varvec{r}}}\). The action for the binary dynamics is then given by [215, 316, 459]

$$\begin{aligned} S=S_{{\mathrm{pm}}}+\int dt \sum _{\ell \ge 2}\bigg [-\frac{1}{\ell !}Q_{L}{{{\mathcal {E}}}}_{L} +{{{\mathcal {L}}}}^{{\mathrm{int}}}\bigg ], \end{aligned}$$
(21)

where \(S_{{\mathrm{pm}}}=\int dt \, \left[ (\mu /2) v^2+\mu M/r\right] \) describes the orbital motion of point-masses and \({{{\mathcal {L}}}}^{{\mathrm{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 [75, 165, 211, 215, 256, 299, 300, 302, 317, 318, 387, 466, 493, 563, 564], 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 [117, 215, 300, 317, 460]

$$\begin{aligned} {{{\mathcal {L}}}}^{{\mathrm{int}}}=\frac{1}{ 2\ell ! \lambda _\ell \omega _{0\ell }^2} \left[ {\dot{Q}}_{L} {\dot{Q}}^{L} - \omega _{0\ell }^2 Q_{L} Q^{L} \right] . \end{aligned}$$
(22)

Here, the quantities \(\omega _{0\ell }\) denote the f-mode frequencies and we have neglected contributions from other quadrupolar \(\ell =2\) modes. The parameters \(\lambda _\ell \) are the tidal deformability coefficients. They are defined by considering the adiabatic limit, where the NS’s internal time scales \(\tau ^{{\mathrm{int}}}\sim \omega _{0\ell }^{-1}\sim \sqrt{R^3/m}\) are fast compared to the time scale of variations in the tidal field \(\tau _{{\mathrm{orb}}}\sim \sqrt{r^3/M}\), since

$$\begin{aligned} Q_{L}^{{\mathrm{adiab}}}=-\lambda _{\ell } \, {{{\mathcal {E}}}}_L. \end{aligned}$$
(23)

The tidal parameters \(\lambda _\ell \) are related to the NS’s dimensionless tidal Love numbers \(k_\ell \) [357] by \(\lambda _\ell =\frac{2}{(2\ell -1)!!} k_\ell R^{2\ell +1} \) and depend on the NS’s interior structure. For data analysis applications, it is useful to define the dimensionless tidal deformabilities

$$\begin{aligned} \varLambda _\ell =\frac{\lambda _\ell }{m^{2\ell +1}}. \end{aligned}$$
(24)

For adiabatically induced multipoles, \(dQ_L/dt=0\), the internal Lagrangian is only the elastic potential energy associated with the deformation \({{{\mathcal {L}}}}^{{\mathrm{int}}}_{{\mathrm{adiab}}}=- Q_{L} Q^{L}/(2\ell ! \lambda _\ell )\). In this limit the mass quadrupole finite size effects on the dynamics depend only on the orbital variables and \(\lambda _\ell \):

$$\begin{aligned} S_{{\mathrm{adiab}}}=S_{{\mathrm{pm}}}+\int dt \bigg [\frac{\lambda _\ell }{2\ell !}{{{\mathcal {E}}}}_{L}{{{\mathcal {E}}}}^{L} \bigg ], \end{aligned}$$
(25)

with \({{{\mathcal {E}}}}_L=-m_{{\mathrm{c}}} \partial _L r^{-1}\).

4.2.2 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 re-parametrizations, and replacing all time derivatives by covariant derivatives along the center-of-mass worldline [505]:

$$\begin{aligned} S = S_{{\mathrm{pm}}}+\int d\sigma \, \left\{ - \frac{z}{2} {{{\mathcal {E}}}}_{\mu \nu } Q^{\mu \nu }+ {z \over 4 \lambda z^2 \omega _{02}^2} \left[ \frac{D{Q}_{\mu \nu }}{d\sigma }\frac{D {Q}^{\mu \nu }}{d\sigma } - z^2\omega _{02}^2 Q_{\mu \nu } Q^{\mu \nu } \right] \right\} .\nonumber \\ \end{aligned}$$
(26)

Here, total derivatives are denoted by \(\frac{D}{d\sigma }=u^\beta \nabla _\beta \), where \(\nabla _\alpha \) 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 effective-field-theoretical approach [198, 245, 505] 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

$$\begin{aligned} S = S_{{\mathrm{pm}}}+\int d\sigma \left[ - \frac{z}{2} {{{\mathcal {E}}}}_{ij} Q^{ij} + \frac{z}{4 \lambda z^2 \omega _{02}^2} \left( {\dot{Q}}_{ij} {\dot{Q}}^{ij} -z^2\omega _{02}^2Q_{ij} Q^{ij} \right) + L_\text {FD} \right] .\nonumber \\ \end{aligned}$$
(27)

Here, the term \(L_{{\mathrm{FD}}}\) describes relativistic frame-dragging effects that arise kinematically when expressing \((D Q^{\mu \nu }/d \sigma )^2\) in terms of \( (d{Q}_{ij}/dt)^2 \). The term \(L_\text {FD}\) describes the coupling of the orbital angular momentum to the angular momentum (or spin) associated with the quadrupole

$$\begin{aligned} S_Q^{i}=\frac{1}{\lambda \omega _{02}^2} \epsilon _{ijm}\left[ Q^{kj} {\dot{Q}}^{m}{}_k-Q^{km} {\dot{Q}}^{j}{}_k\right] . \end{aligned}$$
(28)

To proceed with computing the dynamics from (27) and the GWs requires explicit expressions for the various quantities such as \({{{\mathcal {E}}}}_{ij}\) and z and the frame-dragging terms appearing in the action. These have been computed in post-Newtonian (PN) theory [78, 505, 543], see also [335] for the effective field theory description, in the test-particle limit [80] or the gravitational self-force formalism [76, 186, 408, 490]. An alternative approach based on an affine model has also been studied [211, 372].

4.3 Computation of tidal deformability

We have seen that the action summarizing the dynamics involves the characteristic information about matter properties through the 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 [210, 300]. 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.

4.3.1 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:

$$\begin{aligned}&\nabla ^2 U=-4\pi \rho , \qquad \frac{\partial \rho }{\partial t}+\nabla \cdot (\rho {{\varvec{v}}})=0, \nonumber \\&\frac{\partial v^i}{\partial t}+({{\varvec{v}}}\cdot \nabla )v^i=-\frac{\partial ^i p}{\rho }+\partial _i U +a^i_{{\mathrm{ext}}}. \end{aligned}$$
(29)

Here, \(a^i_{{\mathrm{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 \(\rho _0\), and gravitational potential \(U_0\). Next, we consider the effect of external static tidal field \({{{\mathcal {E}}}}_L\). This causes the NS matter to adjust to a new static configuration with an exterior gravitational potential given by

$$\begin{aligned} U_{{\mathrm{total}}}=U+U_{{\mathrm{ext}}}=\frac{m}{r}+\sum _{\ell =2}^\infty \sum _{m=-\ell }^\ell Y_{\ell m} \left[ \frac{Q_{\ell m}}{r^{\ell +1}} -\frac{1}{(2\ell -1)!!} {{{\mathcal {E}}}}_{\ell m} r^\ell \right] \end{aligned}$$
(30)

The linear response defining \(\lambda \) in (23) can equivalently be written as

$$\begin{aligned} Q_{\ell m}=-\lambda _\ell {{{\mathcal {E}}}}_{\ell m} \end{aligned}$$
(31)

for each \(\ell \)th multipole and azimuthal quantum number m. To obtain \(\lambda _\ell \), 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+\delta p\), \(\rho =\rho _0+\delta \rho \), \(U=U_0+\delta U\), and an external acceleration \({{\varvec{a}}}_{{\mathrm{ext}}}={\varvec{\nabla }} U_{{\mathrm{tidal}}}\). The fluid perturbation can be represented by a Lagrangian displacement \(\xi (x,t)\), which maps a fluid element at position x in the unperturbed NS to \(x+\xi (x,t)\) in the perturbed star. Expanding the Euler equation to linear order in the perturbations leads to the differential equation

$$\begin{aligned} \frac{d^2 \xi ^i}{dt^2}=-\frac{\partial _i \delta p}{\rho _0}+\frac{\partial _i p_0}{\rho _0^2}\delta \rho +\partial _i \delta U+a^i_{{\mathrm{ext}}}. \end{aligned}$$
(32)

This expression is used in general for calculating stellar oscillations. Here we are interested in specializing to static perturbations with \(({{\dot{\xi }}}=0)\). For a barotropic EOS relation \(p=p(\rho )\) we can eliminate \(\delta p\) from (32) in favor of \(\delta \rho \), combine all terms that involve \(\delta \rho \) into a total derivative and integrate. We then decompose the perturbations as \(\delta \rho =f(r)Y_{\ell m}(\theta ,\phi )\) and \(\delta U_{{\mathrm{tot}}}=H(r)Y_{\ell m}(\theta ,\phi )\), substitute into the linearized Poisson’s equation, and after manipulations obtain the differential equation for H(r) in the region \(r\le R\):

$$\begin{aligned} H^{\prime \prime }+\frac{2}{r}H^\prime -\frac{\ell (\ell +1)}{r^2}H =-{4\pi }\left( \frac{1}{\rho _0}\frac{dp_0}{d\rho _0}\right) ^{-1}H. \end{aligned}$$
(33)

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\propto r^\ell \) for \(r\rightarrow 0\). For \(r>R\) the exterior solution is

$$\begin{aligned} H^{{\mathrm{ext}}}=\frac{Q_{\ell m}}{r^{\ell +1}} -\frac{1}{(2\ell -1)!!} {{{\mathcal {E}}}}_{\ell m} r^\ell =-{{{\mathcal {E}}}}_{\ell m}\left[ \frac{\lambda _\ell }{r^{\ell +1}} +\frac{1}{(2\ell -1)!!} r^\ell \right] . \end{aligned}$$
(34)

To extract the Love numbers we eliminate \({{{\mathcal {E}}}}_{\ell m}\) by considering the logarithmic derivative

$$\begin{aligned} y(r)=\frac{r H'(r)}{H(r)}. \end{aligned}$$
(35)

We solve for \(k_\ell \) by using (34) in (35) and matching the result for y(R) obtained from the interior and exterior solutions at the star’s surface:

$$\begin{aligned} k_\ell =\frac{\ell -y(R)}{2\left[ \ell +1+y(R)\right] |} \end{aligned}$$
(36)

For an incompressible \(n=0\) polytrope the exact solution is \(k_\ell =3/(4\ell -4)\).

4.3.2 Calculation of tidal deformability in GR

The definition of \(\lambda _\ell \) from Eq. (23) is general and holds with the relativistic definitions of \({{{\mathcal {E}}}}_L\) and \(Q_L\). The computation of \(\lambda _\ell \) 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.

For the equilibrium solution we use the metric of a spherically symmetric, isolated nonspinning NS solution expressed in the form [258]

$$\begin{aligned} ds_0^2 =g_{\mu \nu }^{(0)}dx^\mu dx^\nu = -e^{\nu (r)} dt^2 + e^{\gamma (r)} dr^2 + r^2 (d\theta ^2 + \sin ^2 \theta d\varphi ^2) . \end{aligned}$$
(37)

The NS matter is described by the stress–energy tensor \(T_{\mu \nu } = (\rho + p) u_{\mu } u_{\nu } + p g_{\mu \nu }\), where \(u^\mu \) is the fluid’s four-velocity and \(\rho \) the energy density. In the NS’s rest frame, the normalization condition \(u_\mu u^\mu = -1\) implies that \(u^{\mu } = (e^{-\nu /2}, 0, 0,0)\). The Einstein field equations and energy–momentum conservation \(\nabla _\nu T^{\mu \nu }=0\) determine the NS structure

$$\begin{aligned} G_{\mu \nu }=8\pi T_{\mu \nu } , \qquad \nabla _\nu T^{\mu \nu }=0. \end{aligned}$$
(38)

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_{\mu \nu }=g_{\mu \nu }^{(0)}+h_{\mu \nu }\). For our purposes we can decompose the perturbation as [171, 279, 462, 520]

$$\begin{aligned} h_{\mu \nu }dx^\mu dx^\nu =\sum _{\ell ,m}\bigg [-e^{\nu }H_0^{\ell m}dt^2+2H_1^{\ell m}dtdr+e^\gamma H_2^{\ell m}dr^2+r^2K^{\ell m}d\varOmega ^2\bigg ]Y^{\ell m}\nonumber \\ \end{aligned}$$
(39)

To compute tidal Love numbers we consider static perturbations where \(H_0\), \(H_2\), and K depend only on r. The perturbations to the stress–energy tensor are likewise decomposed into spherical harmonics. Substituting these decompositions into Eqs. (38), and extracting only the pieces that are linear in the perturbations in all the components leads to various relations between the different perturbation functions. Ultimately this yields a master equation for \(H\equiv H_0^{\ell m}\) [81, 155, 266]

$$\begin{aligned} 0= & {} \frac{d^2 H}{dr^2} + \left\{ \frac{2}{r} + e^{\gamma } \left[ \frac{2M}{r^2} + 4 \pi r (p-\rho ) \right] \right\} \frac{dH}{dr} \nonumber \\&+ \left\{ e^\gamma \left[ - \frac{\ell (\ell +1)}{r^2} + 4 \pi (\rho + p) \frac{d\rho }{dp} + 4 \pi (5 \rho + 9 p) \right] - \left( \frac{d\nu }{dr} \right) ^2 \right\} H .\qquad \quad \end{aligned}$$
(40)

The initial condition at the center is that \(H \propto r^\ell \) to ensure regularity of the solution.

Outside the star, the metric perturbation has the general form \(H_\ell = a_{\ell }^{Q} {Q}_{\ell 2}(x)+a_{\ell }^{P} {P}_{\ell 2}(x) \), where \(x \equiv r/M -1\) and \({P}_{\ell 2}(x)\) and \({Q}_{\ell 2}(x)\) are the normalized associated Legendre functions of the first and second kinds respectively. They are normalized such that for \(x \rightarrow \infty \), \({P}_{\ell 2}(x) \sim x^{\ell }\) and \({Q}_{\ell 2}(x) \sim x^{-(\ell + 1)}\). The constants \(a_{\ell }^{P}\) and \(a_{\ell }^{Q}\) are determined by matching the logarithmic derivative of the interior and exterior solutions,

$$\begin{aligned} y_\ell \equiv \frac{r}{H_\ell } \frac{dH_\ell }{dr} , \end{aligned}$$
(41)

at the NS surface [266]. Comparing with the definition of \(Q_\ell \) and \({{{\mathcal {E}}}}_\ell \) in the asymptotic metric (13) and the definition of the Love numbers shows that

$$\begin{aligned} (2\ell -1)!!\varLambda _\ell = -\left. \frac{{P}_{\ell 2}'(x)-C y_{\ell } {P}_{\ell 2}(x)}{{Q}_{\ell 2}'(x)- C y_{\ell } {Q}_{\ell 2}(x)}\right| _{x=1/C-1}, \end{aligned}$$
(42)

For the dominant quadrupolar effect the explicit expression is

$$\begin{aligned} \varLambda \mid _{\ell =2}&= \frac{16}{15} (1-2 C)^2[2+2C(y-1)-y] \nonumber \\&\quad \times \left\{ 2C[6-3y+3C(5y-8)] \right. \nonumber \\&\quad +4C^3[13-11y+C(3y-2)+2C^2(1+y)] \nonumber \\&\quad \left. +3(1-2C)^2[2-y+2C(y-1)] \ln (1-2C) \right\} ^{-1} , \end{aligned}$$
(43)

where y is evaluated at the surface \(r=R\). For an incompressible star with \(\rho =\text {const}\) or \(p=K \rho ^{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 [155]. From Eq. (41) in the interior we obtain \(y^{{\mathrm{in}}}(R)\). The density discontinuity introduces a correction to y just outside the star \(y^{{\mathrm{out}}}\) given by

$$\begin{aligned} y^{{\mathrm{out}}}_{{\mathrm{incompressible}}} = y^{{\mathrm{in}}}- 3. \end{aligned}$$
(44)

This is valid for any \(\ell \) and (44) is to be substituted for y when evaluating \(\varLambda \) from (42).

A number of studies have computed examples of \(\varLambda \) for a wide range of EOS models and delineated interpretations of its physics content such as symmetry energy, phase transitions, and multi-body interactions [16, 25, 29, 30, 59, 74, 86, 112, 113, 121, 129, 168, 204,205,206, 213, 223, 239, 246, 254, 267, 270, 282, 293, 305, 307, 319, 328, 329, 343,344,345,346,347,348,349,350,351, 353,354,355, 362, 363, 365, 378, 385, 402, 430, 431, 433,434,435, 446, 457, 458, 503, 526, 536, 546, 548, 556, 561, 565,566,572]. Tidal parameters can also distinguish exotic objects, matter halos, and modified gravity [49, 50, 60, 108,109,110, 115, 118, 166, 189, 196, 197, 241, 251, 318, 319, 374,375,376, 379, 380, 406, 407, 417,418,419, 443, 444, 448, 449, 488,489,491, 534, 535, 559, 560]. Substantial recent interest has been attracted by quasi-universal relations [552, 553] that link dimensionless parameters characterizing various global properties of the NS in an approximately EOS-independent way [24, 98, 116, 119, 120, 123, 128, 261, 309, 329, 371, 423, 427,428,429, 465, 480, 501, 505, 551, 552].

4.4 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 [27, 216, 270, 300, 317, 437, 448, 493, 528, 529, 561]. Effects from nonlinear mode coupling have also been studied [200, 323, 550]. In GR, besides the tidal couplings to the gravito-electric fields \({{{\mathcal {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 \( {{{\mathcal {B}}}}_L=3 \epsilon _{<a_1 jk}C^{jk\,}_{\; \; a_20; a_3 \cdots a_\ell >}/2(\ell +1)(\ell -2)!\) where \(\epsilon _{ijk}\) is the completely antisymmetric permutation pseudotensor. The induced current moments \({{{\mathcal {S}}}}_L\) appear in the time–space part of the asymptotic metric in a local asymptotic frame, see Refs. [208, 216, 323, 324, 420] for more details. Similar to the tidal deformability coefficients \(\lambda _\ell \), the relation between \({{{\mathcal {S}}}}_L\) and \({{{\mathcal {B}}}}_L\) is characterized by a set of gravitomagnetic Love numbers \(\sigma _\ell =- {{{\mathcal {S}}}}_L / {{{\mathcal {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 [15, 182, 234, 284, 321, 325, 421, 422] and can significantly impact the tidal excitation of quasi-normal modes [188, 230, 270, 361], in addition to the effects from spin-induced multipoles [307, 337, 338, 341, 342, 369, 441]. As discussed in the previous section, eccentricity also greatly enriches matter phenomena and mode excitations [127, 242, 540, 558]. 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 [212, 225, 291, 313, 371, 424, 497, 535], see Ref. [497] 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  [70, 203, 280, 373, 393, 438, 439, 441, 508, 512, 515].

4.5 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 [87]. The radiative post-Minkowski solution can be formally written in a radiative coordinate system \((u,{{\varvec{X}}})\) where \(u=T-|X|/c\) is the retarded time, as

$$\begin{aligned} h_{ij}^{{\mathrm{rad}}}(u, X)= & {} \frac{4G}{D c^2} {{{\mathcal {P}}}}_{ijab}(N) \sum _{\ell =2}^\infty \frac{1}{c^\ell \ell !}\bigg [N_{L-2}U_{abL-2}(u)\\&-\frac{2\ell }{c(\ell +1)}N_{dL-2}\epsilon _{de(a}V_{b)eL-2}(u)\bigg ] + O(D^{-2}) \nonumber \end{aligned}$$
(45)

Here, \({{{\mathcal {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 non-local in time, and complicated in general. At the leading Newtonian order, there is a simple relation to the source quadrupole \(U_{ij}^{{\mathrm{Newt}}}=d^2 Q^{{\mathrm{T}}}_{ij}/du^2\), where the total quadrupole moment of the system is the sum of the orbital and NS contributions \(Q_{ij}^{{\mathrm{T}}}=\mu 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\) [87].

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 lowest-order matter signatures in GWs. The energy of the binary is \(E=E_{{\mathrm{pm}}}+E_{{\mathrm{tidal}}}\), where the subscript “pm” denotes point-masses. Specializing to circular orbits \(\ddot{r}={\dot{r}}=0\) and \({{\dot{\phi }}}=\varOmega \), \(\ddot{\phi }=0\), and to perturbatively solving the dynamics described by the action (25), and eliminating r in terms of the frequency variable \(x=(M\varOmega )^{2/3}\) we find the following tidal contribution in the limit of adiabatic tides [543]

$$\begin{aligned} E_{{\mathrm{tidal}}}(x)= -\frac{1}{2}\mu x\left[ -9\frac{m_B}{m_A} \frac{\lambda _{\ell =2,A}x^5}{M^5}-65\frac{m_B}{m_A} \frac{\lambda _{\ell =3,A}x^7}{M^7}\right] +(A\leftrightarrow B). \end{aligned}$$
(46)

Here, the subscripts AB 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 [541]

$$\begin{aligned} P_{{\mathrm{GW}}}^{{\mathrm{tidal}}}=\frac{32\mu ^2}{5M^2} x^{5/2}\left[ \left( \frac{18M}{m_A}-12\right) \frac{\lambda _{\ell =2,A }x^5}{M^5}+\frac{80m_B}{m_A}\frac{\lambda _{\ell =3,A}x^7}{M^7}\right] . \end{aligned}$$
(47)

By requiring that \(P_{{\mathrm{GW}}}\) be balanced by a change in the energy E of the binary one can derive the evolution equations

$$\begin{aligned} \frac{d\phi }{dt}=\frac{x^{3/2}}{M}\ \qquad \qquad \frac{dx}{dt}=\frac{-{P}_{{\mathrm{GW}}}}{dE/dx} \end{aligned}$$
(48)

There are several ways to solve for \(\phi \) in a PN approximation. For example, one can numerically solve Eqs. (48) for \(\phi (t)\) and x(t) after first expanding the ratio \({P}_{{\mathrm{GW}}}/({dE/dx})\) about \(x=0\) to the consistent PN order. These waveforms are known as TaylorT4 approximants as review in Ref. [104], where the point-mass terms are also given explicitly. The adiabatic quadrupolar tidal corrections that add linearly to the point-mass contributions:

$$\begin{aligned} \frac{d{x}}{dt}\mid _{{\mathrm{tidal}}}=\frac{32}{5}\frac{m_B}{M^7}\lambda _{\ell =2,A} \, x^{10}\left[ 12\left( 1+11\frac{m_B}{M}\right) \right] +(A\leftrightarrow B). \end{aligned}$$
(49)

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. [544]. 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. [141] and leads to a signal of the form

$$\begin{aligned} {\tilde{h}}(f)={{{\mathcal {A}}}}f^{-7/6}\exp \left[ i\left( \psi _{{\mathrm{pm}}}+\psi _{{\mathrm{tidal}}}\right) \right] , \end{aligned}$$
(50)

where f is the GW frequency, \({{{\mathcal {A}}}}\propto {{{\mathcal {M}}}}^{5/6}/D\), where \({{{\mathcal {M}}}}\) is the chirp mass \({{{\mathcal {M}}}}=\nu ^{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 \({{{\mathcal {A}}}}\), where higher PN order corrections to the amplitude enter. The point-mass phase \(\psi _{{\mathrm{pm}}}\) to the current best knowledge for nonspinning binaries is given e.g. in Eq. (3.18) of [104]. The leading order tidal contributions to \(\psi _{{\mathrm{tidal}}}\) can be expressed as:

$$\begin{aligned} \delta \psi _{{\mathrm{tidal}}}=\frac{3}{128\nu x^{5/2}}\left[ -\frac{39}{2}{\tilde{\varLambda }}x^5-\frac{4000 }{9}\frac{x^7}{M^7}\left( \varLambda _{\ell =3,A} m_A^6 m_B+\varLambda _{\ell =3, B} m_B^6 m_A\right) \right] ,\nonumber \\ \end{aligned}$$
(51a)

where

$$\begin{aligned} {\tilde{\varLambda }}=\frac{16}{13}\left[ \left( 1+\frac{12 m_B}{m_A}\right) \frac{ \varLambda _{\ell =2,A}m_A^5}{M^5}+\left( 1+\frac{12 m_A}{m_B}\right) \frac{ \varLambda _{\ell =2,B}m_B^5}{M^5}\right] \end{aligned}$$
(51b)

The parameter \({{\tilde{\varLambda }}}\) plays an analogous role as the chirp mass. For equal-mass BNSs \({\tilde{\varLambda }}\) reduces to \(\varLambda \) of the individual neutron stars. Other combinations of the two parameters \(\varLambda _A\) and \(\varLambda _B\) are also in use and have advantages in different contexts, e.g. to characterize the dominant effect in the conservative dynamics [156], or to improve the measurability [554, 555]. 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 mode-frequency-dependent and regularized contributions [483]:

$$\begin{aligned} \delta \psi _{{\mathrm{tidal}}}^{{\mathrm{f-modes}}}\approx & {} -\frac{\left[ 10 \sqrt{3}\pi -27-30\log (2)\right] }{96\nu }\frac{\varLambda _{\ell =2,A}m_A^4}{M^6 \omega _{02,A}^2}\left( 155-147\frac{m_A}{M}\right) x^{11/2}\nonumber \\&-\frac{1875\left[ 5-6\log (2)\right] }{16}\frac{\varLambda _{\ell =3,A}m_A^5}{M^7 \omega _{03,A}^2}x^{15/2} \end{aligned}$$
(52)

Finally, we give the full analytical knowledge to the highest complete PN order \(O(\alpha ^5 \epsilon ^{3/2})\), i.e., 1.5PN corrections to tidal effects for the quadrupolar adiabatic phasing [164]:

$$\begin{aligned} \delta \psi _{{\mathrm{tidal}}} = - \varLambda _{\ell =2,A} c^A_{{\mathrm{Newt}}} x^{5/2} \left( 1 + c^A_1 x + c^A_{3/2} x^{3/2} \right) + A \leftrightarrow B , \end{aligned}$$
(53)

with the dimensionless EOB tidal parameter \(\kappa _A\) (defined below) and \(x(M \omega )=(M \omega /2)^{2/3}\). The individual coefficients \(c_i^A\) are

$$\begin{aligned} c^A_{{\mathrm{Newt}}}= & {} -\frac{3m_A^3}{8M^3} \left( \frac{m_A}{m_B}+12\right) , \end{aligned}$$
(54a)
$$\begin{aligned} c^A_1= & {} -\frac{5 \left( 260 \frac{m_A^3}{M^3}-2286 \frac{m_A^2}{M^2}-919 \frac{m_A}{M}+3179\right) }{336 \left( 11 \frac{m_A}{M}-12\right) } , \end{aligned}$$
(54b)
$$\begin{aligned} c^A_{3/2}= & {} -\frac{5}{2} \pi , \end{aligned}$$
(54c)

and similarly with \(A \leftrightarrow B\). Partial knowledge is available at higher orders. The \(O(\epsilon ^{5/2})\) terms are known, however, the lower order \(O(\epsilon ^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):

$$\begin{aligned} c_2^A= & {} \left[ 5 \left( 67702048 \frac{m_A^5}{M^5}-223216640\frac{m_A^4}{M^4}+337457524 \frac{m_A^3}{M^3} \right. \right. \nonumber \\&\left. \left. -141992280 \frac{m_A^2}{M^2}+96008669 \frac{m_A}{M}-143740242\right) \right] / \nonumber \\&\left[ 3048192 \left( 11 \frac{m_A}{M}-12\right) \right] , \end{aligned}$$
(55)
$$\begin{aligned} c_{5/2}^A= & {} -\frac{\pi \left( 10232 \frac{m_A^3}{M^3}-7022 \frac{m_A^2}{M^2}+22127 \frac{m_A}{M}-27719\right) }{192 \left( 11 \frac{m_A}{M}-12\right) } \end{aligned}$$
(56)

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. [93, 103, 307, 338, 339, 368, 444, 445]. The leading-order contribution to the GW phase from the spin-induced quadrupole is given by [306]:

$$\begin{aligned} \delta \psi ^{{\mathrm{spin-quad}}}=\frac{30}{128\nu ^{4/5}(\pi {{{\mathcal {M}}}}f)^{1/3}} \left[ -5\kappa _A\chi _A^2\frac{m_A^2}{M^2}+(A\leftrightarrow B)\right] \end{aligned}$$
(57)

Figure 9 illustrates the various contributions to the GW phase discussed above.

Fig. 9
figure 9

Example of matter effects in the GW phasing for an inspiral starting from 20 Hz. The waveforms were aligned over a window of 30 Hz in the beginning and the parameters correspond to the DD2 EOS. Depending on parameters, the spin-quadrupole effects (red curves) can dominate over tidal effects (blue plus cyan curves), especially at low frequencies

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=(\pi M f)^{2/3}\), which corresponds to \(O(\alpha ^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 Sect. 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 [101, 102] 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.

5 Tidal effective-one-body models for binary inspirals

The EOB approach [101, 102] (see also the review [105, 144, 157, 160]) combines results from the PN theory, valid for arbitrary mass ratios, with information about the strong-field effects in the test-particle limit. In addition, the model makes use of a variety of other theoretical inputs, for example from scattering and gravitational self-force calculations. 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 test-mass 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 \(\nu \in [0,1/4]\) and thus deform the spacetime away from those of a BH. In addition, non-geodesic corrections that depend on higher-than-quadratic 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 [31, 32, 42,43,44, 77].

The dissipative sector of current EOB models is based on a concise factorized representation of the gravitational waveforms  [69, 149, 153, 381, 394, 399, 413, 415, 513, 514] 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. [33, 45,46,47, 94, 134, 412, 414, 415, 513, 514], and one by T. Damour, A. Nagar, and their collaborators, e.g. [40, 70, 143, 149, 150, 154, 158, 158, 159, 161,162,163, 381, 391, 393, 396, 398, 399]. 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 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, 36, 67, 68, 76, 78, 156, 269, 400, 505]. 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 [314, 395]. 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).

Fig. 10
figure 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

5.1 Brief summary of EOB models for nonspinning binaries

The conservative dynamics within the EOB approach are described by the Hamiltonian

$$\begin{aligned} H_{{\mathrm{EOB}}}=M\sqrt{1+2\nu \left( H_{{\mathrm{eff}}}-1\right) }-M, \end{aligned}$$
(58)

where the effective Hamiltonian \(H_{{\mathrm{eff}}}\) is that of a particle of mass \(\mu \). This energy mapping also emerges in quantum-electrodynamics [99], scattering theory [145, 146, 541, 543], and from explicit comparisons of the action variables [101]. For nonspinning binaries the effective Hamiltonian is given by

$$\begin{aligned} H_{{\mathrm{eff}}}^2=p_{r*}^2+A\left( 1+\frac{p_{\phi }^2}{r^2}+2 \left( 4-3\nu \right) \frac{p_{r*}^4}{\nu r^2}\right) . \end{aligned}$$
(59)

Here, \(p_\phi \) is the azimuthal angular momentum and \(p_{r*}=p_r/\sqrt{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_{{\mathrm{pm}}}+A_{{\mathrm{tidal}}}\) and \(D=D_{{\mathrm{pm}}}\). The different choices for the point-mass potentials \(A_{{\mathrm{pp}}}\) are either Eqs. (A1)–(A2h) of Ref. [505] with the calibration parameter determined most recently in [94] or the resummation from [22, 67, 400]. The potential D is either taken from Eq. (A4) of Ref. [505] or [400].

The EOB equations of motion are

$$\begin{aligned} \begin{aligned} \frac{d r}{dt}&=\frac{A}{\sqrt{D}} \frac{\partial H_{{\mathrm{EOB}}}}{\partial p_{r*}}, \qquad \frac{d p_{r*}}{dt}=-\frac{A}{\sqrt{D}} \frac{\partial H_{{\mathrm{EOB}}}}{\partial r}+{{{\mathcal {F}}}}_r , \\ \frac{d \phi }{dt}&=\frac{\partial H_{{\mathrm{EOB}}}}{\partial p_\phi }, \qquad \frac{d p_\phi }{dt}= {{{\mathcal {F}}}}_\phi , \end{aligned} \end{aligned}$$
(60a)

The factors of \(A/\sqrt{D}\) appear because \(p_{r*}\) and r are not canonically conjugate variables. The gravitational radiation reaction forces \({{{\mathcal {F}}}}_\phi \) and \({{{\mathcal {F}}}}_r\) are constructed from

$$\begin{aligned} {{{\mathcal {F}}}}_\phi =-\frac{1}{\varOmega _\phi }\mathcal {\dot{E}}_{{\mathrm{rad}}}, \quad {{{\mathcal {F}}}}_r=\frac{p_{r*}}{p_\phi }{{{\mathcal {F}}}}_\phi , \end{aligned}$$
(61)

where the energy flux \(\mathcal {\dot{E}}_{{\mathrm{rad}}}\) is given by

$$\begin{aligned} \mathcal {\dot{E}}_{{\mathrm{rad}}} =\frac{\varOmega _\phi ^2}{8\pi } \sum _{\ell =2}^8\sum _{m=0}^\ell m^2 |h_{\ell m}^{{\mathrm{F}}}|^2. \end{aligned}$$
(62)

The sum is only over positive m since \(|h_{\ell -m}^{{\mathrm{F}}}|=|h_{\ell m}^{{\mathrm{F}}}|\). The factorized EOB waveforms for point masses have the form [148, 149, 153]

$$\begin{aligned} {h_{\ell m}^{{\mathrm{F}}}}_{{\mathrm{pp}}}=h_{\ell m}^{{\mathrm{Newt}}}S_{\ell m}T_{\ell m}\rho _{\ell m}^\ell N_{\ell m}. \end{aligned}$$
(63)

Each of these terms depends on the EOB canonical variables and the parameters of the binary, see e.g. [94, 400] for the explicit up to date expressions.

5.1.1 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 [78, 156, 543]

$$\begin{aligned} A^{{\text {tidal}}}= & {} -{{{\mathcal {A}}}}_{\ell =2}^{(A)} {\hat{A}}^{(A)} -{{{\mathcal {A}}}}_{\ell =3}^{(A)} \bigg [ 1 + \left( \frac{15}{2}\frac{m_A}{M}-2\right) r^{-1} \nonumber \\&+ \left( \frac{110}{3}\frac{m_A^2}{M^2}-\frac{311}{24}\frac{m_A}{M}+\frac{8}{3}\right) r^{-2} \bigg ] +(A \leftrightarrow B). \end{aligned}$$
(64)

Here \({{{\mathcal {A}}}}_\ell ^{(A)} \) are the Newtonian \(\ell \)th multipolar tidal potentials given by

$$\begin{aligned} {{{\mathcal {A}}}}_\ell ^{(A)}=\frac{2 m_B}{m_A}k_\ell ^{(A)}\frac{R_A^{2\ell +1}}{r^{2\ell +2}}. \end{aligned}$$
(65)

The 2PN corrections to the tidal potential \({\hat{A}}^{(A)}\) from Eqs. (6.6) and (6.18) of Ref. [78] are given by

$$\begin{aligned} {\hat{A}}_{ {\mathrm{adPN}}}^{(A)}= 1 + \frac{5m_A}{2M} r^{-1} + \left( \frac{337}{28}\frac{m_A^2}{M^2}+\frac{1}{8}\frac{m_A}{M}+3\right) r^{-2} \end{aligned}$$
(66)

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_{\ell m}^{{\mathrm{tidal}}}\) to the GW modes so that

$$\begin{aligned} h_{\ell m}^{{\mathrm{F}}}={h_{\ell m}^{{\mathrm{F}}}}_{{\mathrm{pp}}}+h_{\ell m}^{{\mathrm{tidal}}}. \end{aligned}$$
(67)

The explicit results for the adiabatic tidal terms \(h_{\ell m}^{{\mathrm{tidal}}}\) are given in Eqs. (A14)–(A17) of Ref. [164] for \(\ell \le 3\).

5.1.2 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. [269, 505]. 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_\ell \) multiplied by a frequency-dependent enhancement factor such that \(k_\ell \rightarrow k_\ell {\hat{k}}^ {\mathrm{dyn}}_{\ell }\) with

$$\begin{aligned} {\hat{k}}^{{\mathrm{dyn}}}_{\ell }{=}a_\ell {+}b_\ell \bigg [ \frac{\omega _{0\ell }^2 }{\omega _{0\ell }^2{-}(m\varOmega ^2)}{+}\frac{\omega _{0\ell }^2}{2\sqrt{\epsilon _m}{{\hat{t}}} \varOmega ^\prime (m\varOmega )^2}{+}\frac{\sqrt{\pi }\omega _{0\ell }^2}{\sqrt{3}\sqrt{\epsilon _m}(m\varOmega )^2} {{{\mathcal {Q}}}}_{\ell m} \bigg ] \end{aligned}$$
(68)

with \(\varOmega ^\prime =3/8\), and \(\varOmega =M^{1/2}r^{-3/2}\). The quantity \({{{\mathcal {Q}}}}_{\ell m}\) is given by

$$\begin{aligned} {{{\mathcal {Q}}}}_{\ell m}=\cos (\varOmega ^\prime {{\hat{t}}}^2)\left[ 1+2 {\mathrm{F}}_{{\mathrm{S}}}\left( \frac{\sqrt{3}}{\sqrt{4\pi }}{{\hat{t}}} \right) \right] -\sin (\varOmega ^\prime {{\hat{t}}}^2)\left[ 1+2 {\mathrm{F}}_{{\mathrm{C}}}\left( \frac{\sqrt{3}}{\sqrt{4\pi }}{{\hat{t}}} \right) \right] , \end{aligned}$$
(69)

where the functions \({\mathrm{F}}_{{\mathrm{S}}}\) and \({\mathrm{F}}_{{\mathrm{C}}}\) are Fresnel sine and cosine integrals respectively using the conventions in Mathematica. The quantities \({{\hat{t}}}\) and \(\epsilon _m\) are defined as

$$\begin{aligned} {{\hat{t}}}=\frac{8}{5\sqrt{\epsilon _m}}\left( 1-\frac{r^{5/2} \omega _{0\ell }^{5/3}}{m^{5/3}M^{5/6}}\right) , \qquad \epsilon _m=\frac{256 \mu M^{2/3}\omega _{0\ell }^{5/3}}{5m^{5/3}}. \end{aligned}$$
(70)

In Eq. (68) a body label AB on the quantities \(\omega _{0\ell }\), \({{\hat{t}}}\), \(\epsilon _m\), and \({{{\mathcal {Q}}}}_{\ell m}\) is implied. For each \(\ell \)-multipole only \(m=\ell \) contributes in Eq. (68) because the effect of modes with \(m<\ell \) has already been taken into account as adiabatic contributions. For the quadrupole and octupole multipole moments the coefficients are given by \((a_2, a_3)=(\frac{1}{4},\frac{3}{8})\) and \((b_2,b_3)=(\frac{3}{4}, \frac{5}{8})\).

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_\ell \) in \(h_{\ell m}^{{\mathrm{tidal}}}\) in Eq. (67) by an effective function \(k_\ell \rightarrow k_\ell \ {\hat{k}}_{\ell \ {\mathrm{dyn}}}^{{\mathrm{diss}}}\). For \(\ell =2\) this function is given by

$$\begin{aligned} {\hat{k}}_{2\ {\mathrm{dyn}}}^{(A)\ {\mathrm{diss}}}=\frac{\omega _{02}^2+{\hat{k}}_{2\ {\mathrm{dyn}}}^{(A)}\left( \omega _{02}^2+6 \frac{m_B}{M}\varOmega ^2\right) }{3\varOmega ^2\left( 1+2 \frac{m_B}{M}\right) }, \end{aligned}$$
(71)

where \({\hat{k}}_{2\ {\mathrm{dyn}}}^{(A)}\) is the enhancement function for body A in the conservative dynamics from Eq. (68).

5.1.3 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\) [76, 186]. These results have been augmented in Refs. [22, 67, 76, 400] by a term \(\propto m_A^2/M^2\) that would describe currently unknown second-order self-force effects. Specifically, in this model Eq. (64) is employed with the potential \({\hat{A}}^{(A)}\) given by

$$\begin{aligned} {\hat{A}}_{ {\mathrm{adGSF^+}}}^{(A)}= & {} 1+\frac{3r^{-2}}{\left( 1-r_{{\mathrm{LR}}}/r\right) } +\frac{m_A}{M}\frac{a_1^{{\mathrm{GSF}}}(r^{-1}) }{ \left( 1-r_{{\mathrm{LR}}}/r\right) ^{7/2}}\nonumber \\&+\frac{m_A^2}{M^2} \frac{a_2^{{\mathrm{2GSF}}}(r^{-1}) }{ \left( 1-r_{{\mathrm{LR}}}/r\right) ^{p}} +(A\leftrightarrow B) \end{aligned}$$
(72)

In Eq. (72), the coefficient \(a_1^{{\mathrm{GSF}}}(r^{-1}) \) is obtained from Eqs. (7.24)–(7.27) of Ref. [76]. In the model of Refs. [22, 67, 400] the choices for the unknown parameters are \(a_2^{{\mathrm{2GSF}}}=337r^{-2}/28\) and \(p=4\). The radius of the light ring \(r_{{\mathrm{LR}}}\) is obtained from the conservative EOB dynamics by solving

$$\begin{aligned} {\tilde{A}}(r_{{\mathrm{LR}}})-\frac{r_{{\mathrm{LR}}}}{2}{\tilde{A}}^\prime (r_{{\mathrm{LR}}})=0, \end{aligned}$$
(73)

with the potential \({\tilde{A}}=A_{{\mathrm{pp}}}+A^{{\mathrm{tidal}}}_{{\mathrm{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 point-particle sector such as the factor \(N_{\ell m}\) in the waveform modes of Eq. (63), a different resummation of the potentials \(A_{{\mathrm{pm}}}\) and D, and the arguments of factors in \(h_{\ell m}^{{\mathrm{F}}}\) which involve powers of \(v_\phi =(\partial H_{{\mathrm{EOB}}}/\partial p_\phi )^{-2/3}\varOmega \) are evaluated for circular orbits instead of \(v=\varOmega ^{1/3}\), Padé-resummation of higher-order hereditary terms in \(h_{\ell m}^{{\mathrm{F}}}\), and a different calibration of the BBH model to NR.

6 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

$$\begin{aligned} \phi (M \omega ) = \phi _{{\mathrm{pm}}} (M \omega ) + \phi _{{\mathrm{SO}}}(M \omega ) + \phi _{{\mathrm{SS}}}(M \omega ) + \phi _{{\mathrm{T}}}(M \omega ) + \cdots \ , \end{aligned}$$
(74)

where \(\phi _{{\mathrm{pm}}}\) denotes the point-mass, \(\phi _{{\mathrm{SO}}}\) the spin-orbit, \(\phi _{{\mathrm{SS}}}\) the spin-spin, and \(\phi _{{\mathrm{T}}}\) the tidal phase contribution. We will focus in the following on an accurate representation of the tidal phase \(\phi _{{\mathrm{T}}}\) that is required for a proper modelling of BNS systems. However, in addition and as pointed out in Sect. 2, also the spin–spin interaction, \(\phi _{{\mathrm{SS}}}\), or higher order spin-contributions incorporate information about the EOS due to the presence of the quadrupole moment. Therefore, \(\phi _{{\mathrm{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

$$\begin{aligned} \psi (M \omega ) = \psi _{{\mathrm{pm}}} (M \omega ) + \psi _{{\mathrm{SO}}}(M \omega ) + \psi _{{\mathrm{SS}}}(M \omega ) + \psi _{{\mathrm{T}}}(M \omega )+ \cdots \ . \end{aligned}$$
(75)

In the following subsections, we will discuss different phenomenological approaches to describe \(\phi _{{\mathrm{T}}}\) and \(\psi _{{\mathrm{T}}}\).

6.1 NRTidal-approximants

We start our discussion with the newest NRTidal version (NRTidalv2) described in Ref. [182]. Earlier versions use similar principles and assumptions and therefore will not be discussed to avoid repetition, but we refer to [173, 175, 184] for more details.

The main idea of the NRTidal [175, 184] approach is to provide a closed-form approximation for the tidal phase contribution \(\phi _{{\mathrm{T}}}\) and \(\psi _{{\mathrm{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 [182]. 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. Sect. 4. Since the current NRTidal models use only leading and next-to-leading order mass-ratio effects, we restrict the parameters \(c_1^{A,B},c_{3/2}^{A,B},c_2^{A,B},c_{5/2}^{A,B}\) in Eqs. (54c) to their equal-mass values. Thus, the superscripts A and B are discarded in the following.

Under this assumption, an effective representation of tidal effects during the inspiral is obtained using

$$\begin{aligned} \phi _T (x) = - \kappa _{{\mathrm{eff}}}^T \frac{13}{ 8 \nu } x^{5/2} P_{{\mathrm{NRTidalv2}}}(x) \ , \end{aligned}$$
(76)

where \(P_{{\mathrm{NRTidalv2}}}(x)\) is given by the Padé approximant

$$\begin{aligned} P_{{\mathrm{NRTidalv2}}}(x) = \frac{1 + n_1 x+ n_{3/2} x^{3/2} + n_2 x^2 + n_{5/2} x^{5/2} + n_3 x^3}{1+ d_1 x+ d_{3/2} x^{3/2}+ d_2 x^2} \end{aligned}$$
(77)

and

$$\begin{aligned} \kappa _{{\mathrm{eff}}}^T = \frac{3}{16} {\tilde{\varLambda }}. \end{aligned}$$
(78)

To enforce consistency with the analytic PN knowledge [Eqs. (53)–(54c)], some of the parameters in Eq. (77) are fixed:

$$\begin{aligned} n_1= & {} c_1 + d_1, \end{aligned}$$
(79a)
$$\begin{aligned} n_{3/2}= & {} \frac{c_1 c_{3/2}-c_{5/2}-c_{3/2} d_1 + n_{5/2}}{c_1}, \end{aligned}$$
(79b)
$$\begin{aligned} n_{2}= & {} c_2 + c_1 d_1 + d_2 , \end{aligned}$$
(79c)
$$\begin{aligned} d_{3/2}= & {} - \frac{c_{5/2}+c_{3/2} d_1 - n_{5/2}}{c_1}. \end{aligned}$$
(79d)

with

$$\begin{aligned} c_{1}&= \frac{3115}{624}, \qquad c_{3/2} = -\frac{5 \pi }{2}, \end{aligned}$$
(80a)
$$\begin{aligned} c_2&= \frac{28024205}{1100736}, \qquad c_{5/2} = -\frac{4283 \pi }{312}. \end{aligned}$$
(80b)

The remaining parameters are fitted to high-resolution NR data

$$\begin{aligned} n_{5/2}&= 312.48173, \qquad n_3 = -342.15498, \end{aligned}$$
(81a)
$$\begin{aligned} d_1&= -20.237200, \qquad d_2 = -5.361630; \end{aligned}$$
(81b)

cf. Sect. 3.

To obtain the tidal corrections in the frequency domain, the NRTidal approach uses the stationary phase approximation (SPA); e.g., [164],

$$\begin{aligned} \frac{\text {d}^2 \psi _T(\omega )}{\text {d} \omega ^2} = \frac{1}{\omega } \frac{\text {d} \phi _T(\omega )}{\text {d} \omega } . \end{aligned}$$
(82)

Eq. (82) is solved numerically and the final result is again approximated by a Padé approximant

$$\begin{aligned} \psi _T (x) = - \kappa _{{\mathrm{eff}}}^T \frac{39}{ 16 \nu } x^{5/2} {\tilde{P}}_{{\mathrm{NRTidalv2}}}(x) \end{aligned}$$
(83)

with

$$\begin{aligned} {\tilde{P}}_{{\mathrm{NRTidalv2}}}(x) = \frac{1 + {\tilde{n}}_1 x+ {\tilde{n}}_{3/2} x^{3/2} + {\tilde{n}}_2 x^2+ {\tilde{n}}_{5/2} x^{5/2} + {\tilde{n}}_3 x^3}{1+ {\tilde{d}}_1 x+ {\tilde{d}}_{3/2} x^{3/2}+ {\tilde{d}}_2 x^2} . \end{aligned}$$
(84)

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):

$$\begin{aligned} {\tilde{n}}_1= & {} {\tilde{c}}_1 + {\tilde{d}}_1, \end{aligned}$$
(85a)
$$\begin{aligned} {\tilde{n}}_{3/2}= & {} \frac{{\tilde{c}}_1 {\tilde{c}}_{3/2}-{\tilde{c}}_{5/2}- {\tilde{c}}_{3/2} {\tilde{d}}_1 + {\tilde{n}}_{5/2}}{{\tilde{c}}_1}, \end{aligned}$$
(85b)
$$\begin{aligned} {\tilde{n}}_{2}= & {} {\tilde{c}}_2 + {\tilde{c}}_1 {\tilde{d}}_1 + {\tilde{d}}_2 , \end{aligned}$$
(85c)
$$\begin{aligned} {\tilde{d}}_{3/2}= & {} - \frac{{\tilde{c}}_{5/2}+{\tilde{c}}_{3/2} {\tilde{d}}_1 - {\tilde{n}}_{5/2}}{{\tilde{c}}_1}, \end{aligned}$$
(85d)

with

$$\begin{aligned} {\tilde{c}}_{1}&= \frac{3115}{1248}, \qquad {\tilde{c}}_{3/2} = -\pi , \end{aligned}$$
(86a)
$$\begin{aligned} {\tilde{c}}_2&= \frac{28024205}{3302208}, \qquad {\tilde{c}}_{5/2} = - \frac{4283 \pi }{1092}. \end{aligned}$$
(86b)

The remaining, fitting coefficients are

$$\begin{aligned} {\tilde{n}}_{5/2}&= 90.550822, \qquad {\tilde{n}}_3 = -60.253578, \end{aligned}$$
(87a)
$$\begin{aligned} {\tilde{d}}_1&= -15.111208, \qquad {\tilde{d}}_2 = 8.0641096. \end{aligned}$$
(87b)

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 [5, 485] or in determining if a GW source is gravitationally lensed [416].

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:

$$\begin{aligned} {\tilde{A}}_T^{{\mathrm{NRTidalv2}}} = - \sqrt{\frac{5 \pi \nu }{24}} \frac{9 M^2 }{D_L} \kappa ^T_{{\mathrm{eff}}} x^{13/4} \frac{1 + \frac{449}{108} x + \frac{22672}{9} x^{2.89}}{1+d \ x^4}, \end{aligned}$$
(88)

with \(d= 13477.8\).

To truncate the waveform after the moment of merger, a Planck taper [377] is added at the merger frequency

$$\begin{aligned} \omega _{{\mathrm{mrg}}} = \omega _0 \sqrt{\frac{1}{q}} \frac{1+ n_1 \kappa _{{\mathrm{eff}}}^T + n_2 (\kappa _{{\mathrm{eff}}}^T)^2}{1+d_1 \kappa _{{\mathrm{eff}}}^T + d_2 (\kappa _{{\mathrm{eff}}}^T)^2}, \end{aligned}$$
(89)

with \(n_1=3.354 \times 10^{-2}, n_2 = 4.315 \times 10^{-5}, d_1 = 7.542 \times 10^{-2}, d_2 = 2.236 \times 10^{-4}\). The parameter \(\omega _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, \(\omega _{{\mathrm{mrg}}}\), depending on the exact implementation and NRTidal model version. The final amplitude is given by

$$\begin{aligned} {\tilde{A}} = ({\tilde{A}}_{{\mathrm{BBH}}} + {\tilde{A}}_T^{{\mathrm{NRTidalv2}}}) \ {\tilde{A}}_{{\mathrm{Planck}}}; \end{aligned}$$
(90)

we note again that \({\tilde{A}}_T^{{\mathrm{NRTidal}}}=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. [190].

6.2 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. (i)

    Kawaguchi et al. employ hybrid waveforms combining NR and SEOBNRv2T-waveforms [513] in the frequency domain and fits directly the frequency domain phase;

  2. (ii)

    Kawaguchi et al. use a model function of the form

    $$\begin{aligned} \psi _T= & {} \frac{3 x^{5/2}}{128 \nu } \left( -\frac{39}{2} {\tilde{\varLambda }} (1+a \varLambda ^{2/3} x^p) \right) \nonumber \\&\left( 1+ {\tilde{c}}_1 x + {\tilde{c}}_{3/2} x^{3/2} + {\tilde{c}}_2 x^2 + {\tilde{c}}_{5/2} x^{5/2} \right) \end{aligned}$$
    (91)

    to include non-linear tidal effects instead of higher order PN effects as in Eq. (84);

  3. (iii)

    The model is valid up to a frequency of \(1000\,{\mathrm{Hz}}\) to ensure that no postmerger frequency contributions might effect the inspiral approximant.

In addition, the model includes a tidal amplitude correction of the form

$$\begin{aligned} \begin{aligned} {\tilde{A}}_T^{{\mathrm{Kawaguchi}}}&= \sqrt{\frac{160 \pi \nu }{27}} \frac{M^2}{D_L} \ {\kappa ^T_{{\mathrm{eff}}}} x^{-7/4} \\&\quad \times \left( -\frac{27}{16}x^5-\frac{449}{64}x^6-4251 x^{7.890}\right) . \end{aligned} \end{aligned}$$
(92)

6.3 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, i.e., tidal contributions for PN models are smaller than for the phenomenological models using the same tidal deformability. 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 \(1000\,{\mathrm{Hz}}\) 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 \({\tilde{\varLambda }} = 392.1\), systematic waveform uncertainties seem to be negligible; see also [404]. 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
figure 11

Comparison of the tidal phase of different PN and phenomenological waveform approximant as discussed in the main text; cf. [182] for further details

7 Analyzing BNS signals

As emphasized in Sect. 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.

7.1 Estimating source properties

The first analysis to establish the possibility of using a Bayesian approach to estimate BNS source properties was presented in [170] 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. Reference [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 [460]. 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 [537, 538]. The information from individual parameters characterizing the source is encoded in posterior probability distributions

$$\begin{aligned} p(\varvec{\theta }|{\mathcal {H}}_s,d) = \frac{p(\varvec{\theta }|{\mathcal {H}}_s)p(d|\varvec{\theta }, {\mathcal {H}}_s)}{p(d|{\mathcal {H}}_s)}. \end{aligned}$$
(93)

Here, \(\varvec{\theta }\) represents the parameter set and \({\mathcal {H}}_s\) is the hypothesis that a GW signal depending on the parameters \(\varvec{\theta }\) 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, \chi _A, \chi _B, \theta , \phi , \iota , \psi , D_L, t_c, \varphi _c\}\), and in addition the tidal deformability components \(\varLambda _A\) and \(\varLambda _B\). The likelihood of obtaining a signal h(t) in data stream \(d(t) = h(t) + n(t)\), is proportional to

$$\begin{aligned} p(d|\varvec{\theta },{\mathcal {H}}_s) \propto \exp {\left[ -\frac{1}{2}(d-h|d-h)\right] }. \end{aligned}$$
(94)

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

$$\begin{aligned} p(d|{\mathcal {H}}_{{\mathrm{EOS}}_1}, I) = \int p(d|{\mathcal {H}}_{{\mathrm{EOS}}_1},\varvec{\theta },I) p(\varvec{\theta }|{\mathcal {H}}_{{\mathrm{EOS}}_1},I) d\varvec{\theta } \end{aligned}$$
(95)

By computing \(p(d | {\mathcal {H}}_{{\mathrm{EOS}}_2}, I)\) for another EOS, the ratio between their evidences

$$\begin{aligned} B^1_2 = \frac{p(d|{\mathcal {H}}_{{\mathrm{EOS}}_1}, I)}{p(d | {\mathcal {H}}_{{\mathrm{EOS}}_2}, I)} \end{aligned}$$
(96)

Ref. [170] looked into combining multiple sources to distinguish between the EOSs they follow and concluded that \(\sim 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 post-merger 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. [95, 97, 122, 529]. Besides extracting the source parameters, future detections might also allow us to extract EOS-insensitive quasi-universal relations from GW data [479].

7.2 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 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 point-mass 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 \(\varLambda _i\) between [0,5000] and estimated \({\tilde{\varLambda }}\) by effectively using a prior uniform on \({\tilde{\varLambda }}\). 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 \({\tilde{\varLambda }}\) were found to be stronger for the narrower spin priors. With wider spin priors, the component masses were estimated to lie between 1.00 \(\mathrm {M}_\odot \) and 1.89 \(\mathrm {M}_\odot \). With narrower priors on spins, the component masses were found to lie between 1.16 \(\mathrm {M}_\odot \) and 1.60 \(\mathrm {M}_\odot \). The estimate on the upper bound of \({\tilde{\varLambda }}\) was 630 for wider spin priors, whereas within 90% credible interval, an estimate of \(300^{+420}_{-230}\) was made. The posteriors on \({\tilde{\varLambda }}\) from the high spin priors using the four different waveform models mentioned above is shown in Fig. 12.

Fig. 12
figure 12

Posterior probability distributions on \({\tilde{\varLambda }}\) obtained from the waveforms TaylorF2 (blue, dotted line), SEOBNRv4_ROM_NRTidal (magenta, dashed line), IMRPhenomD_NRTidal (green, dash-dotted line), and IMRPhenomPv2_NRTidal (red, solid line) by analyzing the GW170817 signal. The data has been taken from the data release of [12]. The plot shows the results from imposing the high-spin priors

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(\rho )\) was used. Radii estimates of \(R_A = 10.8^{+2.0}_{-1.7}\,\hbox {km}\) and \(R_B=10.7^{+2.1}_{-1.5}\,\hbox {km}\) at 90% credible intervals were obtained for the primary and secondary components by imposing the EOS-insensitive relations. Further constraining bounds of \(R_A=11.9^{+1.4}_{-1.4}\,\hbox {km}\) and \(R_B=11.9^{+1.4}_{-1.4}\,\hbox {km}\) were obtained from the method of parametrization of the EOS, in addition to imposing the maximum mass of a NS to be \(> 1.97 \mathrm {M}_\odot \). Reference [111] 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 [168] 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.

Reference [290] 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 \({\tilde{\varLambda }}\) 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 [417] or that one or two components have been BHs can not be ruled out, e.g., [253, 310]. 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 [135, 221, 391]

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]. Reference [477] 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 \(\sim 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 unequal-mass 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 non-tidal 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. Figure 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. [478] 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 \(\chi > 0.2\). In addition, [478] 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^\circ \) 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.

Fig. 13
figure 13

Posterior probability distributions on \({\tilde{\varLambda }}\) when differing BBH baselines, the tidal description is always NRTidal. The injected value is shown by the vertical black line and the 90% credible intervals are marked for each model. The injected sources are equal-mass (left panel) with \(m_1=m_2=1.375\)\({M}_{\odot }\) and unequal-mass (right panel) with \(m_1 = 1.68\)\({M}_{\odot }\) and \(m_2 = 1.13\)\({M}_{\odot }\). The injected waveform model for both sources is with IMRPhenomPv2_NRTidal

While the above works focused on tidal estimates from the inspiral only, Ref. [190] 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 [380], 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.

Fig. 14
figure 14

Posterior probability distributions on \({\tilde{\varLambda }}\) for different tidal descriptions keeping the same baseline, IMRPhenomPv2. The injected value is shown by the vertical black line and the 90% credible intervals are marked for each waveform with dashed lines. The injected sources are equal-mass (left panel) with \(m_1=m_2=1.375\)\({M}_{\odot }\) and unequal-mass (right panel) with \(m_1 = 1.68\)\({M}_{\odot }\) and \(m_2 = 1.13\)\({M}_{\odot }\). The injected waveform model for both sources is with IMRPhenomPv2_NRTidal

7.4 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 \sim 5\,\hbox {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]. However, the power spectrum of the detectors over the duration of each BNS signal may no longer be stationary, and, in addition, the rotation of the Earth during the period, while a BNS signal lasts in the detectors’ bandwidths, may also bring in further subtleties in GW data analysis. Reference [463] 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.

8 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 review key aspects in the field of NR (Sect. 3), PN theory (Sect. 4), the EOB formalism (Sect. 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.