Analyzing Gravitational Waves with General Relativity

TL;DR

This paper reviews the theoretical methods used in gravitational wave analysis, emphasizing the role of approximation techniques like post-Newtonian methods in modeling signals for detectors such as LIGO, Virgo, and LISA.

Contribution

It provides a comprehensive overview of analytic approximation methods in general relativity, highlighting their application in gravitational wave detection and modeling compact binary systems.

Findings

Post-Newtonian approximation accurately models early inspiral phases.

Dimensional regularization solves divergence problems in calculations.

Finite size effects influence gravitational wave emission from neutron stars.

Abstract

After a short review of prominent properties of gravitational waves and the newly born gravitational astronomy, we focus on theoretical aspects. Analytic approximation methods in general relativity have played a crucial role in the recent discoveries of gravitational waves. They are used to build theoretical template banks for searching and analyzing the signals in the ground-based detectors LIGO and Virgo, and, further ahead, space-based LISA-like detectors. In particular, the post-Newtonian approximation describes with high accuracy the early inspiral of compact binary systems, made of black holes or neutron stars. It mainly consists of extending the Einstein quadrupole formula by a series of relativistic corrections up to high order. The compact objects are modelled by point masses with spins. The practical calculations face difficult problems of divergences, which have been solved thanks to the dimensional regularization. In the last rotations before the merger, the finite size effects and the internal structure of neutron stars (notably the internal equation of state) affect the evolution of the orbit and the emission of gravitational waves. We describe these effects within a simple Newtonian model.

Equations

$$h={\mathcal{F}}_{+}{h}_{+}+{\mathcal{F}}_{\times }{h}_{\times },$$

$$\left(\begin{array}[]{l}h_{+}\\[5.69046pt] h_{\times}\end{array}\right)=\frac{2G}{c^{4}d}\left(\begin{array}[]{l}\frac{p^{i}p^{j}-q^{i}q^{j}}{2}\\[5.69046pt] \frac{p^{i}q^{j}+p^{j}q^{i}}{2}\end{array}\right)\biggl{\{}\frac{\mathrm{d}^{2}Q_{ij}}{\mathrm{d}t^{2}}\bigl{(}t-d/c\bigr{)}+\mathcal{O}\left(\epsilon_{\text{PN}}\right)\biggr{\}}+\mathcal{O}\left(\frac{1}{d^{2}}\right)\,.$$

$$Q_{ij}=\int\mathrm{d}^{3}\mathbf{x}\,\rho\,\Bigl{(}x^{i}x^{j}-\frac{1}{3}\mathbf{x}^{2}\Bigr{)}\,.$$

$$\frac{\mathrm{d}E}{\mathrm{d}t}=-F^{\text{GW}},$$

$$F^{\text{GW}}\equiv \frac{G}{5c^5}\left\{\frac{{\mathrm{d}}^3{Q}_{ij}}{\mathrm{d}{t}^3}\frac{{\mathrm{d}}^3{Q}_{ij}}{\mathrm{d}{t}^3}+\mathcal{O}\left({ϵ}_{\text{PN}}^2\right)\right\}.$$

$$\frac{\mathrm{d}{J}_i}{\mathrm{d}t}=-\frac{2G}{5c^5}{\epsilon }_{ijk}\left\{\frac{{\mathrm{d}}^2{Q}_{jl}}{\mathrm{d}{t}^2}\frac{{\mathrm{d}}^3{Q}_{kl}}{\mathrm{d}{t}^3}+\mathcal{O}\left({ϵ}_{\text{PN}}^2\right)\right\}.$$

$$\frac{\mathrm{d}{P}_i}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{G}_i}{\mathrm{d}t}$$

$$\left(\begin{array}[]{l}h_{+}\\[5.69046pt] h_{\times}\end{array}\right)=\frac{2Gm\nu}{c^{2}d}\left(\frac{Gm\omega}{c^{3}}\right)^{2/3}\left(\begin{array}[]{l}\bigl{(}1+\cos^{2}\!i\bigr{)}\cos(2\phi)\\[5.69046pt] \bigl{(}2\cos i\bigr{)}\sin(2\phi)\end{array}\right)\,,$$

$$\frac{\dot{\omega }}{{\omega }^2}=\frac{96\nu }{5}{\left(\frac{Gm\omega }{c^3}\right)}^{5/3}.$$

$$x=\left(\frac{Gm\omega}{c^{3}}\right)^{2/3}\,,\qquad\Theta=\frac{\nu c^{3}}{5Gm}\bigl{(}t_{\text{c}}-t\bigr{)}\,.$$

$$x$$

$$\varphi$$

$$E_{\text{ADM}}=m{c}^2-\frac{G{m}^2\nu }{2r}+{\int }_{-\infty }^t\mathrm{d}{t}^{\prime }{F}^{\text{GW}}(t^{\prime }),$$

$$E_{\text{ADM}}=M_{\text{c}}{c}^2+{\int }_{-\infty }^{t_{\text{c}}}\mathrm{d}t{F}^{\text{GW}}(t),$$

$$\Delta E^{\text{GW}}=(m-M_{\text{c}})c^2={\int }_{-\infty }^{t_{\text{c}}}\mathrm{d}t{F}^{\text{GW}}(t)=\frac{G{m}^2\nu }{2r_{\text{c}}}.$$

$$\phi=\phi_{0}-\frac{x^{-5/2}}{32\nu}\sum_{p}\Bigl{(}\varphi_{\text{$p$PN}}+\varphi_{\text{$p$PN}}^{(l)}\ln x\Bigr{)}\,x^{p}\,,$$

$${\phi }_{\text{−1PN}}$$

$${\phi }_{\text{0PN}}$$

$${\phi }_{\text{0.5PN}}$$

$${\phi }_{\text{1PN}}$$

$${\phi }_{\text{1.5PN}}$$

$${\phi }_{\text{2PN}}$$

$${\phi }_{\text{2.5PN}}^{(l)}$$

$${\phi }_{\text{3PN}}$$

$$+\left(-\frac{15737765635}{12192768}+\frac{2255}{48}{\pi }^2\right)\nu +\frac{76055}{6912}{\nu }^2-\frac{127825}{5184}{\nu }^3,$$

$${\phi }_{\text{3PN}}^{(l)}$$

$${\phi }_{\text{3.5PN}}$$

$${\phi }_{\text{4PN}}$$

$${\phi }_{\text{4PN}}^{(l)}$$

$${\phi }_{\text{4.5PN}}$$

$$+\left[\frac{1492917260735}{1072963584}-\frac{2255}{48}{\pi }^2\right]\nu -\frac{45293335}{1016064}{\nu }^2-\frac{10323755}{1596672}{\nu }^3)\pi ,$$

$${\phi }_{\text{4.5PN}}^{(l)}$$

$${\phi }_{\text{1.5PN}}^{\text{SO}}$$

$${\phi }_{\text{2.5PN}}^{\text{SO}(l)}$$

$${\phi }_{\text{3PN}}^{\text{SO}}$$

$${\phi }_{\text{3.5PN}}^{\text{SO}}$$

$$+\left(-\frac{170978035}{387072}+\frac{2876425}{5376}\nu +\frac{4735}{1152}{\nu }^2\right)\frac{\delta m}{m}{\Sigma }_z],$$

$${\phi }_{\text{4PN}}^{\text{SO}}$$

$$m_a\frac{\mathrm{d}{v}_a^i}{\mathrm{d}t}=G\sum _{b\ne a}\left[m_a{m}_b\frac{\partial }{\partial y_a^i}\left(\frac{1}{r_{ab}}\right)+\frac{1}{2}\left(m_a{q}_b^{jk}+m_b{q}_a^{jk}\right)\frac{{\partial }^3}{\partial y_a^i\partial y_a^j\partial y_a^k}\left(\frac{1}{r_{ab}}\right)\right],$$

Full text

Français: Après une brève revue des propriétés importantes des ondes gravitationnelles et de la nouvelle astronomie gravitationnelle, nous nous concentrons sur les aspects théoriques. Les méthodes d’approximation analytiques en relativité générale ont joué un rôle crucial dans les récentes détections d’ondes gravitationnelles. Elles sont utilisées pour créer des banques de modèles (patrons) théoriques qui servent à rechercher et analyser les signaux dans les détecteurs au sol LIGO et Virgo et, plus tard, les détecteurs dans l’espace de type LISA. En particulier, l’approximation post-newtonienne décrit avec grande précision le spiralement initial des systèmes binaires compacts de trous noirs ou d’étoiles à neutrons. Elle consiste principalement à étendre la formule du quadrupôle d’Einstein par une série de corrections relativistes jusqu’à un ordre élevé. Les objets compacts sont modélisés par des masses ponctuelles avec spins. Les calculs pratiques font face à des problèmes difficiles de divergences, qui ont été résolus grâce à la régularisation dimensionnelle. Dans les dernières orbites proches de la fusion, les effets de taille finie et de structure interne des étoiles à neutrons (notamment l’équation d’état interne) affectent l’évolution de l’orbite et l’émission des ondes gravitationnelles. Nous décrivons ces effets dans le cadre d’un modèle newtonien simple.

keywords: Gravitational waves; compact binary systems; post-Newtonian theory.

I Gravitational waves and the new astronomy

Paramount breakthroughs in Astronomy and fundamental Physics occurred with the discovery of gravitational waves (GWs) generated by the orbital motion and merger of compact binary systems, made of black holes or neutron stars Abbott et al. (2016, 2017). A new window of observations of our Universe opened up, radically different and complementary from that of the traditional astronomy, essentially based on electromagnetic (EM) waves. The salient properties of GWs shape the key features of the new “Gravitational Astronomy”:

• •

GWs are produced by the overall, “bulk” motion of large masses at relativistic speeds (close to the speed of light $c$), in contrast to EM waves which are in general composed of the incoherent superposition of photons emitted by the atoms and molecules composing the source Thorne (1987). As a result the wavelength of GWs is in general much larger than the size of the source, and there is a deep analogy between GWs and ordinary sound waves. However, in contrast to sound waves, GWs propagate in vacuum. They are ripples in the Riemannian curvature of space-time, which is the fundamental dynamical entity in general relativity (GR).

• •

GWs propagate almost without alteration through the densest regions of the Universe, and have thus the potential of carrying information from very far away — probably up to the first instants after the Big Bang. This is due to the weakness of the gravitational interaction as compared to other forces, and to the fact that GWs cannot be screened by any type of matter field. Indeed the charge associated with the gravitational interaction is the mass, which is always equal to the inertial mass by the equivalence principle, and therefore is always positive. Actually the positivity of the mass-energy of an arbitrary system (involving ordinary bodies and black holes), constitutes an important and difficult theorem in GR Schoen and Yau (1979); Witten (1981).

• •

GWs emitted by coalescing compact binary systems contain the information about their distance Schutz (1986). In this respect these systems constitute an analogue of the standard candels of EM-based astronomy (like Cepheid variables and Type Ia supernovas), and can rightly be called “standard sirens”. However in the case of the GW sirens there is no need for calibrating the distance scale; the calibration is automatically done by GR. One can thus measure the Hubble-Lemaître cosmological parameter $H_{0}$ with GWs, independently of the traditional EM measurements. 222A resolution of the XXX${}^{\text{th}}$ general assembly of the international astronomical union (IAU) recommended that the expansion of the Universe be referred to as the “Hubble-Lemaître law”.

• •

The gravitational astronomy is one of “precision”. It is possible to measure with high precision the parameters of compact binary systems (masses and spins) by direct comparison with a solution of the purely gravitational two-body problem in GR. See Fig. 1 which shows the signal of the first binary black hole event, directly analyzed with the GR prediction. So far no deviation from GR has been observed. For compact sources most of the non gravitational effects, which usually plague the dynamics of ordinary systems (magnetic fields, presence of an interstellar medium, etc.), are in general dominated by the gravitational force. 333We shall discuss an exception in Sec. IV: The internal structure and the non-gravitational equation of state of neutron stars do affect the GW signal close to and during the final merger.

• •

This highlights the crucial role played by analytic approximation methods (reviewed in Secs. II and III) and also numerical calculations, since they permit an accurate description of the two-body problem in GR, without which the full information contained into the signal could not be extracted. The theoretical solution of the problem of motion and radiation is used in the form of accurate GW templates, which are correlated to the observed signal using the technique of matched filtering Buonanno and Sathyaprakash (2015).

• •

The new astronomy is also “fundamental”. As far as we know, our gravitational theory is fundamental, and GR may be valid in a large range of energies, perhaps up to the Planck scale. Thus GW observations have a lot of implications for fundamental physics. With GWs one can confront GR with alternative theories of gravity (such as scalar-tensor theory, massive gravity theory, etc.) and one can test fundamental principles such as the equivalence principle. One can also question the standard model of cosmology $\Lambda$-CDM, including the great mysteries of contemporary physics constituted by cold dark matter (CDM), the cosmological constant $\Lambda$ and dark energy Barack et al. (2018).

• •

Last but not least: The multi-messenger aspect of the gravitational astronomy, i.e., its synergy with other vectors of information, most importantly EM waves but maybe also neutrinos in the future. For instance, the joint observation of GWs from a binary neutron star event and of a gamma ray burst (GRB) permitted to show that the speed of GWs is equal to $c$ with a precision $\sim 10^{-15}$. This ruled out a series of alternative theories of gravity. But of course, the multi-messenger astronomy has outstanding implications in astrophysics, such as refining the model of GRBs, understanding the explosions of kilonovas and the mechanisms for the production of heavy elements.

Consider an isolated (with finite spatial support) source of GWs. Let $m$ be the mass of the source, $r$ its size and $\omega\sim 2\pi/P$ the typical angular frequency of oscillations of the source, with $P$ the typical period. We have $\omega\sim v/r$ with $v$ the typical internal velocity. We suppose that the source is self gravitating, so the gravitational force is responsible for the dynamics and the GW generation. The wavelength of the emitted GW is $\lambda\sim cP/2$, with a factor $1/2$ inserted to take into account the quadrupolar nature of GWs, whose frequency is $f\sim\omega/\pi$. Note that $r/\lambdabar\sim 2v/c$ (posing $\lambda=2\pi\lambdabar$), so that the source is much smaller than one wavelength of the emitted GW, since $\epsilon_{\text{PN}}\sim v/c\ll 1$ in the non relativistic (post-Newtonian) regime. For the self gravitating source we typically have $Gm\sim r^{3}\omega^{2}$, indeed this is exactly Kepler’s third law $Gm=a^{3}\omega^{2}$ in the case of a Newtonian binary system (with $a$ the semi-major axis of the orbit). Furthermore, if the source is compact its size is of the order of the Schwarzschild radius, $r\sim 2Gm/c^{2}$, and the GW frequency scales inversely proportional to the mass: $Gmf\sim c^{3}/(\pi\sqrt{8})$. But the mass decreases because the GW extracts energy from the source, hence the frequency of the GW increases: This is the famous “chirp” of GWs, which we shall compute with high precision in Sec. III.

II Successes with the Einstein quadrupole formula

The “precision” gravitational astronomy requires inputs from the theory side. Searching and analyzing GW signals that are well predicted by GR is made using the technique of matched filtering, which cross correlates the detector output with our best prediction of the expected signal, called the template. The template is weighted (in the Fourier domain) by the power spectral density of the noise in the detector. It depends on a set of trial parameters describing the source’s model (such as masses and spins), and that are measured in the process. As GR is a complicated non-linear theory, there is no hope of finding an exact solution of the Einstein field equations, but tremendous progresses have been made with the development of perturbative and approximation methods in GR, notably the post-Newtonian (PN) approximation. Conjointly with analytic developments, continued efforts in numerical relativity led to the computation of the final merger of binary black holes and other GW sources like supernova explosions.

The first analytic computation of GWs is the Einstein quadrupole formula Einstein (1916, 1918), valid at the dominant “Newtonian” order in a PN expansion, with the small PN parameter being the slowness estimate $\epsilon_{\text{PN}}\sim v/c$, ratio of a typical internal velocity in the source and the speed of light. Originally derived for matter sources with negligible self gravity (hence the source’s oscillations producing GWs have a non gravitational origin), the formula was later shown to be still valid for weakly self-gravitating sources, such as a Newtonian binary system Landau and Lifshitz (1971). The GW amplitude is characterized by two tensorial polarization modes, traditionally denoted $h_{+}$ and $h_{\times}$, that are transverse to the direction of propagation $\bm{n}=(n^{i})$ (with $i=1,2,3$ and $\bm{n}^{2}=1$), pointing from the GW source towards a far-away detector. The detector is sensitive to a certain linear combination of the two polarizations,

[TABLE]

where the “form factors” $\mathcal{F}_{+}$ and $\mathcal{F}_{\times}$ depend on the direction and orientation of the source with respect to the local frame of the detector. The polarizations are defined as the projection of the waveform along two polarization vectors $\bm{p}$ and $\bm{q}$ in the plane orthogonal to $\bm{n}$, and forming with it an orthonormal right-handed triad $(\bm{n},\bm{p},\bm{q})$. The quadrupole formula gives the polarizations at large distance $d$ from the source (and at retarded time $t-d/c$) as

[TABLE]

The quadrupole moment $Q_{ij}$ of the source is just, at the leading approximation, the usual mass type moment of the Newtonian mass density $\rho$ in the source,

[TABLE]

The total energy $E$ of the matter source decreases because of the GW emission, and this is controlled by the “flux-balance” equation

[TABLE]

where the GW flux in the right-hand side is given at the leading PN approximation (i.e., “Newtonian” order in the radiation field) by

[TABLE]

Witness the factor $c^{-5}$ in front of the flux, which shows that the corresponding radiation reaction effect in the matter equations of motion is actually of order 2.5PN, namely $\mathcal{O}(\epsilon_{\text{PN}}^{5})$.444By order $n$PN we refer to a small post-Newtonian term of the order of $\epsilon_{\text{PN}}^{2n}$. Similarly there is a quadrupole flux-balance equation for the angular momentum $J_{i}$, given by

[TABLE]

The laws of motion of a relativistic conservative system (neglecting the GW emission) admit ten Noetherian invariants associated with the symmetries of the Poincaré group. In addition to the energy $E$ and angular momentum $J_{i}$, there is the linear momentum $P_{i}$ and the invariant of the center of mass $G_{i}$. The latter CM invariant is associated with the invariance under Lorentz boosts. When the GWs are turned on, all the invariants obey some flux-balance equations. In addition to (4)–(6) we have

[TABLE]

where $Q_{ijk}$ is the Newtonian mass octupole moment and $D_{ij}$ is the current type quadrupole moment. Notice that (7) represent subdominant radiation reaction effects of order 3.5PN. The flux of linear momentum $P_{i}$ is well known, as it is responsible for the “gravitational recoil” of the source by GW emission (see, e.g., Blanchet et al. (2005a)). However, strangely enough, the expression of the flux of center-of-mass position $G_{i}$ in Eq. (7b), has only been computed and recognized recently Kozameh et al. (2018); Nichols (2018); Blanchet and Faye (2018).

The equations (4)–(6) for energy and angular momentum give the evolution of the orbital parameters (semi-major axis and eccentricity) of the compact binary system under GW emission Peters and Mathews (1963); Peters (1964). An average over the orbital period is applied, so as to consider the secular evolution of the orbit on a radiation reaction time scale much longer than the orbital period. The first success of the quadrupole formula has been that it works perfectly when accounting for the observed decay of the orbital period of the Hulse-Taylor binary pulsar Hulse and Taylor (1975). This test represented the first quantitatively precise proof of the existence of GWs Taylor et al. (1979); Taylor and Weisberg (1982); Damour and Taylor (1991). Nevertheless, since $\epsilon_{\text{PN}}\sim 10^{-3}$ is very small for binary pulsars, the quadrupole formula is not expected to yield any deviation with respect to observations in the regime of binary pulsars.555But see the controversial debate on this point at the time of the binary pulsar Ehlers et al. (1976); Walker and Will (1980).

Even more impressive, a second success occurred recently, because the GW signals from binary black hole and neutron star events can be reasonably well interpreted with the quadrupole formula. In the case of black hole binaries this is truly remarkable because $\epsilon_{\text{PN}}\sim 0.5$ in the last rotations. Take the example of the first event GW150914, shown in Fig. 1. We are observing the signal at high frequency, close to the final merger, so the orbit has been circularized by radiation reaction — a consequence of the balance equation for angular momentum (6). For a binary system modeled by two point masses, and moving on a circular orbit, the GW polarizations (2) become

[TABLE]

where $\nu=m_{1}m_{2}/(m_{1}+m_{2})^{2}$ denotes the symmetric mass ratio between the two compact objects, $m=m_{1}+m_{2}$ is the total mass, and $i$ is the inclination angle of the binary’s orbital plane with respect to the plane of the sky (the polarization vector $\bm{p}$ pointing by convention towards the “ascending node”). In the quadrupole approximation the phase of the signal is $\phi^{\text{GW}}=2\phi$, where $\phi=\int\omega\,\mathrm{d}t$ is the orbital phase and $\omega$ the angular frequency; the signal frequency is usually denoted $f^{\text{GW}}=\omega/\pi$.

For the purposes of detection and subsequent data analysis, the most important information provided by GR is the time evolution of the phase and frequency, that are computed from the energy balance prescription (4). For circular orbits neither the angular momentum balance equation (6) nor the averaging procedure are necessary. Both $E$ and the flux $F^{\text{GW}}$ are only functions of the orbital frequency $\omega$, and the flux is readily computed from Eq. (5) in the case of a Newtonian system of two point masses on a circular orbit. Hence the balance equation becomes an ordinary differential equation for the frequency:

[TABLE]

By integrating this equation one successively obtains the orbital frequency as a function of time, and the orbital phase as a function of frequency. For convenience we use the dimensionless variables

[TABLE]

Note that $x$ can be seen as a small PN parameter of the order of $\mathcal{O}(\epsilon_{\text{PN}}^{2})$; and $t_{\text{c}}$ denotes the instant of coalescence, at which the distance between the particles formally vanishes and the frequency diverges. With those notations, and with $\phi_{0}$ denoting an initial constant phase, we find

[TABLE]

These formulas, together with (8), describe the “chirp” of GWs at the lowest approximation, i.e., the way the frequency, phase and amplitude of the signal increase, untill some point at the onset of the merger of the compact objects, at which the approximation is no longer valid.

Inspection of Eqs. (8)–(11) shows that the GW signal depends on one combination of the two masses only, called the chirp mass and given by $\mathcal{M}=m\,\nu^{3/5}$. Roughly speaking there are two observables, the amplitude $h$ and the frequency chirp $\dot{\omega}$, from which one can determine at once the chirp mass $\mathcal{M}$ and the distance $d$. Actually things are more complicated because in a given detector, the measured amplitude is a certain linear combination (1) of the two polarizations, which also depends on the direction of the source and its orientation with respect to the detector. We need thus several detectors to determine these extra angles, and also the inclination angle $i$ in Eq. (8). Note that the measured mass is the redshifted one, $\mathcal{M}=(1+z)\mathcal{M}_{\text{source}}$, where $z$ is the cosmological redshift of the source ($z\sim 0.1$ in the case of GW150914) and $\mathcal{M}_{\text{source}}$ is the actual chirp mass of the binary. As for the measured distance, it is exactly the “luminosity distance” $d\equiv d_{\text{L}}$ used by cosmologists, who refer to the chirping binaries as standard sirens. Finally, the quadrupole formula gives consistent estimates (even for GW150914 !) for the mass, the distance, the maximal amplitude of the signal, and the number of orbital cycles from the entry frequency of the detector’s band till the merger, with the proviso that we suplement the quadrupole formula with an information from full GR, namely that the merger occurs at a separation of the order of the mass. Then, of course, the merger itself can only be described numerically.

In GR there is a notion of the total energy contained in the space-time, including both the contribution from the matter sources and that of the gravitational field: This is the so-called Arnowitt-Deser-Misner (ADM) mass-energy Arnowitt et al. (1962), which is exactly conserved. Let us apply this notion to the problem of the coalescence of two black holes, assuming a very crude model, in which the binary’s orbit is merely Newtonian and circular, and the radiation field is described by the quadrupole formula (4). At very early times $t\to-\infty$, the black holes were moving almost freely on quasi-hyperbolic orbits, and later formed a gravitationally bound system by GW emission, which then spiralled in till the merger. At any time the ADM energy is

[TABLE]

where the quadrupole flux is given by (5) — this equation is an integrated version of Eq. (4), which represented the time-variation of the Bondi mass Bondi et al. (1962). Since initially the black holes were infinitely far apart ($r\to\infty$), the conserved ADM energy is just equal to $mc^{2}$, while after the merger it also equals

[TABLE]

where $M_{\text{c}}$ is the mass of the black hole formed by the coalescence, occuring at the instant $t_{\text{c}}$. Hence we see, using the constancy of the ADM mass, that the energy which has been radiated away during the process is

[TABLE]

That is, it is equal to the mechanical binding energy of the binary system at the instant of coalescence $t_{\text{c}}$, which occurs when the separation is of the order of the mass, say $r_{\text{c}}\sim\frac{2Gm}{c^{2}}$. In the case of GW150914 the final mass $M_{\text{c}}$ has been measured, and about three solar masses have been released in the form of GWs in a few tens of second. This corresponds to a total power of the order of $10^{49}\,\text{W}$, well consistent with the result of sophisticated numerical simulations of the merger of two black holes in GR. This power is huge, but notice that it is only about a thousandth of $c^{5}/G$, which represents the natural general relativistic scale for a power, i.e., the Planck scale which happens in this case not to depend on Planck’s constant $\hbar$. Thus we may imagine that the gravitational astronomy could some day bring us a surprise, with the discovery of a GW source that is even more powerful than a binary black hole system.

III Analytic approximation methods for computing the GW chirp

Though the quadrupole formula is very useful, it is not sufficient when we want to perform precise calculations, and it becomes inoperational in the final merger phase of the compact objects. As we have seen, it does not permit to measure the two masses $m_{1}$ and $m_{2}$ separately, but only the chirp mass $\mathcal{M}$. The degeneracy over the masses is removed by including the relativistic PN corrections $\sim\mathcal{O}(\epsilon_{\text{PN}}^{2})$ in the quadrupole formula (4). In addition, the PN corrections depend on the spins of the compact objects (i.e., their intrinsic classical angular momenta). Taking into account the spins in the GR templates is important, and their measurement represents a valuable astrophysical information.

At the beginning of the building of the LIGO/Virgo detectors in the end of the 1980s, it was thought that the quadrupole formula is sufficient for detecting and analyzing binary neutron star coalescences. But in the early 1990s it was realized Cutler et al. (1993); Blanchet and Schäfer (1993) that neutron star binaries spend thousands of cycles in the band of the detectors, and that the GW templates should be able to monitor the signal with a precision of a fraction of a cycle over the entire bandwidth, say $\frac{\delta\phi}{2\pi}\lesssim 0.1$ where $\phi$ is the orbital phase (11b). It was then estimated that PN corrections in the quadrupole formula must be developed up to at least the daunting 3PN order $\sim\epsilon_{\text{PN}}^{6}$. Also recall that the quadrupole formula, as seen as a small radiation reaction contribution in the dynamics of the source, is itself a 2.5PN effect $\sim\epsilon_{\text{PN}}^{5}$ relatively to the Newtonian acceleration. This shows the highly relativistic character of compact binary systems observed in the LIGO/Virgo band as compared to the binary pulsar, for which the PN corrections in the orbital $\dot{P}$ are negligible. At the time of the binary pulsar only the 1PN correction to the quadrupole formula was known Epstein and Wagoner (1975); Wagoner and Will (1976); Blanchet and Damour (1989); Blanchet and Schäfer (1989). Furthermore, it was also realized that the PN corrections are important not only for the precise off-line analysis of the signals once they are detected, but also for the on-line process of detection Cutler et al. (1993).

The coalescence signal of two compact objects can be decomposed into three phases, see Fig. 2:

• •

The early inspiral during which the frequency and amplitude of the signal chirp with time (the chirp was discussed and defined in Sec. I). During this phase the signal is universal, i.e., it does not depend on the nature of the compact objects, be they black holes or neutron stars or more exotic objects like boson stars. The signal depends only on the masses and the spins. The PN approximation constitutes the ideal tool for describing the inspiral phase. For low mass compact binaries such as double neutron stars, the detectors are mostly sensitive to the inspiral phase, and the currently known analytical PN templates are accurate enough for detection and analysis, at least for moderate spins. Thus, the data analysis of neutron star binary events like GW170817 Abbott et al. (2017) is based on the PN templates.

• •

The merger phase, when the dynamics undergoes a transition from adiabatic inspiral to some unstable plunge, followed by the rapid collapse of the two objects to form a black hole. Numerical relativity (NR) has succeeded in the 2005s to compute the merger of two black holes Pretorius (2005); Baker et al. (2006); Campanelli et al. (2006). A non trivial point (which did not seem to be obvious some years ago Brady et al. (1998)), is that the overlap between the PN and NR regimes exists and is quite significant.666Because of prohibitive computing times the NR calculations are limited to a few tens of orbits before the merger, and will likely never be competitive with the PN approximation when monitoring tens of thousands of cycles in the early inspiral. The important issue of matching the PN and NR waveforms has been solved using several techniques Boyle et al. (2008); Pan et al. (2010). One is the hybrid inspiral-merger-ringdown (IMR or IMR-Phenom), which consists of introducing between the PN and NR domains of validity an overlapping time interval that is parametrized in a phenomenological way Santamaria et al. (2010); Ajith et al. (2011). The other technique recasts the actual PN two-body equations of motion and radiation into a simpler effective one-body (EOB) form Buonanno and Damour (1999); Damour and Nagar (2011). The EOB dynamics is described in a non perturbative way, which permits to extend the domain of validity of the PN approximation. A variant of EOB called EOB-NR is matched to the NR results. Both the IMR and EOB waveforms are extensively used in the LIGO/Virgo data analysis of the recent binary black hole events.

• •

The “ringdown”, when the newly formed black hole, which is highly deformed due to the nonlinear dynamics of the collision, relaxes to a stationary configuration given by the Kerr solution — the unique stationary rotating black hole in GR, depending only on the mass and the spin. The perturbed black hole emits quasi-normal mode radiation, and the modes can be analyzed by comparing with black hole perturbation theory Berti et al. (2006). A test of the “no-hair” theorem for black holes can be implemented by looking at the presence of an abnormal quadrupole moment endowed by the black hole, which would be independent from the mass and the spin. Alternatively this can be viewed as a test of the existence of a new, exotic form of matter alternative to black holes.

We now review the state-of-the-art on PN approximations applied to the inspiral of two compact objects. The first problem is that of the equations of motion, and has been solved up to the 4PN order $\sim\epsilon_{\text{PN}}^{8}$ for non-spinning compact bodies. 777See Blanchet (2014) for an history of the problems of motion and radiation, and for references to previous PN approximations. Different methods have been used, with equivalent results: The Hamiltonian formalism in ADM coordinates Jaranowski and Schäfer (2012, 2013); Damour et al. (2014), and the Fokker action of GR in harmonic coordinates Bernard et al. (2016); Bernard et al. (2017a, b); Marchand et al. (2018). In addition partial results have been obtained with the effective field theory Foffa and Sturani (2013); Galley et al. (2016). The second problem is the one of the GW field, and of course, this is that problem whose solution is directly used by LIGO/Virgo. Here the state-of-the-art is 3.5PN order $\sim\epsilon_{\text{PN}}^{7}$ beyond the result of the quadrupole formula Blanchet and Schäfer (1989, 1993); Blanchet (1995); Blanchet et al. (1995a); Blanchet (1998a); Blanchet et al. (2002a, b); Blanchet and Iyer (2004); Blanchet et al. (2004), and the 4.5PN term is also known Marchand et al. (2016). To reach this result a cocktail of approximation methods in GR called “MPM-PN” has been used:

• •

In a first stage we control the gravitational field generated by an isolated matter system in the exterior zone of the system. A non-linearity or post-Minkowskian (PM) expansion is combined with a multipolar (M) expansion parametrized by some sets of source multipole moments, yielding the most general solution of the Einstein field equation in the exterior zone Blanchet and Damour (1986); Blanchet (1987); Blanchet and Damour (1992); Blanchet et al. (1995b). In particular this solution recovers the Bondi-Sachs formalism Bondi et al. (1962); Sachs (1962) for the asymptotic structure of radiative fields at infinity from the matter source. This is the MPM part of the method.

• •

The MPM solution is matched to the PN field in the near and interior zones of the source. This is achieved by a matching equation, within a specific variant of the theory of matched asymptotic expansions. The matching is performed up to any PN order and yields unique expressions for the multipole moments of the source Blanchet (1995, 1998b) as well as for the radiation reaction contributions in the inner PN metric Poujade and Blanchet (2002); Blanchet et al. (2005b). This completes the MPM-PN approach.

• •

The MPM-PN solution is applied to systems of compact bodies treated as point masses (possibly with spins), using delta functions singularities. The model thus entails ultra-violet (UV) divergences, that are cured by means of dimensional regularization. Here we borrow dimensional regularization from quantum field theory and use it as a powerful regularisation scheme in classical GR. In addition, independently of the model of point particles, there are infra-red (IR) divergences in the general formalism, treated by a variant of dimensional regularization.

We report the most complete results concerning the PN corrections in the orbital phase — crucial for both processes of detection and subsequent parameter analysis. Extending the “Newtonian” result (11b) we have:

[TABLE]

where the sum runs over the successive PN corrections: $x^{p}\sim\mathcal{O}(\epsilon_{\text{PN}}^{2p})$, with $p$ being an integer or half integer. Some of the PN terms involve the logarithm of $x$, which we indicate by adding a superscript $(l)$ to the PN parameter. Up to the 3.5PN order we have Blanchet et al. (1995a); Blanchet (1998a); Blanchet et al. (2002a, b); Blanchet and Iyer (2004); Blanchet et al. (2004) 888Here $\gamma_{\text{E}}$ is the (probably irrational/transcendental) Euler’s constant.

[TABLE]

The PN parameters have been factorized by the dominant quadrupolar effect, so that $\varphi_{\text{0PN}}=1$ by definition. Notice the term \varphi_{\text{-1PN}} which corresponds to a dipolar effect; this term would appear in scalar-tensor theory but is absent in GR, hence \varphi_{\text{-1PN}}=0. The 1.5PN term is especially interesting as it corresponds to the dominant nonlinear tail effect — backscattering of linear quadrupolar GWs onto the space-time curvature generated by the mass of the source. In scalar-tensor theory there would be a dipolar tail term at the 0.5PN order Bernard (2018a, b), but again this effect is absent in GR. The observational limits on the measurement of the PN parameters by LIGO/Virgo are shown in Fig. 3.

Going beyond 3.5PN order, the 4PN parameter has not yet been computed by PN theory, but its leading term in the small mass ratio limit $\nu\to 0$ is known from black hole perturbation theory Tagoshi and Sasaki (1994); Tanaka et al. (1996):

[TABLE]

We expect that PN theory will be able to fully confirm this result, as well as of course, to provide all the mass ratio corrections $\mathcal{O}(\nu)$ therein. Finally the complete 4.5PN parameter has been derived by PN theory, and is due to an iterated nonlinear tail effect Marchand et al. (2016):

[TABLE]

We emphasized that in the matched filtering analysis of GWs, it is important to take into account the effects of spins. Here we report the spin-orbit (SO) coupling contributions, which are linear in the two spins and result from the coupling with the orbital angular momentum. They are known up to the 4PN order Bohé et al. (2013); Marsat et al. (2013):

[TABLE]

where the mass difference is denoted by $\delta m=m_{1}-m_{2}$, $\bm{S}_{1}$ and $\bm{S}_{2}$ are the two individual spins, and $S_{z}$ and $\Sigma_{z}$ are the projections of the particular combinations $\bm{S}=\bm{S}_{1}+\bm{S}_{2}$ and $\bm{\Sigma}=\frac{m}{m_{2}}\bm{S}_{2}-\frac{m}{m_{1}}\bm{S}_{1}$ perpendicular to the orbital plane, i.e., parallel to the orbital angular momentum. The quadratic spin-spin (SS) coupling contributions to the orbital phase are also known Bohé et al. (2015).

IV Influence of the internal structure of compact bodies

The PN parameters have been obtained within the so-called “pole-dipole” model, which approximates the rotating compact body as a point mass with a spin, but neglects the effect of the finite size and the internal structure of the body, such as the internal velocity field and the type of equation of state. In particular, the quadrupolar tidal deformation of the body is ignored. Neutron stars have a strong internal gravity so it is very difficult to deform them. We expect that they should be distorted by the gravitational field of the companion only in the last orbits before the merger Kochanek (1992); Bildsten and Cutler (1992). On the other hand, the numerical computation of the merger of two neutron stars shows that it is strongly dependent on the internal structure and the (unknown) equation of state Shibata and Taniguchi (2011); Faber and Rasio (2012). Therefore a legitimate question to ask, is whether and at which PN order the internal structure of extended compact objects influences the orbital phase evolution. Here we answer this question by means of a simple Newtonian model for the tidal interaction between extended bodies without spins during the inspiral phase, at the lowest quadrupolar level. 999See Refs. Flanagan and Hinderer (2008); Hinderer (2008); Damour and Nagar (2009); Binnington and Poisson (2009); Bini et al. (2012); Favata (2014) for entries in the literature. The Newtonian equations of motion of $N$ extended spinless bodies ($a,b=1,\dots,N$), to linear order in the quadrupole moments, are

[TABLE]

where $m_{a}$ are the masses, and we denote the position and velocity of the center of mass of the bodies by $\bm{y}_{a}(t)$ and $\bm{v}_{a}(t)=\mathrm{d}\bm{y}_{a}/\mathrm{d}t$, with the Euclidean separation between centers of mass being $r_{ab}=|\bm{y}_{a}-\bm{y}_{a}|$. The quadrupole moments of the bodies, supposed to be made of a perfect fluid, read

[TABLE]

with $\mathrsfs{V}_{a}$ the volume of the body, $\bm{z}_{a}=\mathbf{x}-\bm{y}_{a}(t)$ the distance between a generic point $\mathbf{x}$ inside the body and the center of mass, $\rho_{a}=\rho(\bm{y}_{a}+\bm{z}_{a},t)$ the Newtonian mass density of the body, $\rho(\mathbf{x},t)$ being the usual Eulerian density. The mass-centred condition reads

[TABLE]

The conserved energy of the $N$-body system is the sum of the internal (Newtonian) energies $e_{a}$ and of the orbital contributions, including the quadrupole effects:

[TABLE]

where we have introduced the tidal field acting on body $a$ and due to the other bodies $b\not=a$:

[TABLE]

Posing $\bm{w}_{a}=\mathrm{d}\bm{z}_{a}/\mathrm{d}t$ for the internal velocity field of body $a$, $\Pi_{a}=\Pi(\bm{y}_{a}+\bm{z}_{a},t)$ for the specific internal energy satisfying the thermodynamical relation $\mathrm{d}\Pi=-P\,\mathrm{d}(1/\rho)$ (with $P$ the pressure), and $u_{a}$ for the internal self-gravity given by the Poisson integral over the volume of the body, we have

[TABLE]

The coupling of the quadrupole moment $q_{a}^{ij}$ with the external tidal field $\mathcal{E}_{a}^{ij}$ of the other bodies implies a variation of the internal energy given by

[TABLE]

We consider the case where the quadrupole moment is induced by the tidal field of the other bodies. To linear order, we can introduce a coefficient $\lambda_{a}$ characterizing the deformability (or “polarizability”) of the body under the influence of the external field, such that

[TABLE]

The “response” coefficient $\lambda_{a}$ depends on the internal structure of the body, and is commonly given as $\lambda_{a}=\frac{2}{3G}\,k_{a}r_{a}^{5}$ in terms of the radius $r_{a}$ of the body and the mass-type quadrupolar Love number $k_{a}\equiv k_{a}^{(2)}$ (see for instance Hinderer (2008)). In fact it will be more convenient to characterize the internal structure of the body by the dimensionless parameter

[TABLE]

In the case of the induced quadrupole moments (27), the total energy of the system becomes

[TABLE]

Consider a compact binary system ($N=2$) moving of an exact circular orbit. From Eq. (27) we see that the two quadrupole moments face each other, and remain constant along the circular orbit. The equation of the relative motion reduces to $\mathrm{d}\bm{v}/\mathrm{d}t=-\omega^{2}\bm{x}$, where $\bm{x}=\bm{y}_{1}-\bm{y}_{2}$ and $\bm{v}=\mathrm{d}\bm{x}/\mathrm{d}t$ are the relative position and velocity (with $r\equiv r_{12}$). We find from (20) the orbital frequency

[TABLE]

We pose $X_{a}=m_{a}/m$ so that $X_{1}X_{2}=\nu$ is the symmetric mass ratio, denote $\gamma=\frac{Gm}{rc^{2}}$ and employ the notation (28). In turn the conserved energy (29) for circular orbits reduces to

[TABLE]

As the effect of the deformation of the bodies computed here is purely “Newtonian”, we see that the $c$’s we have introduced into our definitions naturally cancel out in Eqs. (30) and (31).

The above dynamics is conservative, i.e., we have neglected the dissipative radiation reaction effect on the orbit. This effect is taken into account when we impose the flux-balance equation (4). Again there is no need to impose the angular momentum balance equation (6) for circular orbits. The total quadrupole moment of the system is the sum of the orbital one and of the intrinsic moments of the bodies, given by (21):

[TABLE]

Plugging this into the flux formula, computing the time derivatives using the equations of motion including the contributions from the quadrupole moments, see Eq. (30), and keeping only the terms linear in these quadrupoles, yields the flux (still for exact circular orbits) as

[TABLE]

At this stage we reexpress the invariants $E$ and $F^{\text{GW}}$ in terms of the orbital frequency $\omega$ instead of the separation distance $r$ using Eq. (30). The interest of doing this in GR (e.g., when doing relativistic PN calculations), comes from the fact that the separation $r$ depends on the choice of the coordinate system, while the orbital frequency $\omega$ is invariantly defined in a large class of coordinate systems. Recalling the definition for the invariant dimensionless PN parameter $x$ in (10), we obtain

[TABLE]

At this stage we can already draw a firm conclusion: The effect of the internal structure of non-spinning bodies is proportional to $x^{5}$, and is thus comparable to a relativistic effect occuring at the 5PN order. Recall though that we computed this effect using merely Newton’s law of gravity. Of course the latter estimate is just formal, but we expect it to be physically correct in the case of compact bodies. But the numerical value of the coefficient involving the $\Lambda_{a}$’s is to be taken into account. For instance $\frac{Gm_{a}}{c^{2}r_{a}}\sim 0.15$ for neutron stars (hence $x^{5}\sim 8\,10^{-5}$ at the merger), and the numerical estimates of the Love numbers for neutron stars are of the order of one or say, a tenth Damour and Nagar (2009); Binnington and Poisson (2009). Therefore the deformability parameters (28) for compact bodies should be of the order of $\Lambda_{a}\sim 1000$, depending of course on the equation of state, as shown in Fig. 4.

As in Sec. II the phase and frequency evolution follow from (4), where both $E$ and $F^{\textit{GW}}$ have been computed for the conservative dynamics in Eqs. (34). This approximation is justified as we are interested in the secular, adiabatic evolution of the orbit over a radiation reaction time scale. We need thus to evaluate the secular variation of the energy $E$, which we immediately find from (34a) to be

[TABLE]

Combining this with the flux (34b), we get an ordinary differential equation for $x$. It turns out to depend on the following combination of the two deformability parameters:

[TABLE]

so normalized that in the case of two identical neutron stars (with the same mass, $X_{1}=X_{2}=\frac{1}{2}$, and the same equation of state) it reduces to $\tilde{\Lambda}=\Lambda_{1}=\Lambda_{2}$. We obtain the quadrupole finite size effect due to the internal structure on the frequency and phase evolution, extending the point-mass results given by (11), as

[TABLE]

These corrections should be added linearly to the purely gravitational PN corrections presented in Sec. III. For compact bodies the effect appears at the very small order 5PN, but for neutron stars becomes relatively large at the merger, of the order of a radian in the phase. Remarkably, it has been possible to put a bound on the tidal deformability of neutron stars with the recent binary neutron star event GW170817, and to infer a constraint on several possible equations of state for the nuclear matter inside neutron stars, see Fig. 4.