The physics of gravitational waves

Enrico Barausse
SISSA, Via Bonomea 265, 34136 Trieste, Italy and INFN Sezione di Trieste and
IFPU - Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy

These lecture notes collect the material that I have been using over the years for various short
courses on the physics of gravitational waves, first at the Institut d’Astrophysique de Paris (France),
and then at SISSA (Italy) and various summer/winter schools. The level should be appropriate for
PhD students in physics or for MSc students that have taken a first course in general relativity. I try
as much as possible to derive results from first principles and focus on the physics, rather than on
astrophysical applications. The reason is not only that the latter require a solid understanding of the
physics, but it also lies in the trove of data that are being uncovered by the LIGO-Virgo-KAGRA
collaboration after the first direct detection of gravitational waves [1]. Any attempt to summarize
such a rich and fast changing landscape and its evolving astrophysical interpretation is bound to
become obsolete before the ink hits the page.
arXiv:2303.11713v1 [gr-qc] 21 Mar 2023

Contents

I. Prerequisites 2

II. The propagation and generation of gravitational waves 3

A. Linear perturbations on flat space 3
B. Linear perturbations on curved space 5
C. Linear perturbations on flat space: a scalar-vector-tensor decomposition 5
D. Generation of gravitational waves: a first derivation of the quadrupole formula 9
E. Dimensional analysis 11

III. Post-Newtonian expansion 11

A. The motion of massive and masseless bodies 11
B. The Einstein equations 12
C. A more rigorous derivation of the quadrupole formula 14

IV. Local flatness and the equivalence principle 16

A. The local flatness theorem and Riemann normal coordinates 16
B. Fermi Normal Coordinates 18

V. The stress energy tensor of gravitational waves 20

A. The gravitational contribution to the mass of a compact star 23

VI. The inspiral and merger of binary systems of compact objects 24

A. Geodesics in Schwarzschild and Kerr 26
B. A qualitative description of the inspiral and merger 29

VII. The post-merger signal 33

A. Scalar perturbations of non-spinning black holes 33
B. Tensor perturbations of non-spinning black holes 34
C. Tensor perturbations of spinning black holes 36

VIII. The detection of gravitational waves 37

A. The response of a gravitational wave detector: the low frequency limit 38
B. A geometric interpretation of the polarizations 40
C. The response of a gravitational wave detector: the transfer function 41

IX. Gravitational wave data analysis 44

A. Gaussian noise and power spectral density 45
1. Detection in the presence of noise 47
B. The signal-to-noise ratio for inspiraling binaries 48
C. Parameter estimation 49
 2

Acknowledgements 51

References 51

I. PREREQUISITES

These notes assume familiarity with Einstein’s equations, which in units with G = c = 1 can be written
as
Gµν = 8πTµν , (1)
µν µν
with G the Einstein tensor and T the matter stress energy tensor. It may be useful to recall that
because of the Bianchi identity ∇µ Gµν = 0, the stress energy tensor satisfies the “conservation” equation
∇µ T µν = 0 on shell.
For a perfect fluid, in the (− + ++) signature that we will use throughout these notes, the stress energy
tensor takes the form
T µν = (ρ + p)uµ uν + pg µν , (2)
µ
where u is the 4-velocity of the fluid element, ρ the energy density and p the pressure. With this ansatz,
the conservation of the stress energy tensor implies
dρ
= −(p + ρ)∇µ uµ , (3)
dτ
with τ the fluid element’s proper time, and
γ µν ∂ν p
aµ = − , (4)
ρ+p
with aµ = uν ∇ν uµ the 4-acceleration and
γ µν = g µν + uµ uν (5)
µ
the projector on the hypersurface orthogonal to u . Let us recall that Eq. 3 simply encodes the con-
servation of energy, while Eq. 4 generalizes the Newtonian Euler equation. In particular, for p = 0 the
relativistic Euler equation reduces to the geodesic equation aµ = 0.
It is worth recalling that the stress energy tensor can be defined in terms of the functional derivative
of the matter action, i.e.
2 δSpp
T µν = √ . (6)
−g δgµν
For a point particle of mass m, the action is simply given by
Z
Spp = −m dτ , (7)

where the integral is along the trajectory. By varying this action with respect to the trajectory, one
obtains the geodesic equation aµ = 0, while the functional derivative with respect to the metric yields
m   uµ uν
µν
Tpp ~
= √ δ (3) ~x − X(t) , (8)
−g ut
~
with X(t) the trajectory. This stress energy tensor can be mapped into that of a perfect fluid with p = 0
(dust). The same clearly applies to a collection of point particles. We can therefore conclude that for
point particles, the geodesic equation is implied by the conservation of the stress energy tensor.
Finally, let us recall that the “divergence” or “focusing” of neighboring geodesics is described by the
geodesic deviation equation. More precisely, given a family of geodesics xµ (v, τ ) labeled by a parameter
v and the proper time τ , let us introduce the separation vector v µ = (∂xµ /∂v)δv joining two neighboring
geodesics with parameters v and v + δv. This vector then satisfies the geodesic deviation equation
D2 v µ µ
= Rαβγ uα uβ v γ , (9)
dτ 2
µ
with Rαβγ the Riemann tensor and D/dτ the covariant derivative along the four velocity.
R
Exercise 1: From the action of a point particle, S = −m dτ , derive the geodesics equations by
varying with respect to the trajectory, and the stress energy tensor by varying with respect to the metric.
 3

II. THE PROPAGATION AND GENERATION OF GRAVITATIONAL WAVES

A. Linear perturbations on flat space

Let us start by considering generic vacuum perturbations hµν of a flat background spacetime. The
perturbed spacetime’s metric at linear order is therefore given by

gµν = ηµν + hµν , g µν = η µν − hµν , (10)

where indices are meant to be raised and lowered with the Minkowski metric ηµν . We will now derive
the (linearized) Einstein equations for the metric perturbation hµν .
Let us first recall that the gauge group of GR is given by diffeomorphisms, i.e. coordinate transforma-
tions. Under a transformation x̃µ = x̃µ (x), the metric transforms as

∂xα ∂xβ
g̃µν (x̃(x)) = gαβ (x) . (11)
∂ x̃µ ∂ x̃ν
If the coordinate transformation is “infinitesimal”, x̃µ = xµ + ξ µ with ∂ν ξ µ  1, a Taylor expansion of
Eq. 11 implies that the metric perturbation transforms as

h̃µν = hµν − Lξ ηµν = hµν − ∂µ ξν − ∂ν ξµ . (12)

where Lξ is the Lie derivative along the vector field ξ.

Let us use this gauge freedom to impose the Lorenz gauge condition

∂µ h̄µν = 0 , (13)

where we have introduced the trace-reversed metric perturbation h̄µν = hµν − 21 hηµν , with h = ηµν hµν
the trace. Note that this condition is also known as de Donder gauge condition, or also as harmonic
gauge condition (as it can be easily proven that it is equivalent to 2xν = 0, where 2 = η µν ∂µ ∂ν is the
flat space d’Alembertian and the xν are scalar functions defining the coordinates). In this gauge, it is
straightforward to show that the linearized Ricci tensor is
1
Rµν = − 2hµν . (14)
2
The Einstein equations can be written in the two equivalent forms
 
1 1
Rµν − Rgµν = 8πTµν ⇔ Rµν = 8π Tµν − T gµν . (15)
2 2

From the second expression, the linearized equations can therefore be written as
 
1
2hµν = −16π Tµν − T gµν ⇔ 2h̄µν = −16πTµν , (16)
2

respectively in terms of hµν and h̄µν . In vacuum (Tµν = 0), both quantities satisfy the homogeneous
wave equation

2h̄µν = 2hµν = 0. (17)

From this, it can already be seen that the metric perturbations (i.e. gravitational waves) on flat space
travel at the speed of light.
Let us now further investigate how many independent components/propagating degrees of freedom the
metric perturbation has. In principle, a symmetric rank-2 tensor has 10 independent components. The
Lorenz gauge condition is a vector equation and removes 4 of them, and hence one would expect 6 degrees
of freedom. However, let us note that the Lorenz gauge condition does not fix completely the gauge. In
fact, let us consider a metric perturbation hµν that respects the Lorenz condition, and let us perform an
infinitesimal change of coordinates. According to Eq. 12, one has

h̃µν = hµν − ∂µ ξν − ∂ν ξµ ,
(18)
h̃ = h − 2∂µ ξ µ ,
 4

and thus
¯
h̃µν = h̄µν − ∂µ ξν − ∂ν ξµ + ηµν ∂α ξ α . (19)

We therefore see that

¯
∂µ h̃µν = ∂µ h̄µν − 2ξ ν . (20)

Clearly, if one starts with a perturbation in the Lorenz gauge, any gauge related to the original one by
harmonic generators (i.e. ones such that 2ξ µ = 01 ) still satisfies the Lorenz gauge condition. In other
words, the Lorenz gauge is defined up to a harmonic gauge generator.2
Let us now exploit this residual gauge freedom to simplify the (trace reversed) metric perturbation
h̄µν in vacuum. Since the latter satisfies 2h̄µν = 0, it can be decomposed, without loss of generality, in
planes waves:
α
h̄µν (x) = Aµν eikα x + c.c. , (21)

where c.c. denotes the complex conjugate, Aµν is a constant polarization tensor, and k µ is a null-vector
(kα k α = 0). Similarly, the (residual) gauge generator must satisfy 2ξ µ = 0, and thus
α
ξ µ (x) = B µ eikα x + c.c. , (22)

where again kα k α = 0 and B µ is a constant vector. Imposing now the condition 13, one finds that

k µ Aµν = 0, (23)

i.e. the wavevector must belong to the null space of the polarization tensor.
Considering now for simplicity a wave propagating along the z axis, i.e.
∂ ∂
k= + , (24)
∂t ∂z
Eq. 23 yields

Atν = Azν . (25)

Evaluating the transformation given by Eq. 19 for µ = t and ν = i, one also finds
¯ α α
h̃ti = h̄ti − ∂t ξi − ∂i ξt = h̄ti − ikt Bi eikα x − iki Bt eikα x , (26)
¯
from which it follows that we can choose Bi such that h̃ti = 0. Moreover, the condition 25 for ν = i
¯
implies Azi = 0 (since h̃ti = 0 and thus Ati = 0). Evaluating the same condition for ν = t, one has
instead Att = Azt = 0. Looking the gauge transformation for the trace,
α
h̃ = h − 2∂µ ξ ν = h − 2ikµ B µ eikα x , (27)

one can finally see that by an appropriate choice of B t one can set h̃ = 0.
In light of the above, the residual gauge freedom allows one to write the polarization tensor as

0 0 0 0
 
0 h+ h× 0
(Aµν ) =  . (28)
0 h× −h+ 0
0 0 0 0

The conclusion is that gravitational waves propagating on flat space have only two independent transverse
polarizations, i.e. there exist only two propagating degrees of freedom.

1 This terminology derives from the fact that the d’Alembertian is the Minkowskian generalization of the (Euclidean)
Laplacian ∇2 = δ ij ∂i ∂j . In vector analysis, a function that satisfies the Laplace equation ∇2 f = 0 is called harmonic.
This is because the eigenfunctions of the Laplacian on the sphere are called “spherical harmonics”, since they are the
higher dimensional analog of the Fourier basis, consisting of sines and cosines, which describes harmonic motion.
2 This is also obvious because as mentioned above, the harmonic gauge condition can also be written as 2xµ = 0.
 5

B. Linear perturbations on curved space

Let us now generalize the previous calculation to a generic curved background. Again, the spacetime
metric is given by the background metric gµν and a perturbation hµν , and the trace reversed perturbation
can be defined as h̄µν = hµν − 21 hgµν , with h = hµν g µν . Indices are understood to be raised and lowered
with the background metric. Choosing the Lorenz gauge condition ∇µ h̄µν = 0, where ∇ is the covariant
derivative compatible with the background metric gµν , and linearizing the Einstein equations, one finds [2]

2h̄αβ + 2Rµανβ h̄µν + Sµανβ h̄µν = −16πT αβ , (29)

where 2 = g µν ∇µ ∇ν , Rµνρσ is the (background) Riemann tensor and
Sµανβ = 2Gµ(α gβ)ν − Rµν gαβ − 2gµν Gαβ . (30)
In vacuum Rµν = Tµν = Gµν = Sµανβ = 0, and Eq. 29 becomes

2h̄αβ + 2Rµανβ h̄µν = 0. (31)

Note in particular the coupling between gravitational waves and the background’s Riemann tensor, which
affects their propagation (i.e. gravitational waves can be scattered by the curvature and can therefore
also travel inside an observer’s lightcone).
Let us again count the physical (propagating) degrees of freedom in vacuum. After fixing the gauge,
hµν still has, in principle, 6 independent components. However, proceeding like in the flat background
case, one easily finds that in vacuum the Lorenz gauge leaves a residual gauge freedom, namely one can
always perform a gauge transformation with a harmonic generator (i.e. a generator ξ µ obeying 2ξ µ = 0)
and still preserve the harmonic gauge condition. Under a gauge transformation, the metric perturbation
transforms as
h̃µν = hµν − Lξ gµν = hµν − ∇µ ξν − ∇ν ξµ , (32)

hence the trace transforms according to h̃ = h − 2∇µ ξ µ . Let us first try to set the trace h̃ = 0 using the
residual gauge freedom of the Lorenz gauge. That would require choosing ξ µ such that
h
∇µ ξ µ − =0 (33)
2
throughout the whole spacetime.
To see if this requirement is compatible with the Lorenz gauge condition’s residual freedom, let us
note that the generators of the latter must obey 2ξ µ = 0 and are therefore completely characterized
by the initial conditions to this wave equation, i.e. ξ µ (t = 0, xi ) and ∂t ξ µ (t = 0, xi ). It is not a priori
obvious that by choosing these initial conditions properly, Eq. 33 can be satisfied in the whole spacetime.
However, taking a d’Alembertian of Eq. 33, an involved calculation using the trace of Eq. 31 and 2ξ µ = 0
yields
 
h
2 ∇µ ξ −µ
= 0. (34)
2
Therefore, for Eq. 33 to be satisfied in the whole spacetime, we just need to impose ∇µ ξ µ − h/2 =
∂t (∇µ ξ µ − h/2) = 0 at t = 0. This can be attained by choosing the initial conditions ξ µ (t = 0, xi ) and
∂t ξ µ (t = 0, xi ) characterizing the residual gauge freedom [3].
In conclusion, also in curved spacetime the trace of the metric perturbations can be set to zero in
the Lorenz gauge. However, showing that only two non-zero components (h+ and h× ) survive, like in
Minkowski space, is in general not possible. We will shed light on this fact in the next section, where we
will do perturbation theory in a slightly different way, by exploiting a scalar-vector-tensor-decomposition
of the metric perturbation. As we will see, there will still be just two propagating degrees of freedom for the
gravitational field, but additional non-propagating potentials will be present, including and generalizing
the Newtonian potential.

C. Linear perturbations on flat space: a scalar-vector-tensor decomposition

To gain more insight on the degrees of freedom of the metric perturbation, let us go back to the case of a
flat background spacetime. Let us introduce a book-keeping parameter   1 and write gµν = ηµν + hµν .
 6

Moreover, let us describe matter by the perfect fluid stress energy tensor

T µν =  [(ρ + p)uµ uν + pg µν ] . (35)

One can then split the components of the metric perturbation according to their transformation prop-
erties under spatial rotations. For instance, htt is a scalar under rotations, hti a vector, hij a tensor.
Moreover, we can perform a Helmholtz decomposition of the vector into the gradient of a scalar plus a
divergenceless vector (i.e. a curl), and similarly decompose the tensor into two scalars, a divergenceless
vector and a transverse traceless tensor. As a result, one has [3–5]

htt = 2φ,
with ∂i β i = 0,
hti = ∂i γ + βi
(36)
 
1 1
hij = Hδij + ∂(i εj) + ∂i ∂j − δij ∇2 λ + hTT
ij ,
3 3
with ∂i εi = 0 and ∂i hTTij = 0 = hTTii .

Here, spatial indices are understood to be raised and lowered with the Euclidean metric δij ; φ, γ, H
and λ are scalars under spatial rotations; β and ε are divergenceless vectors and hTT is a transverse
(∂i hTTij = 0) and traceless tensor. One can easily verify that the number of degrees of freedom of this
decomposition is correct. For instance, hti has three independent components, which correspond to γ
(one degree of freedom) and βi (two degrees of freedom). Similarly, hij has 6 independent components,
which correspond to the scalars λ and H (one degree of freedom each), the divergenceless vector εi (two
degrees of freedom) and the transverse traceless tensor hTTij (two degrees of freedom). Moreover, one
can verify that these decompositions are uniquely defined (up to boundary conditions). For instance, to
obtain γ we can compute

∂i hti = ∇2 γ + ∂i βi = ∇2 γ, (37)

and we can formally invert this expression to obtain

γ = ∇−2 (∂i hti ) . (38)

Here ∇−2 is the inverse of the (Euclidean) Laplacian, which is well defined if boundary conditions are
given for γ (e.g. it is reasonable to assume that γ decays “fast” at large distances to preserve asymptotic
flatness) and which can be expressed explicitly in terms of a Green function (see below). Once γ is
determined, βi can be computed as

βi = hti − ∂i ∇−2 (∂i hti ) .

 
(39)

Similarly, one can show that the decomposition of hij is well defined and unique by computing hii , ∂i hij
and ∂i ∂j hij .
A similar decomposition can be performed on the stress energy tensor [3],

Ttt = ρ,
Tti = ∂i S + Si with ∂i S i = 0,
(40)
 
1 2
Tij = pδij + ∂(i σj) + ∂i ∂j − δij ∇ σ + σij ,
3
with ∂i σ i = 0 and ∂i σ ij = 0 = σ ii ,

where ρ, S, p, σ are scalars, Si and σi divergenceless vectors and σij a transverse traceless tensor.
Similarly, the generator ξ µ of infinitesimal coordinate transformations can be expressed as

ξ t = A,
(41)
ξ i = ∂i C + Bi , with ∂i B i = 0

with A and C scalars and Bi a divergenceless vector.

 7

By using this decomposition for the generator in Eq. 12, we obtain [3]

φ̃ = φ − ∂t A , (42)
β̃i = βi − ∂t Bi , (43)
γ̃ = γ − A − ∂t C , (44)
H̃ = H − 2∇2 C , (45)
λ̃ = λ − 2C , (46)
ε̃i = εi − 2Bi , (47)
h̃TT
ij = hTT
ij . (48)

First, let us notice that hTT

ij is gauge invariant. Moreover, one can remove two scalars and one diver-
genceless vector by a suitable choice of A, C and Bi . For instance, one can choose to remove γ, λ and εi ,
so that
htt = 2φ,
hti = βi , (49)
1
hij = Hδij + hTT
ij .
3
This particular choice is called Poisson gauge, and unlike the Lorenz gauge, it fixes completely the
coordinates at linear order (i.e. there is no residual gauge freedom).
Alternatively, one can construct particular combinations of the scalar and vector potentials that are
gauge invariant (recall that hTT
ij is already gauge invariant). These are Bardeen’s gauge invariant vari-
ables [5], i.e.
1
ψ = −φ + ∂t γ − ∂t2 λ ,
2
1 2


θ= H −∇ λ , (50)
3
1
Σi = βi − ∂t εi ,
2
which reduce respectively to −φ, H/3 and βi in the Poisson gauge. Thus, using the Poisson gauge is
exactly equivalent to using Bardeen’s gauge invariant variables.
We can now express the linearized Einstein equations in terms of the Bardeen variables (or alternatively
compute them in the Poisson gauge). For the Einstein tensor we obtain [3]

Gtt = −∇2 θ,
1
Gti = − ∇2 Σi − ∂i ∂t θ,
2 (51)
 
1 TT 1 1 2
Gij = − 2hij − ∂(i ∂t Σj) − ∂i ∂j (2ψ + θ) + δij 2
∇ (2ψ + θ) − ∂t θ ,
2 2 2

with 2 the flat d’Alembertian. Decomposing also the right-hand side (i.e. the stress energy tensor), we
find that for the tt-component of the equations one has

Gtt = 8πTtt ⇔ ∇2 θ = −8πρ. (52)

From the ti-components one then has

 
1
0 = Gti − 8πTti = (−∂i ∂t θ − 8π∂i S) + − ∇2 Σi − 8πSi . (53)
2

where the first term in round brackets is the scalar part and the second is the (divergenceless) vector
part. Since the Helmholtz decomposition of a vector is unique, in order for this equation to be satisfied
both terms must vanish, i.e.
(
∂t θ + 8πS = 0 ,
1 2 (54)
2 ∇ Σi + 8πSi = 0 .
 8

The same procedure can be applied to the spatial components:

1
0 = Gij − 8πTij = − (2hTT + 16πσij ) − [∂(i ∂t Σj) + 8π∂(i σj) ]
 2 ij  
θ 1 2 8
− ∂i ∂j ψ + + 8πσ + δij ∇ (2ψ + θ) − ∂t2 θ + π∇2 σ − 8πp , (55)
2 2 3

leading to

2hTT



 ij = −16πσij ,
∂ Σ + 8πσ = 0 ,
t j j
θ (56)
ψ + + 8πσ = 0,
 2 2


 θ

3 2
∇ ψ + 2 − 2 ∂t θ − 12πp = 0 .

Let us consider now the energy conservation and relativistic Euler equations, which we have seen to
follow from the conservation of the matter stress-energy tensor, ∂µ T µν = 0. Decomposing this vector
equations in two scalar equations and one equation for a divergenceless vector in the (by now) usual way,
one gets [3]

∇2 S = ∂t ρ,
3 3
∇2 σ = − p + ∂t S, (57)
2 2
∇2 σi = 2∂t Si .

These equations can be used to simplify the Einstein equations 52, 54 and 56. In particular, as expected
from the Bianchi identify, one can show explicitly that three of those equations (one involving a diver-
genceless vector and two involving scalars) are automatically satisfied on shell (i.e. if the matter stress
energy conservation is enforced). The remaining Einstein equations can then be written as [3]


 ∇2 θ = −8πρ (1 d.o.f.) ,
∇2 ψ = 4π (ρ + 3p − 3∂ S)

(1 d.o.f.) ,
t
2
(58)

 ∇ Σ i = −16πS i (2 d.o.f.’s) ,
2hij = −16πσij

 TT
(2 d.o.f.’s) .

These are six independent equations for the six gauge invariant degrees of freedom of the metric per-
turbation. As can be seen, only the transverse traceless tensor modes (i.e. the gravitational waves) are
propagating, while the scalar and vector modes are not, as they satisfy elliptic (Poisson-like) equations.
The solution of the Poisson equation can be easily written in terms of the Laplacian’s Green function.
Let us recall indeed the distributional identity,
1
∇2 = −4πδ (3) (~x) , (59)
|~x|

which implies that the Laplacian’s Green function is proportional to 1/|~x − ~x0 |. We can then write, for
instance,

ρ (t, ~x0 ) 3 0
Z
θ (t, ~x) = 2 d x , (60)
|~x − ~x0 |

which resembles the Newtonian potential (as we will further discuss later on). Far from a localized source,
we can then approximate |~x − ~x0 | ≈ |~x| = r and write
Z
2
θ (t, ~x) ≈ d3 x0 ρ (t, ~x0 ) . (61)
r

It is tempting to call “mass” the integral of ρ, M = d3 xρ. We can in fact see that this quantity is
R
conserved:
Z Z
d
ρ d3 x = ∂t ρ d3 x, (62)
dt
 9

and using the conservation of the stress-energy tensor (Eq. 57) one has
Z  Z Z
d ~ · ~n d2 S,
ρ d3 x = ∇2 S d3 x = ∇S (63)
dt

where in the last equality we used Gauss’ theorem to reduce the integral to the flux of ∇S ~ through a
~
surface at infinity with normal unit vector ~n. If there is no matter at infinity, ∇S vanishes, proving the
conservation of the “mass”. Solving in a similar way the equation for ψ one gets
Z
1 M̄
ψ (t, ~x) ≈ − d3 x0 (ρ + 3p − 3∂t S) = , (64)
r r

with M̄ = d3 x (ρ + 3p − 3∂t S). Using the conservation of the stress-energy tensor (Eq. 57), we can
R
prove that
Z   Z
M̄ − M = d x 3p − 3Ṡ = −2 ∇2 σ d3 x0 = 0,
3 0
(65)

where in the last step we used again the absence of matter at infinity. This proves that M = M̄ . For Σi
we obtain
Z
4
Σi ≈ Si d3 x0 , (66)
r
and with similar steps one finds that
Z
πi = Si d3 x0 , (67)

which physically describes the linear momentum of the source, is conserved.

In summary, the scalar and vector gauge invariant Bardeen variables present a Newton-like behavior,
i.e.
M
ψ∼ ,
r
M (68)
θ∼ ,
r
πi
Σi ∼ .
r
As will become clearer from the post-Newtonian (PN) formalism, these three non-propagating degrees
of freedom generalize the Newtonian potential (ψ), and encode relativistic effects such as periastron
precession and light bending (θ) and frame dragging (Σi ).

D. Generation of gravitational waves: a first derivation of the quadrupole formula

Let us consider now not the propagation, but the generation of gravitational waves from matter sources,
by solving

hTT
ij = −16πσij . (69)

This will lead us to a first derivation of the quadrupole formula. We will then highlight some shortcomings
of this derivation, which will be amended in section III.
The solution of Eq. 69 can be obtained in terms of retarded potentials. The Green function of the flat
space d’Alembertian 2 = −∂t2 + ∇2 is
1
G (t, ~x) = − δ (t − |~x|) , (70)
4π |~x|

which indeed satisfies the distributional identity

2G (t, ~x) = δ(t)δ (3) (~x) . (71)

 10

The solution can then be written as

σij (t − |~x − ~x0 | , ~x0 ) 3 0
Z
hTT
ij (t, ~x) = 4 d x. (72)
|~x − ~x0 |

In order to find the source σij from the stress-energy tensor, we have to invert the Eq. 40. To this purpose
we can formally define the projector

Pij = δij − ∇−2 ∂i ∂j , (73)

and write
 
1
σij = Pi k Pj l − Pij P kl Tkl . (74)
2

Using the fact that partial derivatives commute and thus ∇−2 ∂i = ∂i ∇−2 , one can show that Eq. 74
implies ∂i σ ij = 0 and σ ii = 0, i.e. Eq. 74 correctly defines the transverse and traceless part of the matter
stress energy tensor. We can therefore write
 
TT −1 −1 k l 1 kl
hij = −16π2 σij = −16π2 Pi Pj − Pij P Tkl . (75)
2

Because the flat d’Alembertian and partial derivatives commute, we can then proceed to write
 
1
hTT
ij = −16π P i
k
P j
l
− P ij P kl
2−1 Tkl =
2
Tkl (t − |~x − ~x0 | , ~x0 ) 3 0
 Z
1
= 4 Pi k Pj l − Pij P kl d x = (76)
2 |~x − ~x0 |
 Z   
4 1 1
= Pi k Pj l − Pij P kl Tkl (t − |~x − ~x0 | , ~x0 ) d3 x0 1 + O ,
r 2 r

where in the last step, besides approximating |~x − ~x0 | with r = |~x|, we have also commuted 1/r with the
projectors, which is appropriate at leading order in 1/r. Moreover, in the last step we can approximate,
up to subleading terms in 1/r,

Pij ≈ δij − ni nj , (77)

where ni = xi /r is a unit vector in the direction of the observer. Defining Pij kl = Pi k Pj l − 21 Pij P kl , one
can write this “Green” formula in more compact form as
Z
4
hTT
ij = P kl
Tkl (t − |~x − ~x0 | , ~x0 ) d3 x0 . (78)
r ij
To go from this equation to the quadrupole formula, one can note that from the conservation of the
stress energy tensor (in flat space) it follows that

∂t2 T tt xi xj = 2T ij + ∂k ∂l T kl xi xj − 2∂k T ik xj + T kj xi .




(79)

Using this equation, and neglecting surface terms that vanish if the source is confined, Eq. 78 then
becomes
Z
2
hTT kl
∂t2 T tt x0k x0l d3 x0 .


ij = P ij (80)
r
Defining the inertia tensor
Z
Iij = d3 x0 ρ x0i x0j (81)

and the quadrupole tensor

1
Qij = Iij − Iδij , (82)
3
 11

where I = I ii , we finally arrive at the “quadrupole formula”

2G
hTT
ij = P kl Q̈kl , (83)
c4 r ij

where ˙ = d/dt and we have reinstated G and c for physical clarity.

While this final result looks reasonable, two key assumptions were used to derive it, namely (i) linear
perturbation theory and (ii) the conservation of the stress energy tensor on flat space. Neither of these
assumptions is justified for a compact binary system, as (i) the spacetime is not a perturbation of
Minkowski space near black holes or neutron stars, and (ii) ∂µ T µν = 0 implies the geodesic equation
in flat space (cf. section I), which in turn implies straight line motion (which clearly cannot describe
quasicircular binary systems). In fact, when applying the Green and quadrupole formulae to binary
systems one gets into paradoxes such as that described in the next exercise. This shows that a better,
more rigorous treatment of gravitational wave generation is needed, which will prompt us to go beyond
linear theory in the next section.

Exercise 2: Consider an equal-mass, Keplerian binary (i.e. a binary with large separation, for
which the laws of Newtonian mechanics are applicable) on a circular orbit on the (x, y) plane, and an
observer (far from the source) along the z axis. Compute the gravitational-wave signal according to the
“Green formula” of the lectures, and according to the quadrupole formula. Show that the amplitudes of
the two predictions differ by a factor 2.

E. Dimensional analysis

Let us try to derive the quadrupole formula from dimensional arguments. A matter source can be
3
the mass dipole Di =
R
characterized by its multipole moments, e.g. the mass monopole M = ρ d x,
i 3
R
ρ x Rd x, the mass quadrupole Qij , the angular momentum (i.e. the first moment of the mass current)
Li = ρ eijk xj v k d3 x etc. Since metric perturbations are dimensionless, one can try to write a monopole
gravitational wave signal using dimensional analysis as h ∼ GM/(rc2 ), a dipole signal as h ∼ GḊ/(rc3 ) ∼
GP/(rc3 ) (where P is the linear momentum), an angular momentum term h ∼ GL̇/(rc4 ). These terms
are zero (or static) because of conservation of mass, linear momentum and angular momentum. The
quadrupole term, again by dimensional analysis, is instead h ∼ GQ̈/(rc4 ). Note that radiation sourced
by the mass monopole and dipole and by the angular momentum can be present beyond GR, because in
that case the mass, linear momentum and angular momentum of matter may not be conserved (due to
exchanges with additional gravitational degrees of freedom different from the tensor gravitons). Similarly,
the static scalar and vector degrees of freedom of Eq. 68 will generally become dynamical beyond GR.

III. POST-NEWTONIAN EXPANSION

In order to assess which of the two expressions for the generation of gravitational waves (the quadrupole
formula or the “Green formula”) is correct, let us take a small detour. We will now study perturbations
of flat space not by expanding in the perturbation amplitude (like we did previously), but in powers of
1/c (with c → ∞). This is known as post-Newtonian (PN) expansion, and will allow us to re-derive the
quadrupole formula in a more rigorous way.

A. The motion of massive and masseless bodies

Let us start by writing the following ansatz for the metric:

 
2φ
g00 = − 1 + 2 ,
c
ωi
g0i = 3 , (84)
c
 
2ψ χij
gij = 1 − 2 δij + 2 ,
c c
 12

where χij is traceless (χii = 0) and we have used Cartesian coordinates xµ = (ct, xi ). Latin indices are
meant to be raised and lowered with the flat spatial metric δij . Note that we have reinstated c (which
was set to 1 in the previous sections), as that will be our book-keeping parameter. Before venturing into
the actual calculation, let us note that the choice of powers of c appearing in Eq. 84 is exactly the one
that will be needed to consistently solve the Einstein equations (i.e. should we choose different powers,
the Einstein equations would set the potentials to zero, or they would not allow for a consistent solution).
However, it is possible to make sense of this ansatz also in a more physical way.
Let us consider a point particle moving in the geometry described by Eq. 84. From R the point particle
action (Eq. 7), one can obtain the Lagrangian (by recalling that by definition S = Ldt). By replacing
therefore Eq. 84 into Eq. 7, one obtains
Z Z
dτ
S = −mc2
p
dt = −mc −gµν ẋµ ẋν dt
dt
Z s
2ωi v i χij v i v j
 
2ψ
= −mc c2 + 2φ − 2 − 1 − 2 v 2 − dt
c c c2
1 v2 φ2 φv 2 ψv 2 v4 vi ω i χij v i v j
Z  
2 φ
≈ −mc 1+ 2 − − 4+ 4 + 4 − 4− 4 − dt, (85)
c 2 c2 2c 2c c 8c c 2c4

where v i ≡ ẋi ≡ dxi /dt, v 2 ≡ δij v i v j , and in the last step we have Taylor expanded in 1/c (for c → ∞).
The Lagrangian for a point particle therefore reads
 2  2
φv 2 ψv 2 v4 vi ω i χij v i v j
 
v φ
L = −mc2 + m −φ +m − − + + + + ... (86)
2 2c2 2c2 c2 8c2 c2 2c2

As can be seen, in the limit c → ∞ the ansatz of Eq. 84 leads to the correct Newtonian limit (note
the appearance of the Newtonian Lagrangian after the irrelevant constant offset term), as well as to
deviations from the Newtonian Lagrangian that are suppressed by O(1/c2 ) relative to the Newtonian
dynamics. These are known as 1PN corrections (where nPN denotes terms that are suppressed by 1/c2n
relative to the leading order Newtonian term).
Note however that this counting is only applicable to massive particles, and not to photons (for which
v ∼ c). In the latter case, the terms ψv 2 /c2 , φv 2 /(2c2 ) and χij v i v j /(2c2 ), which are of order 1PN for a
massive particle, are of order 0PN (i.e. Newtonian order). In more detail, one can write the Lagrangian
from photons as L ∝ dτ /dt, which using Eq. 84 leads to
r
φ 2ωi β i ψβ 2 χij β i β j
L ∝ 1 − β2 + 2 2 − 3
+2 2 −
c c c c2  
p γ 1
= 1 − β 2 + 2 2φ + 2ψβ 2 − χij β i β j + O 3 ,


(87)
2c c
p
with β i = v i /c and γ = 1/ 1 − β 2 . In other words, the bending of light in GR is determined at leading
order by both g00 (φ) and gij (ψ and χij ), but not by g0i (ωi ). The same can be seen, e.g., by using
Eq. 84 into the dispersion relation for a photon, pµ pν gµν = 0, with pµ = (E, pi ) the 4-momentum, and
solving for the energy E to derive the Hamiltonian describing the photon’s motion.

B. The Einstein equations

Let us now compute the potential appearing in the metric from the Einstein equations. To do so, let
us first perform the usual scalar-vector-tensor decomposition on the metric ansatz of Eq. 84:

ωi = ∂i ω + ωiT ,
(88)
 
1
χij = ∂i ∂j − δij ∇ χ + ∂(i χT
2 TT
j) + χij ,
3

where the index “T” identifies divergenceless (i.e. transverse) vector fields and “TT” transverse and
traceless tensors. Let us first adopt the same gauge that we used in linear theory, i.e. the Poisson gauge,
defined by

0 = ∂i ω i = ∂i χij , (89)
 13

which yields ω = χ = χTi = 0. Let us also describe matter as a perfect fluid with stress energy tensor

T µν = (p + ρc2 )uµ uν + pg µν , (90)

with ui /u0 = v i /c [and therefore u0 = 1−(φ−v 2 /2)/c2 +O(1/c4 ) because of the 4-velocity normalization].
Using then the Einstein equations, in which we reinstate c to obtain Gµν = 8πTµν /c4 , one gets the
following equations for the potentials [6]:
 
1
ψ =φ+O 2 , (91)
c
 
1
∇2 ωTi = 4(4πρv i + φ,ti ) + O 2 , (92)
c
 p  v 2  
 2 3 1
∇2 φ = 4π 3 2 + ρ + 2 φ,i φ,i + 8πρ − 2 φ,tt + O 4 , (93)
c c c c c
 
1
∇2 χTijT = O 2 . (94)
c

As a consistency check, note that by taking the divergence of Eq. 92, both sides evaluate to zero: the
left hand side because ωTi is transverse, and the right hand side because of the continuity equation for
the number density, which at leading (Newtonian) order reads ∂t ρ + ∂i (ρv i ) = ∂t ∇2 φ/(4π) + ∂i (ρv i ) = 0.
(This is because the rest mass density and the energy density differ by the internal energy p/[c2 (Γ − 1)],
with Γ the adiabatic index; see e.g. [6].)
One can then write
φ1PN
φ = φN + + ..., (95)
c2
ψ2PN
ψ = φN + 2 + . . . , (96)
c
i i ωi
ωT = ω1PN + 2PN + ..., (97)
c2
χij χij
χij
TT = 2PN
+ 2.5PN
+ ..., (98)
c2 c3
where φN is the Newtonian potential (obtained by solving ∇2 φN = 4πρ), and we have left indicated
the terms that appear at 1PN order and higher in the Lagrangian for massive particles derived in the
previous section. These PN terms can be obtained explicitly by solving Eqs. 91–94 and their higher order
generalizations. In particular, one can show explicitly that the leading order term for χij T T appears at
O(1/c2 ) (2PN order in the Lagrangian for massive particles). This is a conservative term (as it is of even
parity in c, i.e. it is left unchanged by a time reversal). However, a dissipative term appears at O(1/c3 ),
i.e. at 2.5 PN order. This corresponds to the loss of energy and angular momentum to gravitational
waves (see e.g. [6] for details).
One unsightly feature of the Poisson gauge is however apparent from Eq. 93, which features a double
time derivative of φ on the right-hand side. That term corresponds to ∂t S in Eq. 58, i.e. one can
re-express it in terms of the matter density by writing it as −3φ,tt /c2 = −12π∇−2 ρ,tt /c2 + O(1/c4 ) =
−12πS,t /c2 + O(1/c4 ). That requires, however, solving a non-local equation to compute ∇−2 ρ. A better
option is to eliminate the term −3φ,tt /c2 by performing a gauge transformation with generator ξ0 ∝ ∂t X,
where X = −2∇−2 φ is the Newtonian “superpotential” [7] (see also the appendix of [9]). This leads to
the “standard PN gauge”, which is defined exactly as a gauge in which the 1PN spatial metric is isotropic
(i.e. χij is zero at 1PN, which we have seen to be already the case in our Poisson gauge) and in which
no term proportional to φ,tt appears at 1PN in the equation for ∇2 φ [7, 8].
Even more simply, one can do the calculation starting directly in the standard PN gauge, which satisfies
the gauge conditions [7]
1
∂µ hµi − ∂i hµµ = 0 , (99)
2
µ 1 1
∂µ h0 − ∂0 hµµ = − ∂0 h00 , (100)
2 2
where gµν = ηµν + hµν and the indices of hµν are raised and lowered with the Minkowski metric. In this
 14

gauge, the Einstein equations at 1PN order become

φ2 v2
   
Φ2 φ p
∇2 φ − 2 + 4 2 = 4π ρ + 2ρ 2 + 2ρ 2 + 3 2 , (101)
c c c c c
∇2 ψ = 4πρ , (102)
2
∇ ωj = 16πρvj + ∂t ∂i φ (103)

with ∇2 Φ2 = 4πρφ. Solving these equations by using the Green function of the flat Laplacian one gets
the 1PN metric in the standard PN gauge as [7]

φN φ2 Φ1 Φ2 Φ4
g00 = −1 − 2 2
−2 N +4 4 +4 4 +6 4 (104)
c c4 c c c
7 Vi 1 Wi
g0i = − 3 − (105)
2 c 2c3
φN
gij = 1 − 2 2 δij (106)
c

in terms of the PN potentials

ρ(~x0 , t)vi (~x0 , t)

Z
Vi = d3 x0 , (107)
|~x − ~x0 |
ρ(~x0 , t)[~v (~x0 , t) · (~x − ~x0 )](x − x0 )i
Z
Wi = d 3 x0 , (108)
|~x − ~x0 |3
ρ(~x0 , t)v(~x0 , t)2
Z
Φ1 = d3 x0 , (109)
|~x − ~x0 |
ρ(~x0 , t)φN (~x0 , t)
Z
Φ2 = − d3 x0 , (110)
|~x − ~x0 |
p(~x0 , t)
Z
Φ4 = d3 x0 . (111)
|~x − ~x0 |

Note that to obtain this result we have used the relation [7] ∂t ∂i X = Wi − Vi (where ∇2 X = −2φN ,
as defined previously), which follows from the explicit expression for the Newtonian superpotential,
X = d3 x0 ρ(~x0 , t) |~x − ~x0 |, and from the continuity of the number density [∂t ρ = −∂i (ρv i ) at Newtonian
R
order]. While the choice of the standard PN gauge bears no physical significance (the Poisson gauge
has the same physical validity, since observables in general relativity are gauge invariant), Eq. 104 is
the metric usually adopted to describe tests of general relativity in the solar system (e.g. periastron
precession, light bending, Shapiro time delay, lunar laser ranging, frame dragging, etc; see [7] for more
details).

C. A more rigorous derivation of the quadrupole formula

By using now the PN expansion in place of the linear approximation, let us revisit the generation of
gravitational waves from binary systems. This will lead us to re-derive the quadrupole formula in a more
rigorous fashion, which will in turn shed light on the discrepancy between quadrupole formula and “Green
formula”, which we discovered in Exercise 2. As previously mentioned, the problem with the linear theory
derivation of the quadrupole formula is two-fold: the assumption of “weak gravity” (hµν  1) and the
use of the stress energy tensor conservation in flat space. Here, we will fix both of these shortcomings.
To drop the weak gravity assumption, let us start from the full Einstein equations, which we write in
“relaxed form” in the harmonic gauge, defined by

2xα = 0 . (112)

It is important to keep in mind that the coordinates xα , despite the space-time indices, are not vectors
but scalars, as can be seen from their transformation properties under diffeomorphisms. Using this fact,
 15

it is straightforward to see that the condition 112 can be rewritten in terms of the pseudo-tensor3 :
√
H̄ µν = η µν − −gg µν . (113)

Expanding then
1 √ 1 √
2xα = √ ∂µ −gg µν ∂ν xα = √ ∂µ −gg µν δνα =



−g −g
(114)
1 √
−gg µα ∝ ∂µ H̄ µν ,


= √ ∂µ
−g

the condition 112 turns out to be equivalent to

∂µ H̄ µν = 0 . (115)

Note that the quantity H̄ µν becomes the trace reversed metric perturbation at linear order. Indeed, if
g µν = η µν − hµν + O(h2 ), then δg = −h + O(h2 ), and H̄ µν = hµν − 21 hµν + O(h2 ) = h̄µν + O(h2 ). As a
result, at linear order the harmonic gauge condition 115 coincides with the Lorenz gauge condition 13.
In the harmonic gauge, the fully non-linear Einstein equations take the form

2η H̄ µν = −16πτ µν , (116)

where 2η = η µν ∂µ ∂ν is the flat space d’Alembertian operator and

Λµν
τ µν = (−g)T µν + , (117)
16π
with

Λµν = 16π(−g)tµν µα
∂α H̄ νβ − ∂α ∂β H̄ µν H̄ αβ .


LL + ∂β H̄ (118)

Here, tαβ
LL is the Landau-Lifshitz pseudo-tensor,

16π(−g)tαβ
LL
≡ gλµ g νρ H̄,ν
αλ βµ
H̄,ρ (119)
1 β)ν ρµ
+ gλµ g αβ H̄,ρ
λν ρµ
H̄,ν − 2gµν g λ(α H̄,ρ H̄,λ
2
1
+ (2g αλ g βµ − g αβ g λµ )(2gνρ gστ − gρσ gντ )H̄,λ
ντ ρσ
H̄,µ ,
8
which describes the stress-energy of the gravitational field. Because it is a pseudo-tensor, it can be non-
zero in a set of coordinates but vanishing in a different one. This is simply a consequence of the well
known fact that in general relativity the gravitational field can be locally set to zero (by choosing Riemann
normal coordinates where the metric is locally ηµν and the Christoffel symbols vanish, c.f. section IV).
Indeed, there is no way of defining a local energy density for the gravitational field: only the global energy
(or mass) of an asymptotically flat spacetime is well defined in general relativity.
Taking now a (partial derivative) divergence of Eq. 116 and using the condition 115, one obtains the
conservation law ∂µ τ µν = 0, which is equivalent to the equations of motion of matter (c.f. section I).
Note that because we have not made any approximations thus far, these are the fully nonlinear equations
of motion of matter, i.e. unlike in the case of linear theory, we are not assuming straight-line motion or
∂µ T µν = 04 . We can then follow the same procedure of section II D to derive the quadrupole formula,
i.e. we can invert Eq. 116 as

H̄ µν = −16π2−1 τ µν , (120)

3 A pseudo-tensor is an object that transforms as a tensor under linear transformations, but not under more general
coordinate transformations.
4 The fact that ∂µ τ µν = 0 does not imply straight-line motion can be tracked back to the presence of second derivatives of
H̄ µν in Eq. 118: as a result, even though the left hand side of the relaxed Einstein equations is written in terms of the
flat wave operator, wavefronts do not follow straight lines in the eikonal approximation.
 16

which gives in particular the “Green formula”

Z
4  r 
H̄ ij (t, ~x) ≈ τ ij t − , ~x0 d3 x0 . (121)
r c
Using the conservation law ∂µ τ µν = 0 like in section II D, one can write
2
H̄ µν ≈ Q̈ij , (122)
r
where now
Z
Qij = τ tt x0i x0j d3 x0 . (123)

Let us now examine the relation between τ µν and T µν . Reinstating the appropriate powers of c, and
using the PN expanded metric of Eq. 84 and the stress energy tensor for a system of two point particles
[cf. Eq. 8], one finds that [10]
  
1
τ tt = T tt 1 + O 2 ,
c
  
1
τ ti = T ti 1 + O 2 , (124)
c
     
1 1 1
τ ij = T ij + ∂ i φ∂ j φ − δ ij ∂k φ∂ k φ 1+O 2 .
4π 2 c
Therefore, at leading PN order τ tt ≈ T tt , and Eq. 123 reduces to the quadrupole formula derived in
linear theory; however, τ ij cannot be approximated by T ij at leading PN order, i.e. the Green formula
121 does not reduce to the Green formula of linear theory. It follows that if one applies the formulae
derived in linear theory, only the quadrupole formula is correct.

Exercise 3: Show that the additional terms contributing to τ ij in Eq. 124 solve the factor 2 dis-
crepancy between the Green and quadrupole formula found in Exercise 2. [Hint: use the fact that the
solution to ∇2 g(x, y 0 , y 00 ) = |x − y 0 |−1 |x − y 00 |−1 (with ∇2 the Laplacian with respect to x) in the sense
of distributions is g = ln(|x − y 0 | + |x − y 00 | + |y 0 − y 00 |) + constant.]

IV. LOCAL FLATNESS AND THE EQUIVALENCE PRINCIPLE

In the previous section, when discussing the Landau-Lifshitz pseudotensor, we recalled that in general
relativity the gravitational field can always be locally set to zero, i.e. it is always possible to choose a
“Local Inertial Frame” where the gravitational force vanishes (i.e. where the metric is locally flat and the
Christoffel symbols vanish). This can be seen as a manifestation of the equivalence principle of general
relativity. In this section, we will provide a proof of this statement, which will clarify issues such as
that of the non-existence of a covariant stress-energy tensor for the gravitational field. The coordinates
that we introduce will also be useful when deriving the response of a gravitational wave detector in the
following.

A. The local flatness theorem and Riemann normal coordinates

The local flatness theorem states that at any given event (i.e. space-time point) P , there exists a
coordinate system such that gµν |P = ηµν and Γαµν |P = 0 (or equivalently ∂α gµν |P = 0). We will now
provide two proofs of this theorem.
Algebraic proof: Let us start with a system of coordinates {xα } such that P corresponds to xα = 0.
0
Let us then perform a coordinate transformation to some new coordinates {xα } also centered in P :
0 0 0
xα = Aαβ xβ + O(x2 ) ⇔ xα = Aαβ 0 xβ + O(x02 ), (125)
0
where Aαβ and Aαβ 0 are constant matrices (the Jacobian of the transformation and its inverse) that satisfy
0 0 0
Aαµ Aµβ 0 = δ αβ 0 and Aαµ0 Aµβ = δ αβ . (126)
 17

The metric at P transforms as

∂xα ∂xβ
gα0 β 0 |P = gαβ = gαβ Aαα0 Aββ 0 . (127)
∂xα0 ∂xβ 0 P

The matrix A has 16 coefficients, 10 of which can be chosen to set gα0 β 0 |P = ηα0 β 0 . The remaining 6
degrees of freedom correspond to the 6 generators of the Lorentz transformations, which are isometries
of the Minkowski metric.
In order to show that the Christoffel symbols vanish, let us expand the transformation to second order:
0 0 1 0
xα = Aαβ xβ + B αβγ xβ xγ + O(x3 ), (128)
2
0 0
where B αβγ is the Hessian of the transformation. Note that B αβγ has the same symmetries as the
Christoffel symbols, i.e. it is symmetric under the exchange β ↔ γ. As well known, the Christoffel
symbols are not tensors under generic coordinate transformations (otherwise it would not be possible to
set all of them to zero with a choice of coordinates, which is what we are trying to prove), but they
transform according to
0 0
0 ∂xβ ∂xγ ∂xα ∂xα ∂xβ ∂xγ
Γαβ 0 γ 0 = Γαβγ 0 0 − =
β
∂x ∂x ∂x γ α ∂xβ ∂xγ ∂xβ 0 ∂xγ 0 (129)
0 0
= Aαα Aββ 0 Aγγ 0 Γαβγ − B αβγ Aββ 0 Aγγ 0 .
0 0
We can therefore impose Γαβ 0 γ 0 =0 by solving for B αβγ . This equation has a unique solution, because B
and Γ share the same symmetries.
This concludes our first proof of the local flatness theorem. The coordinates where the latter holds
are known as “Riemann Normal Coordinates” (RNCs). We will now give a more geometric proof of the
theorem, which includes a procedure to explicitly construct these coordinates.
Geometric
0
proof: Let us consider a spacetime endowed with coordinates xµ , and explicitly construct
RNCs xµ around an event P . To assign coordinates to a neighboring point P 0 , let us consider the unique
geodesic connecting P and P 0 , and the vector v tangent to this geodesic in P . Let us decompose this
vector onto its components on a tetrad centered in P , i.e. on a basis of four orthogonal unit-norm vectors
{e(α) }α=1,...,4 :

v µ = Ω(α) eµ(α) . (130)

By definition, the tetrad vectors satisfy the orthonormality and completeness relations
e(α) · e(β) ≡ gµν eµ(α) eν(β) = ηαβ ,
(131)
eµ(α) e(α) µ
µ =δ ν,

where the tetrad (bracketed) indices are raised and lowered with the Minkowski metric ηαβ , whereas the
space-time indices are raised and lowered with the space-time metric gµν .
As a working hypothesis, let us choose the new coordinates of the point P 0 to be
0 0
xα = Ω(α ) ∆λ , (132)

with ∆λ = λP 0 − λP , where λ is the affine parameter of the geodesic connecting P and P 0 . First, let
us check that this definition is invariant under a re-parametrization of the geodesic. We know that the
geodesic equation is invariant under affine transformations of the parameter, λ0 = aλ + b. Under this
transformation, ∆λ0 = a∆λ, and
dxµ dxµ dλ vµ
v µ |λ0 = 0
= 0
= . (133)
dλ dλ dλ a
Thus, the coordinates of P 0 remain unchanged:
0
0 Ω(α ) 0
xα |λ0 = a∆λ = xα |λ . (134)
a
To see that the coordinates that we constructed are indeed RNCs, let us first check that the metric
at the event P is given by the Minkowski metric in the new coordinates. To this purpose let us first
 18

compute the Jacobian of the transformation from the old to the new coordinates evaluated at the point
P . From Eq. 132 one has
0
dxα 0
= Ω(α ) . (135)
dλ P

Therefore, one also has

0
dxα α ∂xα dxα ∂xα 0
=v = Ω(µ) eα
(µ) = = Ω(α ) . (136)
dλ P ∂xα0 P dλ P ∂xα0 P

From the arbitrariness of the components Ω(µ) one then gets

∂xα
= eα
(α0 ) , (137)
∂xα0 P

and therefore
∂xα ∂xβ β
gα0 β 0 |P = gαβ |P = gαβ |P eα
(α0 ) e(β 0 ) = ηα0 β 0 . (138)
∂xα0 P ∂x
β0
P

To compute instead the Christoffel symbols at P , let us consider a one parameter family of events0 along
the geodesic connecting P and P 0 . This family has coordinates growing linearly with λ, i.e. xα ∝ λ,
2 α0 2
which implies d x /dλ = 0. Since this one-parameter family is (by definition) a geodesic with affine
parameter λ, one must have
0 0 0
d2 xα 0 dxµ dxν 0 0 0
= −Γαµ0 ν 0 = −Γαµ0 ν 0 Ω(µ ) Ω(ν ) = 0. (139)
dλ2 P
P dλ
P dλ P
P

0
Since the procedure can be repeated for all geodesics originating from P , the components Ω(µ ) are
α0
arbitrary, from which it follows that Γ µ0 ν 0 P = 0, which completes the proof.
We can therefore conclude that around the event P , the metric in RNCs is gµ0 ν 0 = ηµν + O(x0 )2 . It is
possible to prove [11] that the quadratic terms O(x0 )2 are proportional to components of the Riemann
tensor, i.e. those terms, being dimensionless, scale as (x0 /L)2 , with L the curvature radius of the spacetime
(defined from the Riemann tensor).

B. Fermi Normal Coordinates

Let us now slightly modify the idea behind the geometric construction of RNCs to build a set of
coordinates describing the reference frame of an observer in motion along a generic timelike worldline γ
with 4-acceleration aµ . The coordinates, usually referred to as “Fermi normal coordinates” (FNCs) are
defined in a worldtube surrounding the worldline γ.
Let us consider a tetrad e(α) (with α = 1, . . . 4) attached to the wordline γ, with eµ(0) = uµ (i.e. the
“time” leg of the tetrad coincides with the 4-velocity of the worldline). The tetrad is assumed to be
Fermi-Walker transported along γ.5 For any point P 0 within a worldtube surrounding γ, let us consider
the unique (spacelike) geodesic that passes through P 0 and intersects orthogonally the trajectory γ. The
intersection point will be denoted by P . Let us assign to P 0 a new time coordinate
0
x0 = τ , (141)

5 A vector ω µ is said to be Fermi-Walker transported along a worldline with 4-velocity uµ , proper time τ and acceleration
aµ if
Dω µ
= −ων (aµ uν − aν uµ ) , (140)
dτ
with D/dτ the covariant derivative along the four velocity. This transport preserves angles and internal products, and
reduces to parallel transport when aµ = 0. Moreover, it is easy to check that the 4-velocity is always Fermi-Walker
transported, which explains why we can always choose the time leg of our Fermi-Walker transported tetrad to be the
4-velocity.
 19

where τ is the proper time of the worldline at P . As for the spatial coordinates, let us consider, at the
point P , the vector v tangent to the spacelike geodesic linking P and P 0 , and decompose it on the spatial
“triad” e(i) . (Note that the projection on e(0) vanishes, because the spacelike geodesic is constructed to
be orthogonal to γ.) Denoting by s the proper length along this spacelike geodesic, let us therefore write

v µ = Ω(i) eµ(i) , (142)

where the Ω(i) are the projections on the triad, and define the new spatial coordinates of the point P 0 as
0 0
xi = ∆s Ω(i ) . (143)
Proceeding like in the case of RNCs, it is easy to check that this definition is invariant under affine
reparametrizations of the geodesic connecting P and P 0 .
To explore the consequences of these definitions, let0 us compute the Jacobian of the transformation
between the old coordinates xµ and the new FNCs xµ , on the worldline γ. First, note that
∂xα
= uα = eα
(0) (144)
∂x00 γ
0
simply because x0 = τ . Again because of how we defined FNCs, we have
dxα
= v α = Ω(i) eα
(i) , (145)
ds γ

but also
0
dxα ∂xα dxi ∂xα 0
= = Ω(i ) . (146)
ds γ ∂xi0 γ ds γ ∂xi0 γ

Comparing the two, we then obtain

∂xα
= eα
(i0 ) . (147)
∂xi0 γ

Equations 144 and 147 can then be expressed concisely as

∂xµ
= eµ(µ0 ) . (148)
∂xµ0 γ

The metric at point P in FNCs is then

∂xµ ∂xν
gµ0 ν 0 |γ = gµν |γ = eµ(µ0 ) eν(ν 0 ) gµν |γ = ηµ0 ν 0 . (149)
∂xµ0 γ ∂x
ν0
γ

Therefore, FNCs ensure that the metric is Minkowskian on the whole worldline γ.
Let us now see what happens 0for the Christoffel
0
symbols 0 on γ. Since the spacelike geodesic connecting
P and P 0 is parametrized by x0 = const and xi = ∆s Ω(i ) in FNCs, one must have
0 0 0
d2 xµ 0 dxα dxβ 0 0 0
+ Γµα0 β 0 = Γµi0 j 0 γ Ω(i ) Ω(j ) = 0 , (150)
ds2 γ
γ ds
γ ds γ

0 0
and since the components Ω(i ) are generic, Γµi0 j 0 |γ = 0.
0 0
Let us then parametrize the worldline γ in FNCs, i.e. x0 = τ and xi = 0. Since γ is not necessarily a
geodesic, one has
0 0 0
µ0 d2 xµ 0 dxα dxβ 0
a = + Γµα0 β 0 γ = Γµ00 00 . (151)
dτ 2 γ dτ γ dτ γ
γ

Moreover, by inverting Eq. 148 one finds

0
∂xµ 0
= e(µ
µ
)
, (152)
∂xµ γ
 20
0
which allows for computing aµ :
0
0 ∂xµ 0
aµ = aµ = eµ(µ ) aµ , (153)
∂xµ γ

0
which yields a0 = 0 (since the time leg of the tetrad is the 4-velocity, which is orthogonal to the
0 0 (i0 )
acceleration because of the unit-norm condition) and ai = a(i ) = eµ aµ . From Eq. 151 one therefore
0 0 0
obtains Γ000 00 |γ = 0 and Γi00 00 |γ = a(i ) .
Finally,
0
let us consider a space-like unit-norm vector ω orthogonal to the worldline, i.e. in FNCs
ω µ = (0, Ω̄(i) ), with Ω̄(i) = const such that δij Ω̄(i) Ω̄(j) = 1. Its components in the original coordinates
are
∂xµ 0 0
ωµ = ω µ = eµ(µ0 ) ω µ = eµ(i0 ) Ω̄(i) . (154)
∂xµ0 γ

Because the components Ω̄(i) are constant and because the triad e(i) is Fermi-Walker transported along
γ, this vector is also Fermi-Walker transported along γ, and thus
0
Dω µ 0 0 0 0 0 0 0 0
= −(aµ uν − aν uµ )ων 0 = aν 0 ω ν uµ = a(i0 ) Ω̄(i ) δ0µ0 , (155)
dτ
0 0
where we have used the fact that uµ = δ0µ0 . However, by definition one also has
0 0
Dω µ dω µ 0 0 0 0 0
= + Γµα0 β 0 uα ω β = Γµ00 i0 Ω̄(i ) , (156)
dτ dτ
from which, by comparing to Eq. 155, one obtains
0 0
Γ000 i0 = a(i0 ) and Γj00 i0 = 0. (157)

In summary, collecting our findings, the Christoffel symbols on γ in FNCs are

0 0 0 0 0 0
Γµi0 j 0 = 0, Γ000 00 = 0, Γi00 00 = a(i ) , Γ000 i0 = a(i0 ) , Γj00 i0 = 0. (158)

These conditions completely determine ∂µ0 gα0 β 0 on the worldline, from which one can expand the metric
0
in the worldtube up to second order in xi to obtain
0
g00 00 = −1 − 2a(i0 ) xi + O(x0 )2 ,
g00 i0 = O(x0 )2 , (159)
0 2
gi0 j 0 = δi0 j 0 + O(x ) .

This is the metric “felt” by an observer on the trajectory γ. Like in the case of RNCs, the quadratic
remainders O(x0 )2 are actually proportional to 1/L2 , with L the curvature radius [11].
As can be seen, the metric is not locally flat unless the worldline is geodesic, i.e. unless the observer is
in free fall. If that is not the case, the acceleration enters the metric, and specifically g00 00 , as an apparent
Newtonian potential. This potential encodes the apparent forces due to the non-inertial motion of the
observer. On the other hand, if the acceleration vanishes, FNCs generalize RNCs to the entire worldtube.
Finally, note that the form of the metric in FNCs would change if we were to use a different transport
law for the triad e(i) , which can be interpreted as the spatial frame of the observer’s laboratory. If the
triad were not Fermi-Walker transported, additional terms would appear in the metric, corresponding to
more general apparent forces (centrifugal forces, Coriolis forces, etc).

V. THE STRESS ENERGY TENSOR OF GRAVITATIONAL WAVES

Let us now revisit the problem of defining the stress energy tensor of the gravitational field in GR. As
clear from the local flatness theorem, it is always possible to choose a coordinate chart in which the effect
of the gravitational field locally vanishes. In these local coordinates, the metric is therefore flat and local
experiments obey the laws of physics without gravity. Choosing these coordinates, as should be clear
 21

from the geometric construction of FNCs provided in the previous section, corresponds to adopting a
reference frame attached to a free falling observer, e.g. the frame of a free falling elevator. Their existence
is thus a manifestation of the well-known equivalence principle of GR.
Because the effect of the gravitational field can always be made to vanish with a proper choice of
coordinates, it is clear that it is not possible to define a stress energy tensor for the gravitational field,
because if a tensor is non-zero in a set of coordinates, it is non-zero in any other set of coordinates
connected to the first by a non-singular transformation. The most we can aspire to is therefore to build
a stress energy pseudotensor for the gravitational field in GR. We have already encountered this pseu-
dotensor, the Landau-Lifschitz pseudotensor, when deriving the relaxed form of the Einstein equations in
section III C. An alternative derivation exploits the fact that for a field theory defined by a Lagrangian
density L(ψn , ∂µ ψn ), where ψn are fields, the stress energy tensor in flat space can be defined as
X ∂L
T µν = ∂ν ψn − δ µν L , (160)
n
∂(∂µ ψn )

which satisfies the conservation equation ∂µ T µν = 0 because of the Lagrange equations. This stress
energy tensor is in general not symmetric in its indices (when the latter are both covariant or both
contravariant), and therefore needs to be redefined as T µν → T µν + ∂α S αµν , with a suitable tensor S αµν
satisfying S [αµ]ν = 0. As a result of this symmetry, the redefined tensor still satisfies ∂µ T µν = 0.
Because this construction only works for a Lagrangian depending on first derivatives of the fields, for
GR let us start from the Γ-Γ (or Schrödinger) action
Z
1
S=− dx4 L , (161)
16π
√  
L = − −gg αβ Γµαβ Γνµν − Γναµ Γµβν , (162)

which differs from the usual Einstein-Hilbert action by a surface term. One can then define the Einstein-
Schrödinger complex
 
µ 1 ∂L µ
t ν =− √ ∂ν gαβ − δ νL , (163)
16π −g ∂(∂µ gαβ )
√
which satisfies ∂µ [ −g(tµν + T µν )], where T µν is the matter stress energy tensor. Not only is this object
a pseudotensor and not a tensor, but it is also not symmetric in its two indices, when the lower index
is raised. The idea is therefore to redefine it as explained above to achieve symmetry between the two
indices.
In more detail, one can then show [12, 13] that
√
16π −g(tµν + T µν ) = ∂α S µαν (164)

µ[ασ]
where long algebraic manipulations give S µαν = U µαν + ∂σ W µασν , with W ν = 0 and

[µα] 1
U µαν = U ν = √ gνβ ∂σ U µαβσ , (165)
−g

where U µαβσ = (−g)(g µβ g ασ − g αβ g µσ ) is the Landau-Lifschitz superpotential. As a result of the sym-

metries of W µασν , one then has
√
16π −g(tµν + T µν ) = ∂α U µαν . (166)

By manipulating this equation, one can then define a new complex, the Landau-Lifschitz pseudotensor
µν
τLL , which is symmetric in the two indices and which satisfies
µν
16π(−g)(T µν + τLL ) = ∂α ∂β U µανβ . (167)

Since the Landau-Lifschitz superpotential is antisymmetric under exchanges of the first and second indices
and under exchanges of the third and fourth indices, one finally obtains the conservation equation
µν
∂µ [(−g)(T µν + τLL )] = 0 . (168)
 22

Evaluating the Landau-Lifschitz pseudotensor for a flat spacetime with a linear transverse traceless
perturbation (∂µ hµν µ
TT = hTT µ = 0) gives

LL 1 ρσ
ταβ = ∂α hTT
ρσ ∂β hTT , (169)
32π
which can be interpreted as encoding the stress energy of gravitational waves. The same result can be
obtained in an even simpler way by starting from the action for a linear transverse traceless perturbations
on flat space. The latter is derived by replacing gµν ≈ ηµν + hTTµν in the Einstein-Hilbert or Schrödinger
action and then expanding at quadratic order:
Z
1
d4 x ∂α hTT αβ
∂β hρσ


S=− ρσ η TT . (170)
64π
One can then apply the procedure of Eq. 160 to get to Eq. 169.
This derivation of a “stress-energy tensor” for gravitational waves, while correct, can be confusing. We
have indeed started by noticing that defining a (local) stress-energy tensor for the gravitational field is
impossible (because of the equivalence principle and the local flatness theorem) and we have ended up
deriving one for gravitational waves. The only possible explanation is that the “stress-energy tensor” for
gravitational waves is inherently a non-local object. This will become clear from the alternative derivation
that we will now undertake.
Let us start by considering a background spacetime with a small perturbation. Let us also assume
that the perturbation changes on a characteristic time and length scale λ much smaller than the back-
ground’s curvature radius L. This situation is usually referred to as geometric-optics regime and is shown
schematically in Fig. 1. Let us also define an average h. . . i over lengths and times  λ and  L.

FIG. 1: Sketch of a perturbed spacetime in the geometric-optics regime. L is the characteristic length of the background
and λ is the wavelength of the perturbation.

Let us split the spacetime metric as

B
gαβ = gαβ + hαβ + 2 jαβ + O(3 ), (171)
B
where gαβ is unperturbed background metric, while hαβ and jαβ are the first and second order pertur-
bations ( being a small parameter). Let us then consider the vacuum Einstein equations Gαβ = 0, and
expand them in  as
 
(1) (1) (2)
0 = Gαβ [g B ] + Gαβ [h, g B ] + 2 Gαβ [j, g B ] + Gαβ [h, g B ] . (172)

Here, the first term is the Einstein tensor computed with the background metric, the second one gives
the equations of motion for the first order perturbations (c.f. section II B), while the 2 term gives the
equations of motion for the second order perturbations (which are in turn comprised of terms quadratic
(2) (1)
in h, denoted by Gαβ [h, g B ], and terms linear in j, denoted by Gαβ [j, g B ]).
Taking now an average of this equation, one can note that since G(1) [. . . , g B ] is a linear operator (the
linearized Einstein tensor on the background metric), the average commutes with it, giving therefore
(1) (1)
hGαβ [h, g B ]i = Gαβ [hhi, g B ] = 0 (because the first order perturbations are oscillatory and thus average
(1) (1)
to zero, cf. section II B). Similarly, hGαβ [j, g B ]i = Gαβ [hji, g B ], where we are allowing for hji =
6 0.
In fact, the second order equations of motion have the cartoon form 2jαβ = O(h)2 , i.e. second order
 23

perturbations are sourced by products of first order ones. As such, the average of the second order
perturbations cannot be zero as hO(h)2 i will not vanish, in general. The average of Eq. 172 thus yields
 
(1) (2)
0 = Gαβ [g B ] + 2 Gαβ [hji, g B ] + hGαβ [h, g B ]i , (173)

where we have used the fact that hGαβ [g B ]i = Gαβ [g B ] since the background metric varies on scales much
larger than those on which the average is performed. This equation can then be written, by resumming
the Taylor expansion, as

Gµν [g B + 2 hji] = 8πGTµν

GW
, (174)

with
GW 1
Tµν =− hG(2) [h, g B ]i, (175)
8π µν
(2)
To compute this average, one can write Gµν [h, g B ] explicitly, use the fact that covariant derivatives
commute up to terms depending on the Riemann tensor, and show that these terms are subleading in the
geometric-optics regime λ  L. A more detailed calculation can be found in [3] and [14] and and yields
(for transverse traceless perturbations)

GW 1 ρσ
Tαβ = h∇α hTT
ρσ ∇β hTT i. (176)
32π
Remarkably, this is the same result that we derived previously, except for the average, which shows
explicitly that the stress-energy tensor of gravitational waves is a non-local object (i.e. it makes no
sense to define the stress-energy tensor of gravitational waves pointwise, but only on scales larger than
the wavelength). Just as the gravitational force disappears for an observer in a free-falling elevator if
the elevator is much smaller than the Earth (otherwise non-local tidal effects appear, cf. Eq. 159), the
gravitational wave perturbation can be set to zero by going to RNCs locally, but only on scales smaller
than the spacetime curvature radius (which is given by λ for the perturbed spacetime represented in
Fig. 1).

A. The gravitational contribution to the mass of a compact star

In this section we will take a short detour and investigate another example showing that the gravi-
tational field, although it can be set to zero locally, provides a finite contribution to the energy of the
system on non-local scales. We will consider indeed a compact spherically symmetric star, and compute
the contribution of the gravitational field to its mass.
Let us start by modeling the metric as

ds2 = −B(r)dt2 + A(r)dr2 + r2 dΩ2 , (177)

and the matter by a perfect fluid. The Einstein equations and the conservation of the fluid’s stress energy
tensor then yield (assuming asymptotic flatness) the famous Tolman-Oppenheimer-Volkoff equations [15]:
! ! !−1
dp m(r) p 4πr3 p 2m(r)
= −ρ 2 1+ 1+ 1− ,
dr r ρ m(r) r
1 dB 1 dp (178)
=− ,
2B dr p + ρ dr
1
A(r) = ,
1 − 2m(r)
r

where
Z r
m(r) = ρ(r0 )4πr02 dr0 .
0
3
2m(r)
Clearly, the first equation is the relativistic Euler equation, where ρp , 4πr p
m(r) and r are the relativistic
corrections, which disappear in the Newtonian limit c → ∞ (if c is reinstated). As can be seen, the
 24

solution to these equations reduces to the Schwarzschild metric in vacuum (and thus in the exterior
of the star). In particular, the metric in the exterior is given by the Schwarzschild metric with mass
RR
m(∞) = 0 ρ(r0 )4πr02 dr0 , where R is the radius of the star. This mass can be interpreted as the star’s
gravitational mass, as measured by an observer that were to fly satellites far from the star and interpret
their motion with Kepler’s law. More formally, one can show that this mass matches the Arnowitt-Deser-
Misner mass of the spacetime (which in turn matches its Komar mass), see e.g. [16, 17].
The total baryonic mass of the star can instead be obtained from the continuity equation for the
baryonic current j µ = mb nuµ , where mb is the average baryon mass and n is the baryon number density.
The baryonic mass is then
√
Z Z p
Mb = dx3 −gj t = mb 4πr2 A(r)n(r)dr . (179)

This mass corresponds to the sum of the rest masses of all the baryons of the system, but does not include
the internal energy. Since the internal energy density of a fluid is given by the difference ρ − mb n, we
would expect the total mass of the system to be given simply by
Z p
M? = 4πr2 A(r)ρ(r)dr . (180)

This mass, however, still differs from m(∞). To understand this discrepancy one may compute the
difference between the two, reinstate factors of G and c (using e.g. dimensional analysis) and expand it
in orders of 1/c2 . Equivalently, one can expand the difference in the weak gravity limit m(r)/r  1. By
using the fact that A(r) = [1 − 2Gm(r)/(rc2 )]−1 , one then finds

4πr2 ρ(r)Gm(r)
Z
m(∞) − M? = − dr = Uself /c2 , (181)
rc2
Z
m(r)
Uself = −G dm(r) . (182)
r

This shows (i) that the gravitational mass m(∞) is always smaller than the expected value M? , and
(ii) that in the Newtonian limit the difference is given exactly by the contribution of the (Newtonian)
gravitational self-energy, i.e. (in absolute value) the work that one would need to perform against the
gravitational force to destroy the star. Once again, we have found that even though it can be locally set
to zero, the gravitational field does contribute to the mass of an extended object. This contribution is of
the order of 10 − 20% for neutron stars.

VI. THE INSPIRAL AND MERGER OF BINARY SYSTEMS OF COMPACT OBJECTS

In this section we will use the results that we have derived to gain some semi-quantitative understanding
of the physics of binary systems of compact objects (black holes and neutron stars).
Let us first apply the quadrupole formula to a system of two compact objects with masses m1 and m2
on a circular orbit of radius r and orbital frequency Ω, which at lowest (i.e. Newtonian) order is given
by Kepler’s law
r
M
Ω= , (183)
r3
with M = m1 + m2 . The gravitational wave signal predicted by the quadrupole formula, for an observer
at distance D and angle ι with respect to the direction orthogonal to the orbital plane, then reads
 1+cos2 ι 
2 cos 2Ωt cos ι sin 2Ωt 0
hTT
ij = h
 cos ι sin 2Ωt − 1+cos2 ι cos 2Ωt 0 (184)
2
0 0 0

with the amplitude (also referred to as “strain”) being

4µΩ2 r2 4µM 2/3 Ω2/3

h= = , (185)
D D
 25

where we have also introduced the reduced mass µ = m1 m2 /M .

Several comments are worth making here. First, the frequency of the gravitational wave signal is twice
the orbital frequency. This is due to the tensor nature of gravitational waves. By reinstating G and c
(recalling that h must be dimensionless) and computing explicitly the amplitude, one finds e.g. a strain
of ∼ 10−22 for a system of two neutron stars of masses m1 = m2 = 1.4M , orbital period of 10 ms
and distance of 50 Mpc. Similarly, a strain h ∼ 10−21 can be obtained e.g. for an equal mass binary of
black holes of 30 M each, at a distance of 400 Mpc and with the same orbital period of 10 ms. Clearly,
these sources have gravitational wave frequencies ∼ 100 Hz, and as we will see they are detectable by
current ground based interferometers (LIGO-Virgo-KAGRA). Similarly, typical sources in the bands of
LISA (mHz) and pulsar-timing arrays (nHz) are e.g. a binary of massive black holes of 107 M each, with
period of a few hours and distance of few Gpc (h ∼ 10−16 ) or a binary of supermassive black holes of
109 M each, with period of one year and distance ∼ 1 Gpc (h ∼ 10−15 ), respectively. Note how small
the amplitude of these metric perturbations is compared to the amplitude of the metric perturbation at
the surface of the Sun, h ∼ GM /(R c2 ) ∼ 10−6 . Another important observation is that the strain
decays as 1/D. This is why gravitational wave observations allow for exploring the Universe up to high
redshift. Not only do gravitational waves interact very weakly with matter (since the interaction is only
gravitational), but interferometers detect directly the gravitational wave strain h, which decays more
slowly than the electromagnetic fluxes (∝ 1/d2 ) collected e.g. by optical telescopes.
Using now Eq. 185 in the stress energy tensor of gravitational waves, one can get the gravitational
wave flux TtiGW . Like all fluxes, this decays as 1/D2 , but we stress again that interferometers observe
h directly, and not the flux. Integrating the flux on a sphere far from the source, one finds that the
gravitational wave luminosity of a binary (i.e. the energy carried away by gravitational waves per unit
time) takes the simple form
!2 !2
32 G Gm1 m2 v
ĖGW = , (186)
5 c3 r2 c

where we have explicitly reinstated G and c. The energy that gravitational waves remove from the source
must of course come from the system’s kinetic and potential energy. For a circular Keplerian binary, the
sum of kinetic and potential energy is simply given, in the center of mass frame, by
1 2 GM µ 1 GM µ
Etot = µv − =− . (187)
2 r 2 r
By requiring energy conservation (Ėtot = −ĖGW ), one can then obtain an expression of the rate of change
of the separation, ṙ, due to gravitational wave emission. Clearly ṙ < 0, i.e. the binary slowly spirals in
under the backreaction of gravitational waves. This can of course be interpreted as a PN contribution to
the acceleration of the system, i.e.
 1 1  1 
~a12 = ~aN 1 + O 2 + O 4 + O 5 , (188)
c c c
where ~aN is the Newtonian acceleration, O( c12 ) is the 1PN (conservative) correction, O( c14 ) is the 2PN
(conservative) correction and O( c15 ) is the 2.5PN (dissipative) backreaction of gravitational waves. Note
that the O( c15 ) scaling of the last term follows from the factors of 1/c in Eq. 186.
Using Kepler’s law, ṙ can be recast into the rate of change of the orbital angular frequency, Ω. Using
then the relation between gravitational wave frequency and Ω, f = 2forb = Ω/π, one finally obtains
96 1 8/3
f˙ = π (GMc )5/3 f 11/3 , (189)
5 c5
where we have introduced the chirp mass

Mc = M η 3/5 , (190)
µ
where η = M is the symmetric mass ratio. The chirp mass is indeed the quantity that can be most easily
estimated from the gravitational wave signal from inspiraling binaries: as gravitational waves remove
energy and angular momentum, the binary spirals in, the separation decreases, and the frequency of
gravitational waves increases depending on Mc alone (at leading PN order). This expression can also be
recast into an equation for the rate of change of the orbital period. The latter is the quantity monitored in
binary pulsar systems, i.e. systems at least one component of which is a millisecond pulsar. The presence
 26

of the pulsar allows for tracking the period of the binary system with exquisite accuracy, historically
providing for the first time evidence for the existence of gravitational waves [18].
In order to gain more qualitative understanding of the inspiral phase beyond the leading PN order, we
will now make a short detour and recall the most salient features of geodesic orbits in Schwarzschild and
Kerr spacetimes. While this is only applicable to binaries with very small mass ratio q = m2 /m1  1,
the qualitative features that we will discover (e.g. the effect of spins, the plunge) will survive even at
mass ratios q ≈ 1. This is somewhat expected from Newtonian mechanics, where one can map a binary
with arbitrary masses into a particle with the reduced mass µ around a particle with the total mass
M , but it is not at all obvious in GR. Only recently has evidence started accumulating that a similar
mapping between arbitrary binaries and the test-particle limit may exist even in PN theory, although
approximately. This approximate mapping goes under the name of ‘effective-one-body’ model [19].

A. Geodesics in Schwarzschild and Kerr

Let us start by studying geodesics in the Schwarzschild spacetime, whose line element we write in areal
coordinates in the usual form
dr2
 
2 M
ds = − 1 − 2 dt2 + + r2 dΩ2 . (191)
r 1 − 2M/r

Since the metric is static and spherically symmetric, we can look at equatorial geodesics without loss of
generality (i.e. the coordinates can always be chosen to be such that orbits have θ = π2 ).
Let us start with particles having non-zero mass (timelike geodesics). From the existence of the two
Killing vectors ∂t and ∂φ , it follows that the specific 6 energy and angular momentum observed at infinity,
i.e. E = −ut and L = uφ with uµ = dxµ /dτ the four-velocity, are conserved. One can then obtain the
first-order equations
 
2M dt
1− = E,
r dτ
(192)
dφ
r2 = L.
dτ
Moreover, by using these equations in the conservation of the norm uµ uµ = −1, one can obtain a first-
order equation for the radial motion:
 2
1 dr 1 2
+ V (r) = E , (193)
2 dτ 2

with the effective potential being given by

1 M L2 M L2
V (r) = − + 2− 3 . (194)
2 r 2r r
Apart from the first (constant and thus irrelevant) term, this potential includes the Newtonian potential,
the usual Newtonian centrifugal term, but the last term does not appear in Newtonian mechanics. In
fact, it is a 1PN term, as can be seen by reinstating c by dimensional analysis.
As a consequence of this term, the behavior of the potential at small r drastically differs from the
Newtonian one. Unlike the latter, which predicts the existence of stable circular orbits down to arbitrarily
small radii, Eq. 194 predicts the existence of an innermost stable circular orbit (ISCO) at r = 6M . This
can be seen by solving the equations defining circular orbits, E 2 /2 − V (r) = V 0 (r) = 0, and checking
the sign of V 00 (r) to assess stability. As can also be understood by plotting V , for E > 1 only unstable
circular orbits exist. However, for each value of E < 1 (corresponding to bound orbits) two circular orbits
exist, with√ the one at larger radius being
√ stable and the other being unstable. These orbits exist √ only
for L ≥ 2 3M and coincide for L = 2 3M , which corresponds to the ISCO (r = 6M and E = 2 2/3).
The unstable circular orbits lie instead at r < 6M , but they always have r > 3M (they only approach
r = 3M in the ultra-relativistic limit E, L → ∞).

6 By “specific”, we mean “normalized by the particle’s mass”.

 27

By redoing the same analysis for null orbits, one can prove that the radial motion obeys
 2
b2
 
1 dr 2M
= −V ph = 1 − 1 − , (195)
E 2 dλ r2 r

where λ is an affine parameter and b = L/E is usually referred to as impact parameter. An analysis √ of the
effective potential Vph shows that a circular orbit exists only for a critical value of b, namely b = 3 3M .
This circular null orbit (also known as light ring) lies at a radius of r = 3M and is unstable. Its existence
and properties are not only important for the interpretation of the observations by the Event Horizon
Telescope (EHT) [20], but also for the physics of black hole quasinormal modes, as we will see in the
following.
The case of geodesics in a Kerr spacetime is slightly more involved, because one can no longer assume
equatorial motion due to the absence of spherical symmetry. In Boyer-Lindquist coordinates the metric
reads
 
2 2M r Σ
ds = − 1 − dt2 + dr2 + Σ dθ2
Σ ∆
2M a2 r
 
4M ar
2
+ r +a + 2
sin θ sin2 θ dφ2 −
2
sin2 θ dt dφ, (196)
Σ Σ

where a = S/M (with S the spin), M is the mass and

Σ = r2 + a2 cos2 θ, ∆ = r2 − 2M r + a2 . (197)

The Killing vectors ∂t and ∂φ still exist, and imply that the specific energy E = −ut and angular
momentum (in the spin direction) L = uφ must be conserved. However, one more conservation equation
(besides the unit norm condition) is needed to reduce the equations of motion to first order. Fortunately,
the Kerr geometry has a ‘hidden’ symmetry, which can be described by a Killing-Yano tensor. This
symmetry implies the existence of an additional constant of motion, the Carter constant Q, which allows
for writing the equations for timelike geodesics (with particle mass µ) as
 2
dr dt
= Vr (r), = Vt (r, θ),
dλ dλ
 2
dθ dφ
= Vθ (θ), = Vφ (r, θ) , (198)
dλ dλ
with
$4 $2
   
2 2
Vt (r, θ) ≡ E − a sin θ + aL 1 − , (199a)
∆ ∆

2
Vr (r) ≡ E$2 − aL − ∆ r2 + (L − aE)2 + Q ,
 
(199b)
2 2 2 2 2
Vθ (θ) ≡ Q − L cot θ − a (1 − E ) cos θ, (199c)
 2
a2 L

$
Vφ (r, θ) ≡ L csc2 θ + aE −1 − , (199d)
∆ ∆
where we have defined

$2 ≡ r2 + a2 (200)

and the “Carter time” λ by

dτ
≡Σ. (201)
dλ
The most striking difference from the Schwarzschild case is the presence of a non-trivial equation for
the angular variable θ, which describes the precession of the orbital angular momentum of the particle
around the spin of the Kerr geometry. Moreover, the spin parameter a also enters the effective potential
for the radial motion. This has profound implications for the particle’s motion. Focusing for instance
on non-precessing (i.e. equatorial) orbits, we can set θ = π/2 and Q = 0, and we can obtain the radius,
specific energy and specific angular momentum of circular orbits by solving Vr (r) = Vr0 (r) = 0. Like in
 28

9
4.0
8
7 3.5
6 3.0

LISCO/M
rISCO/M

5
2.5
4
3 2.0
2 1.5
1
1.0
1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00

0.95
0.90
0.85
0.80
EISCO

0.75
0.70
0.65
0.60
1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00

FIG. 2: ISCO radius, specific angular momentum and specific energy for equatorial orbits in Kerr.

the Schwarzschild case, this analysis shows the existence of an ISCO, whose properties depend on the
spin as [21]
s
2
EISCO (χ) = 1 − , (202)
3rISCO (χ)
2 h p i
LISCO (χ) = √ 1 + 2 3rISCO (χ) − 2 , (203)
3 3

χp
rISCO (χ) = 3 + Z2 − (3 − Z1 )(3 + Z1 + 2Z2 ) , (204)
|χ|
h i
Z1 = 1 + (1 − χ2 )1/3 (1 + χ)1/3 + (1 − χ)1/3 , (205)
q
Z2 = 3χ2 + Z12 , (206)

where χ = a/M . The parameter χ ranges from -1 to 1, with positive values corresponding to orbits
co-rotating with the Kerr black hole, while negative values correspond to counter-rotating orbits. In the
limit χ → 0, these expressions reduce to those for the Schwarzschild ISCO.
The dependence of these expressions on the spin is a special case of a more general phenomenon
appearing in GR, the dragging of inertial frames or Lense-Thirring effect. Unlike what happens in
Newtonian theory, the spin has a clear impact on the dynamics, “dragging” matter into rotation. In
accordance with this, Eqs. 202–206 predict that as the spin increases in magnitude, prograde orbits (i.e.
ones co-rotating with the Kerr black hole) have smaller and smaller ISCO radii, down to rISCO = M in
the extremal limit χ = 1. 7 As a result, the values of EISCO and LISCO decrease as χ grows from 0 to 1.

7 Note that although the ISCO seems to coincide with the event horizon in this limit, this is simply an artifact of the
Boyer-Lindquist coordinates becoming singular in the extreme limit, as can be seen by computing the proper distance
between the ISCO and the event horizon [21].
 29

4.0 4.0

3.5 3.5

3.0 3.0

phM
rph/M

2.5 2.5

2.0 2.0

1.5 1.5

1.0 1.0
1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00

FIG. 3: Circular photon orbit’s radius and orbital frequency for equatorial orbits in Kerr.

For retrograde orbits (χ < 0) the ISCO instead moves to larger radii as the spin magnitude increases, up
to r = 9M in the extreme limit. As a consequence, the values of EISCO and LISCO increase as χ goes from
0 to -1. These behaviors are shown in Fig. 2.
For a photon in the Kerr geometry, the geodesics equations are instead [21]
 2
dr dt
= Vr (r), = Vt (r, θ),
dλ̃ dλ
e
 2
dθ dφ
= Vθ (θ), = Vφ (r, θ) , (207)
dλ
e dλ
e

e a “time” parameter, b = L/E, q = Q/E 2 , and

with λ

−a2 ∆ sin2 θ + ab(∆ − $2 ) + $4

Vet (r, θ) ≡ (208a)
∆
Vr (r) ≡ ($ − ab) − ∆ (a − b)2 + q ,
2 2
 
e (208b)
Veθ (θ) ≡ a2 cos2 θ − b2 cot2 θ + q, (208c)
2
b a(ab + ∆ − $ )
Veφ (r, θ) ≡ − . (208d)
sin2 θ ∆
Specializing again to equatorial orbits (q = 0, θ = π/2), one finds, like in Schwarzschild, that there exists
an unstable circular photon orbit, whose coordinate radius reads [21]

rph (χ) = 2M {1 + cos[2/3 arccos(−χ)]} . (209)

This is plotted in Fig. 3, together with the orbital frequency Ωph = dφ/dt (obtained from Eq. 208).
As can be seen, the frame dragging once again makes the photon orbit radius decrease with χ, with
rph → M as χ → 1. 8 Like in the Schwarzschild case, not only are circular photon orbits relevant for
EHT observations, but also for the physics of quasinormal modes, as we will see below.

B. A qualitative description of the inspiral and merger

Let us now utilize what we have learned so far to gain some qualitative understanding of the evolution of
a quasi-circular binary of compact objects (black holes or neutron stars). As have seen, gravitational waves
are emitted and carry energy and angular momentum away from the system. As a result, the separation
of the binary decreases and the orbital frequency increases. The rates of change of these quantities, as
we have seen, only depend on a combination of the two masses, the chirp mass, at leading PN order.

8 As for the ISCO, this is due to the coordinates becoming singular in the extreme limit. The proper distance between the
circular photon orbit and the event horizon remains non-zero [21].
 30

FIG. 4: Two waveform approximants for the signal h+ (normalized by the luminosity distance DL and the system’s total
mass M ) produced by a black hole binary with mass ratio 1 : 5, spin parameters 0.5 and 0 (for the heavier and lighter
object respectively, and with the non-vanishing spin initially in the orbital plane) as function of retarded time. The
observer is located 60 degrees away from the orbital angular momentum axis. Taken from Ref. [22]

However, at higher PN orders the gravitational wave fluxes also depend on the individual masses and
spins. PN corrections also appear in the conservative sector, changing e.g. the relation between the
orbital frequency and the system parameters (which is only given by Kepler’s law at Newtonian order,
cf. sections III and VI A), giving rise to precession (if at least one spin is non-zero and misaligned with
the orbital angular momentum, cf. section VI A), etc. Precession (spin-spin and spin-orbit) introduces
modulations in the gravitational waveforms (both in amplitude and phase), as shown for instance in
Fig. 4. This makes measurements of the spin directions possible (at least in principle) with gravitational
wave detectors.
As the binary’s separation shrinks, the system transitions from one circular orbit to the next until it
either reaches the ISCO or the two bodies touch. We have seen that an ISCO exists in the test-particle
limit (i.e. for geodesics), but a similar transition to unstable circular orbits occurs also for comparable-
mass binaries of black holes [23] (neutron stars touch and interact before they reach this effective ISCO).
When the separation reaches the effective ISCO, or when the two bodies touch (in the case of neutron
stars), the binary plunges and merges. The merger phase can only be studied via numerical-relativity
simulations, but at least for black holes the post-merger phase can be understood analytically in terms
of quasi-normal modes, as we will see in the next section. As for neutron stars, the post-merger phase
depends critically on the microphysics (e.g. on the equation of state of nuclear matter) and can only be
predicted via numerical simulations.
The position of the effective ISCO, just like in the test-particle limit, depends critically on the spins
of the two objects. The larger the spins, the more the effective ISCO moves inwards and the longer the
binary emits gravitational waves before plunging. This is exactly the same effect that takes place, in the
electromagnetic case, for geometrically thin, optically thick accretion disks. In the latter, the gas spirals
in on quasi-circular orbits, as it loses energy because of friction with the neighboring gas elements. The
gas potential and kinetic energy is thus converted into heat, and eventually radiated away (e.g. in the
optical, infrared, UV and/or X-ray bands). This process can only continue, however, until the gas reaches
the ISCO of the central black onto which it is accreting. For a gas element with unit mass starting at rest
at infinity, the energy when it reaches the ISCO is EISCO (as given by Eq. 202). By energy conservation,
the radiative efficiency of an accretion disk is therefore η = 1 − EISCO , which is a strong function of spin,
as can be seen from Fig. 5 (top panel).
Computing a similar efficiency for gravitational waves from binary systems is not easy away from
the test-particle limit, as the effective ISCO energy is not known analytically for comparable masses.
However, one can perform numerical relativity simulations, which show in fact exactly the same effect
(qualitatively). Fig. 6, taken from Ref. [26], shows the “trajectories” (in some sense) of two black holes
with spins respectively aligned (right) and antialigned (left) with the orbital angular momentum. It can be
visually seen that the orbit with aligned spins (corresponding to prograde Kerr geodesics) reaches smaller
separations and performs more cycles than that for antialigned spins (corresponding to retrograde Kerr
geodesics). This effect is known in the literature as “orbital hang up”, but it is really a manifestation
of the frame dragging of GR. The same effect can be seen at play in Fig. 5 (bottom panel), adapted
 31

0.40
0.35
0.30
0.25
0.20
0.15
0.10
0.05
1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00

FIG. 5: Top: The radiative efficiency of a geometrically thin, optically thick accretion disk around a Kerr black hole with
spin parameter χ. Note that the maximum spin expected for black holes surrounded by such accretion disks is
χ = 0.998 [24], for which the efficiency is ≈ 32%. Bottom: the fraction of the total mass M = m1 + m2 emitted by a black
hole binary system in gravitational waves, as a function of the average of the projections of the two spins on the orbital
angular momentum axis, (χ1 + χ2 )/2. (Adapted from Ref. [25].)

from Ref. [25], which collected the energy emitted in gravitational waves in various black hole binaries
simulated in the literature, and plotted it as function of a combination of the two spins (projected on
the orbital angular momentum axis). Although the emission efficiency tops at about 10% at high spins,
thus remaining lower than the electromagnetic efficiency shown in the top panel of Fig. 5, the behavior
is qualitatively the same as the latter.
As a final application, let us show that the knowledge that we have gained so far, albeit qualitative,
allows for interpreting the data of GW150914 [1], the first direct gravitational wave detection, and for
concluding that the components of this binary system must be black holes. A spectrogram of this event is
 32

FIG. 6: Numerical relativity simulations of binary black holes with spins |χ1 | = |χ2 | = 0.76 either antialigned (left) or
aligned (right) with the orbital angular momentum. The squares and circles represent the puncture positions every 10 M
(with M = m1 + m2 ) of evolution, and the lines joining them can therefore be thought as the trajectories of the black
holes. The dashed green circles represent the first common apparent horizon.

FIG. 7: A spectrogram of the GW150914 event observed by LIGO Hanford (left) and Livingston (right). The color code
represents the normalized strain amplitude. Figure adapted from Ref. [1].

shown in Fig. 7: the color code represents the gravitational wave strain amplitude in a given bin of time
(x-axis) and frequency (y-axis). One can clearly recognize a “chirp”, i.e. an increase of the gravitational
wave frequency with time, which one can fit with Eq. 189 to obtain Mc ≈ 30M . This in turn implies,
through Eq. 190, that the total mass must be M & 70M . The power’s peak, which one expects to
coincide with the plunge/merger phase, lies at f ∼ 150 Hz. Translating that into an orbital frequency
(dividing by a factor 2), and using Kepler’s law to convert to a separation (assuming for simplicity roughly
equal masses), one finds that the plunge/merger takes place when the binary separation decreases to just
350 km. This is very close to the sum of the Schwarzschild radii of the two objects, GM/c2 & 210 km.
Therefore, the separation at which the plunge happens is comparable to the effective ISCO radius, which
seems to favor the hypothesis that the two objects are black holes. In fact, if the objects were stars, they
would touch, interact and plunge way before reaching the effective ISCO separation, i.e. the two objects
must be very compact. Among compact objects – i.e. ones with Gm/(Rc2 ) = O(1), with m and R the
object’s mass and radius – the only options (in GR) are black holes and neutron stars. Neutron stars,
however, are excluded because they cannot be more massive than 2 − 3M . This leaves black holes as
 33

FIG. 8: Top: the GW150914 data, as observed by the Hanford and Livingstone LIGO detectors. Bottom: the de-noised
signal, reconstructed with template and wavelet techniques, alongside the prediction from numerical relativity simulations.
Figure adapted from Ref. [1]

the only possibility.

VII. THE POST-MERGER SIGNAL

As can be seen even in real strain data (c.f. Fig. 8 for the GW150914 event), after the amplitude peaks
at the merger, the gravitational wave signal from a black hole binary seems to be well described by one (or
more) damped sinusoids. This is in fact what happens, as can be understood rather easily by employing
linear perturbation theory on a Schwarzschild or Kerr background, as we will do in this section. We will
start with the simple toy problem of a test Klein-Gordon field on a Schwarzschild background, and we
will then move to the case of gravitational perturbations on the same geometry. We will finally generalize
the treatment to the Kerr case, which will allow for concluding that indeed the post-merger signal is well
described by a linear superposition of quasi-normal modes of the final (spinning) black hole resulting
from the merger. For more details, we refer the reader to the extensive review [27].

A. Scalar perturbations of non-spinning black holes

Let us start by considering the toy problem of a free scalar field on a Schwarzschild background, i.e.
let us consider the Klein-Gordon equation

g µν ∇µ ∇ν ϕ = 0, (210)

with g µν the contravariant components of the Schwarzschild metric given by Eq. 191. Since the metric is
static and spherically symmetric, it is natural to decompose the scalar field in spherical harmonics (Y`m )
and Fourier modes (e−iωt ) as
X R`m (r)
ϕ= Y`m (θ, φ)e−iωt , (211)
r
`,m
 34

where R`m characterizes the radial profile of the scalar field. By replacing this ansatz in the Klein-Gordon
equation, that reduces to a single ordinary differential equation in the radial coordinate,

d2 R`m
+ ω 2 − V` R`m = 0,


2
(212)
dr∗
where
 r 
r∗ ≡ r + 2M ln −1 , (213)
2M
is the tortoise coordinate and the potential V` the potential is given by
  
2M `(` + 1) 2M
V` (r) ≡ 1 − + 3 . (214)
r r2 r

As can be seen, in the geometric optics (eikonal) limit `  1 this potential reduces to the effective
potential for the radial motion of photons, Eq. 195, if one identifies the impact parameter b = L/E of
photons with `. This is expected, since in the eikonal limit the wavefronts of a scalar field satisfying
the Klein-Gordon equation move along null geodesics of the metric. We also stress that going from the
partial differential equation 210 to the single ordinary differential equation 212 is highly non-trivial, and
depends critically on the choice of the ansatz of Eq. 211.
To solve Eq. 212, we need to impose suitable boundary conditions. As r∗ → ∞, in order for
nothing to enter the system, we need to impose outgoing boundary conditions R`m ∼ exp (iωr∗ ).
Conversely, since nothing can escape the horizon (which corresponds to r∗ → −∞), we need to impose
R`m ∼ exp (−iωr∗ ) there (ingoing boundary conditions). Solving this boundary-value problem, one
obtains a discrete spectrum of complex frequencies ω. The corresponding excitations are referred to
as quasinormal modes (as opposed to normal modes, which have real frequencies). Moreover, one can
check that all frequencies in the spectrum have negative imaginary part, which corresponds to damped
modes (c.f. Eq. 211). This shows that a test scalar field is (linearly) stable on the Schwarzschild geometry.

Exercise 4: Plot the effective potential of Eq. 214 and approximate it qualitatively with a rect-
angular potential. Solve Eq. 212 in the three regions in which this rectangular potential is constant, and
impose appropriate junction conditions at the transition radii and ingoing/outgoing boundary conditions
at the horizon and at infinity. By counting the integration constants, show that the spectrum is discrete
and complex. Solve numerically for a few frequencies in the spectrum.

B. Tensor perturbations of non-spinning black holes

A similar analysis can be performed for the metric perturbation hµν of a Schwarschild spacetime.
To exploit again the fact that the Schwarzschild geometry is static and spherically symmetric, we can
decompose the time dependence in Fourier modes and the angular dependence in scalar, vector and tensor
harmonics. In more detail, htt htr , hrr are scalars on the two-sphere (i.e. under rotations), and can thus
be expanded in the usual (scalar) spherical harmonics Y`m . The cross terms htA and hrA (with capital
Latin letters spanning the two angles θ, φ) are instead vectors on the two sphere. A basis for vectors on
the two-sphere can be obtained by taking gradients of the scalar harmonics,

YAE,`m = ∂A Y`m , (215)

or exterior derivatives of these gradients,

YAB,`m = A B ∂B Y`m , (216)

where AB is the Levi-Civita tensor on the two-sphere and the angular indices are raised and lowered
with the metric of the two-sphere, γAB :

γAB dX A dX B = dθ2 + sin2 θ dφ2 , (217)

√
AB = γ eAB , (218)

with γ = det(γAB ) = sin θ and eθφ = −eφθ = 1, eθθ = eφφ = 0. Because scalar harmonics have parity
(−1)` , YAE and YAB have respectively parity (−1)` and (−1)`+1 . Similarly, the components hAB transform
 35

as a tensor on the two-sphere, and can thus be decomposed in the following basis

E,`m `m 1 ;C
YAB = Y;A;B − γAB Y`m;C ,
2
T ,`m
YAB = Y`m γAB , (219)
B,`m B,`m 1  C E,`m E,`m

YAB = Y(A;B) = A YCB + B C YCA ,
2
where the semicolon denotes covariant derivatives on the two-sphere. Because of the parity of the scalar
E T B
harmonics, YAB and YAB have parity (−1)` , while YAB has parity (−1)`+1 . Note also that by using the
explicit expressions for the scalar spherical harmonics, the vector harmonics YAE,`m and YAB,`m start at
` = 1 (i.e. YAE,00 = YAB,00 = 0), while YAB
E,`m B,`m
and YAB E,`m
start at ` = 2 (i.e. YAB B,`m
= YAB = 0 for ` ≤ 1).
Expanding the metric perturbation in these scalar, vector and tensor harmonics one can then set
X 2M

htt = 1− H0 (r)Y`m (θ, φ)e−iωt (220)
r
`m
X
htr = H1 (r)Y`m (θ, φ)e−iωt (221)
`m
X −1
2M
hrr = 1− H2 (r)Y`m (θ, φ)e−iωt (222)
r
`m

[−h0 (r)YAB,`m (θ, φ) + H0 (r)YAE,`m (θ, φ)]e−iωt

X
htA = (223)
`m

[−h1 (r)YAB,`m (θ, φ) + H1 (r)YAE,`m (θ, φ)]e−iωt

X
hrA = (224)
`m
T ,`m E,`m B,`m
X
hAB = [r2 K(r)YAB (θ, φ) + r2 G(r)YAB (θ, φ) + h2 (r)YAB (θ, φ)]e−iωt , (225)
`m

where H0 , H1 , H2 , h0 , h1 , H0 , H1 , K, G, h2 are free radial functions.

By similarly expanding the generator ξ µ of gauge transformations in scalar (ξ t and ξ r ) and vector (ξ A )
harmonics, one can set h2 = H0 = H1 = G = 0 (Regge-Wheeler gauge [28]). Moreover, without loss of
generality one can set m = 0 when studying perturbations of Schwarzschild. In fact, because of spherical
symmetry the modes with m 6= 0 can be set to zero by rotation φ → φ + const. Note that is also true
for the scalar field considered above (indeed m does not appear in the potential of Eq. 214). In the odd
sector the metric perturbation then becomes

0 0 0 h0 (r) 
 

 0 0 0 h1 (r)  ∂
hodd
µν = sin θ Yl0 (θ)e−iωt , (226)
0 0 0 0  ∂θ
h0 (r) h1 (r) 0 0

whereas in the even sector one has

H0 (r) 1 − 2M


r H1 (r) 0 0
2M −1


heven
 H1 (r) H2 (r) 1 − r 0 0  Yl0 (θ)e−iωt

µν =  (227)

0 0 r2 K(r) 0 
0 0 0 r2 K(r) sin2 θ

with xµ = (t, r, θ, φ). Because of their different parity, the even and odd parity perturbations decouple
when this ansatz is replaced into the Einstein equations. By using the latter, in the odd sector one obtains
the famous Regge-Wheeler equation [28]

d2 Ψ
+ ω2 − V Ψ = 0 .


(228)
dr∗2

with
  
− 2M `(` + 1) 6M
V =V = 1− − 3 (229)
r r2 r
 36

and
 
h1 (r)
− 2M i d
rΨ− .


Ψ=Ψ = 1− , h0 = (230)
r r ω dr∗

Similarly, the even sector is described by the same Eq. 228, but with the Zerilli [29] potential

2M 3λ2 M r2 + λ2 (1 + λ) r3 + 9M 2 (M + λr)
 
+ 2
V =V = 3 1− 2 , (231)
r r (3M + λr)

where λ ≡ (` − 1)(` + 2)/2, and the Zerilli variable Ψ = Ψ+ is defined implicitly by

6M 2 + λ (1 + λ) r2 + 3M λr + dΨ+
K = Ψ + , (232)
r2 (3M + λr) dr∗


iω 3M 2 + 3λM r − λr2 + iωr dΨ+
H1 = Ψ − . (233)
r (3M + λr) (1 − 2M/r) 1 − 2M/r dr∗

Note that H0 can also be obtained from the algebraic relation

   
6M `(` + 1)
(` − 1)(` + 2) + H0 + i M − 2iω r H1
r ω r2
2ω 2 r2 + 2M 2 /r2
 
2M
− (` − 1)(` + 2) + − K = 0, (234)
r 1 − 2M/r

which also follows from the Einstein equations.

A few comments are in order at this stage. First, even though the Regge-Wheeler and Zerilli potentials
are different, they give rise to the same spectrum of quasinormal mode frequencies when ingoing boundary
conditions are imposed at the horizon and outgoing boundary conditions are imposed at infinity, i.e. the
two potentials are isospectral. The resulting frequencies have negative imaginary part, which is consistent
with the Schwarzschild geometry being linearly stable to gravitational perturbations. Moreover, both the
Regge-Wheeler and Zerilli potentials coincide, in the eikonal limit `  1, with the scalar potential of
Eq. 214 and the photon potential of Eq. 195. As discussed above, this is expected and simply amounts to
the fact that gravitational wavefronts travel along null geodesics, but it is also of practical importance.
In fact, in the eikonal limit the peak of the quasinormal mode potential must coincide with the peak of
the photon potential, which lies at the location of the (unstable) circular photon orbit. The real part of
the quasinormal mode frequencies can then be shown to be simply proportional to the orbital frequency
of the circular photon orbit, while the imaginary part turns out to be related to the Lyapunov exponent
of null geodesics near the circular photon orbit, which in turn depends on the curvature of the photon
effective potential near its peak [30, 31]. Intuitively, this means that quasinormal modes can be thought
of as being generated at the circular photon orbit, after which they slowly leak outwards (because the
circular photon orbit is unstable to radial perturbations).

C. Tensor perturbations of spinning black holes

The calculation of the quasinormal modes of spinning black holes is considerably more complicated.
In fact, already for a test scalar field in the Kerr geometry it is not obvious at all that the Klein-Gordon
equation can be reduced to a one-dimensional Schrödinger-like equation. This is definitely not the case
if the scalar is decomposed in spherical harmonics (as in Eq. 211). Still, the existence of a Killing-Yano
tensor for the Kerr geometry, which allows for separating the equations for geodesic motion (which also
regulate the motion of scalar and tensor wavefronts in the eikonal limit), suggests that the equations for
both scalar and tensor perturbations may similarly separate under an appropriate choice of basis on the
two-sphere.
In fact, one can try to solve the Klein-Gordon equation on the Kerr metric in Boyer-Lindquist coordi-
nates (Eq. 196) by separation of variables, i.e.

ϕ = R(r)Θ(θ)eimφ e−iωt , (235)

 37

which yields [32]

  h
∂ ∂R
2 i
+ a2 m2 − 4M ramω + r2 + a2 ω 2 R = Q + m2 + ω 2 a2 ∆R,


∆ ∆ (236)
∂r ∂r
m2
   
1 ∂ ∂Θ
+ a2 ω 2 cos2 θ − Θ = − Q + m2 Θ ,


sin θ 2 (237)
sin θ ∂θ ∂θ sin θ
where Q is a separation constant (defined to reduce to the Carter constant in the eikonal limit). As can
be seen, for a = 0 the second equation reduces to the equation defining associated Legendre polynomials,
i.e. for a = 0, Θ(θ)eimφ reduces to spherical harmonics. For a 6= 0, Eq. 237 defines instead the scalar
spheroidal harmonics.
A similar calculation is possible, albeit much more involved, for metric perturbations in Kerr. We will
provide no proof here, but just state the result. Let us first introduce the Newman-Penrose scalars

Ψ0 = −Cµνλσ lµ mν lλ mσ , (238)
µ ∗ν λ ∗σ
Ψ4 = −Cµνλσ n m n m , (239)

with Cµνλσ the Weyl curvature tensor, and l , n , m , m∗ a (complex) null tetrad defined at each spacetime
point. Note that Ψ0 and Ψ4 can be thought of as describing ingoing and outgoing gravitational wave
signals. If one defines the tensor spheroidal harmonics s Slm [33]

(m + su)2
   
∂ 2 ∂ 2 2 2
(1 − u ) s Slm + a ω u − 2aωsu + s + s Alm − s Slm = 0 , (240)
∂u ∂u 1 − u2

with s = ±2 and u = cos θ, one can then decompose

Z ∞ X
l
1 X
ψ(t , r , θ , φ) = e−iωt eimφ s Slm (θ)Rlm (r)dω , (241)
2π
l=|s| m=−l

where ψ stands for either Ψ0 (in which case s = 2) or ρ−4 Ψ4 [with ρ ≡ −1/(r − ia cos θ) and s = −2].
Here, A`m is a separation constant, which for a = 0 can be computed analytically to be A`m (a = 0) =
`(` + 1) − s(s + 1). This ansatz allows for solving the linearized Einstein equations by separation of
variables, and lead to the master radial equation [33]

∆∂r2 Rlm + 2(s + 1)(r − M )∂r Rlm + V Rlm = 0 . (242)

with
 
2 2 1 2 2 2 2 2 2 2 2


V = 2isωr−a ω − s Alm + (r +a ) ω −4M amωr+a m +2is am(r − M ) − M ω(r − a ) . (243)
∆

Again, this equation can be solved as a boundary value problem with ingoing boundary conditions at the
event horizon and outgoing boundary conditions at infinity. This yields a discrete spectrum of complex
quasinormal mode frequencies, whose imaginary part is again negative for |a| = 6 M , pointing at linear
stability. In the extreme limit a → ±M , the imaginary part goes to zero for some modes, which suggests
that extreme Kerr is only marginally stable.
As can be seen, the master equation defining the quasinormal mode spectrum only depends on the
mass and spin of the black hole, M and a. This is a manifestation of the no-hair theorem of GR, and
can be used to perform consistency tests of the latter [34]. Indeed, if two modes are observed, one can
use the real and imaginary part of the first to obtain M and a, and then use them to predict the real
and imaginary part of the second mode. These predictions can then be compared to the measurements.

VIII. THE DETECTION OF GRAVITATIONAL WAVES

In this section, we will review the principles behind the experimental detection of gravitational waves.
The effort to observe these signals started in the 60s with Weber [35], who pioneered the use of resonant
bars. The idea behind this setup is that gravitational waves change the size of the bar (or any object
for that matter) in a periodic fashion (with frequency given by the wave’s frequency f ). If the latter
coincides with one of the characteristic frequencies of the bar, the system can resonate and the detector’s
 38

FIG. 9: Schematic representation of a Michelson (left) and a Fabry-Perot (right) interferometer (courtesy
Caltech/MIT/LIGO Laboratory). In both, the 45-degrees slab represents a beam splitter (which reflects 50% of the
photons and lets 50% through) and the photodetector is represented by a dot. The two mirrors close to the beam splitter
in the Fabry-Perot design can be thought of as semi-transparent, e.g. a photon will go back and forth many times along
each arm before reaching the photodetector. This “Fabry-Perot cavity” thus increases the effective length of the arms (by
a factor ∼ 300 for LIGO). Note that e.g. LISA will instead be a Michelson interferometer.

response is therefore amplified. A similar idea is behind the more modern interferometric detectors, in
which lasers travel back and forth between mirrors in two or more arms, before being collected at the same
point (a photodetector; c.f. schematic setup in Fig. 9). Because gravitational waves have a non-trivial
angular dependence, the lengths of the two arms change in different fashions. This time-varying length
difference gives rise to a phase difference between the lasers, which can be measured by observing their
interference pattern at the photodetector. Interferometric detectors include for instance the ground-based
LIGO, VIRGO and KAGRA experiments; the next generation ground-based detectors Cosmic Explorer
and Einstein Telescope; and the future space-borne interferometer LISA. In the following, we will derive
in detail the response of an interferometer to a gravitational signal.

A. The response of a gravitational wave detector: the low frequency limit

Let us consider a laser traveling back and forth between two mirrors in free fall. 9 More precisely, let
us assume that the mirrors move along geodesics in a flat spacetime perturbed by a gravitational wave
hTT
µν :

gµν = ηµν + hTT

µν ,
d2 xµ dxα dxβ (244)
2
+ Γµαβ = 0.
dτ dτ dτ
Changing variable from the proper time τ to the coordinate time t, one can compute the mirror’s coor-
dinate acceleration as
! "  −1 # "  −1  −2 µ 2 #
d2 xµ d dxµ dτ d dxµ dt dτ d2 xµ dt dt dx d t dτ
= = = − .
dt2 dt dτ dt dτ dτ dτ dt dτ 2 dτ dτ dτ dτ 2 dt

d2 xµ d2 t
Expressing dτ 2 and dτ 2 by using the geodesic equations, one obtains

d2 xµ
= −Γµαβ ẋα ẋβ + ẋµ Γ0αβ ẋα ẋβ .
dt2

9 For ground detectors such as LIGO and Virgo, the mirrors cannot be thought of as strictly speaking in free fall since
the experiment takes place on Earth. However, the mirrors are isolated from the Earth’s motion and vibrations by a
sophisticated suspension system, and can thus be thought of as effectively in free fall, at least at high frequencies.
 39

where the dot denotes d/dt. Since ẋt = 1, the spatial acceleration is

d2 xi
= − Γitt + 2Γitj v j + Γijk v j v k + v i Γttt + 2Γttj v j + Γtjk v j v k ,



2
(245)
dt
where v i ≡ ẋi . Assuming that the mirror moves at v  1, one finally obtains

d2 xi
≈ −Γitt . (246)
dt2
It is now easy to see that
1 ij
Γitt = δ 2∂t hTT TT


jt − ∂j htt = 0, (247)
2
for a gravitational perturbation in the transverse traceless gauge (c.f. section II A). We have thus reached
the apparently paradoxical result that the mirrors do not move under the effect of a passing gravitational
wave (if they start at rest).
One must not forget, however, that coordinates (and therefore also the coordinate acceleration) have
no physical meaning in GR. What is physically relevant is not the coordinate distance between the free
falling mirrors, but their proper distance (which is independent of the coordinates). It is indeed the proper
distance that determines the light travel time between mirrors, and thus in turn the phase difference (the
laser “fringes”) at the photodetector.
Let us assume for simplicity that the mirrors are on the x axis. Their proper distance is then
Z L Z L q  
√ 1 TT
q
proper TT TT
L = dx gxx = dx 1 + hxx ' L 1 + hxx (t, x = 0) ' L 1 + hxx (t, x = 0) , (248)
0 0 2
where L is the coordinate distance, and we have assumed not only that the metric perturbation is small,
but also that it has wavelength much larger than L. (This assumption is needed in the third step.)
The change δL in proper distance is therefore
δL 1 1 TT i j
' hTT
xx (t, z = 0) = hij u u (249)
L 2 2
where in the last step we have introduced the unit-norm vector u to denote the detector’s arm direction.
In this way, this equation is valid in a general reference frame.
An alternative (and more instructive) derivation of the same result can be obtained by resorting to the
geodesic deviation equation,

D2 v µ µ
= Rαβγ uα uβ v γ , (250)
dτ 2
µ
with Rαβγ the Riemann tensor, D/dτ the covariant derivative along the four velocity, and v µ = δxµ the
separation vector between two neighboring geodesics (the two free falling mirrors). Using FNCs attached
to one mirror, one has

gµν = ηµν + O x2 /λ2 ,

(251)
µ µ 2 2


u = δt + O x /λ (252)

where we have used the fact that the mirrors are in free fall (aµ = 0) and the curvature radius of the
spacetime (Minkowski plus a gravitational wave signal) is given by the signal’s wavelength λ.
Since the mirrors are at rest in these coordinates, v µ = δxµ = (0, Li ). Because the metric is locally
flat, δxi = i

L describes the proper distance between the mirrors, up to corrections of fractional order
O L2 /λ2 . From Eq. 250, one then obtains

d2 Li
= −Ritjt Lj . (253)
dt2
One can now note that the Riemann tensor of a perturbed Minkowski spacetime is a gauge invariant
quantity. This follows simply from the fact that the gauge transformation of a tensor is proportional
to the Lie derivative of the background tensor along the gauge generator. Since the Riemann tensor
vanishes on the Minkowski background, it follows that it is a gauge invariant quantity at linear order.
 40

More explicitly, one can evaluate the Riemann tensor for a perturbed flat spacetime using the results of
section II C, obtaining
1 1
Ritjt = − ḧTT + ∂i ∂j ψ + Σ̇(i,j) − θ̈δij . (254)
2 ij 2
in terms of the gauge invariant variables introduced in that section.
Using now Eqs. 58 and 68, it is clear that all the terms in Eq. 254 decay as 1/r3 (i.e. they correspond
to Netwonian and PN tidal forces), and the only term that survives at large distances from a source is
the first one. Therefore, one has
d2 Li
 
j 1 TT j 1
= −Ritjt L = ḧij L + O , (255)
dt2 2 r3
and integrating this equation one obtains Eq. 249. Note that again, we have implicitly made the assump-
tion that the signal’s wavelength

is much larger than the distance between the two mirrors, because we
have neglected terms O x2 /λ2 in Eqs. 251–252.
The phase difference at the photodetector depends of course on the changes in proper length (i.e. light
travel time) induced by the gravitational signal on the first (δL1 ) and second (δL2 ) arm, i.e. in the low
frequency (large wavelength) limit
δL1 − δL2
= hij
TT Dij = F+ h+ + F× h× , (256)
L
where
1 i j
Dij = (u u − v i v j ) (257)
2
(with ui and v i unit vectors in the arm directions) is the detector tensor and F+ , F× are the pattern
functions (which encode the detector’s response in various directions). The phase difference at the
photodetector is then explicitly
4πν 4πνL
∆φ = (δL1 − δL2 ) = (F+ h+ + F× h× ) , (258)
c c
where ν is the laser frequency. (Notice the presence of a factor 2 due to the round trip of the laser.)

Exercise 5: Consider an interferometer on the (x, y) plane and a gravitational wave coming
from the sky position φ, θ (spherical coordinates). Define the wave plus and cross polarizations with
respect to two unit vectors ex and ey orthogonal to the propagation direction n, and such that the triad
(ex , ey , n) is right-handed. Show that the detector’s pattern functions are
1
cos 2ζ cos2 θ + 1 cos 2φ − sin 2ζ cos θ sin 2φ


F+ =
2
and
1
sin 2ζ cos2 θ + 1 cos 2φ + cos 2ζ cos θ sin 2φ


F× =
2
where ζ < π/2 is the angle between ey and the intersection between the plane of the detector and the
plane orthogonal to the wave propagation direction.

B. A geometric interpretation of the polarizations

Thanks to Eq. 256, we can now provide a geometric interpretation of the two polarizations h+ and h× .
If one considers a ring of particles on the (x, y) plane and a signal traveling along the z-direction, the h+
and h× polarizations deform the ring in a distinct and characteristic way (Fig. 10). Also shown, in Fig. 11,
is the response of a ring to a right-handed or left-handed circularly polarized wave, i.e. one for which
the h+ and h× polarizations have a phase difference of ±π/2 and the same amplitude. A generic linear
polarization corresponds instead to h+ and h× having the same phase (although potentially different
amplitudes). Note that from Eq. 184, it follows that for an observer that sees a circular binary face-
on (ι = 0 or ι = π) the signal is circularly polarized; if instead the observer sees the binary edge-on
(ι = π/2), the polarization is linear. We also notice that additional polarization patterns are present
beyond GR [36]. Those arise because the scalar and vector degrees of freedom of Eq. 68 become dynamical
(c.f. also section II E).
 41

t=0 t=T/4 t=T/2 t=3T/4 t=T

t=3T/4 t=T
t=0 t=T/4 t=T/2

FIG. 10: The response of an originally circular ring of particles to a linearly polarized gravitational wave (h+ or h× ) of
period T traveling in the orthogonal direction into the page.

t=0 t=T/4 t=T/2 t=3T/4 t=T

t=0 t=T/4 t=T/2 t=T

t=3T/4

FIG. 11: The response of an originally circular ring of particles to a circularly polarized gravitational wave (with right or
left helicity) of period T traveling in the orthogonal direction into the page.

C. The response of a gravitational wave detector: the transfer function

Let us now relax the short wavelength assumption (λ  L) that we made when deriving the detector’s
response in the previous section. Let us then consider two free falling mirrors in a spacetime with metric

gµν = ηµν + hTT

µν , (259)

and integrate null geodesics (the lasers) back and forth between them. As derived in section VIII A, the
mirrors do not move (with respect to the coordinates) if the perturbation is written in the transverse and
traceless gauge.
Let us also make the simplifying assumption that hTT µν is given by a plane wave traveling along the
 42

z-axis:
0 0 0 0
 
0 cos 2ψ sin 2ψ 0
hTT
µν = h(t − z) . (260)
0 sin 2ψ − cos 2ψ 0
0 0 0 0

(Recall that a generic linear perturbation can always be decomposed in a superposition of plane waves.)
It is then clear that the spacetime has three Killing vectors
∂
k1 = , (261)
∂x
∂
k2 = , (262)
∂y
∂ ∂
k3 = + . (263)
∂t ∂z
One can then write a tetrad (carried by the mirrors) for the perturbed spacetime as

∂
e(0) = , (264)
  ∂t 
h ∂ ∂
e(1) = 1 − cos ψ + sin ψ + O(h)2 , (265)
2 ∂x ∂y
  
h ∂ ∂
e(2) = 1 + − sin ψ + cos ψ + O(h)2 , (266)
2 ∂x ∂y
∂
e(3) = , (267)
∂z
and write the null wave-vector of the laser as
 
σ = ν e(0) + sin θ(e(1) cos φ + e(2) sin φ) + cos θe(3) . (268)

Here, ν is the laser frequency as measured at the mirror, and θ and φ describe the laser’s direction.
In particular, θ is the angle between the arm of the detector and the gravitational wave’s propagation
direction.
Now, let us recall that the projections of the laser wave-vector on the three Killing vectors are conserved
along the laser’s null geodesics:

−σ · k3 = ν(1 − cos θ) = const , (269)

   
h h
σ · k1 = 1 + (σ · e(1) ) cos ψ − 1 − (σ · e(2) ) sin ψ = const , (270)
2 2
   
h h
σ · k2 = 1 + (σ · e(1) ) sin ψ + 1 − (σ · e(2) ) cos ψ = const . (271)
2 2
p
This implies also conservation of the combination (σ · k1 )2 + (σ · k2 )2 , i.e.
 
p 1/2 h
(σ · k1 )2 + (σ · k2 )2 = (1 + h)(σ · e(1) )2 + (1 − h)(σ · e(2) )2

≈ ν sin θ 1 + cos 2φ = const .
2
(272)
Note that this conserved quantity can be rewritten in more compact form as
   
h Q
ν sin θ 1 + cos 2φ = ν sin θ 1 + = const , (273)
2 2
hTT i j
ij n n
Q= , (274)
1 − cos2 θ
where n = sin θ cos(φ + ψ)∂x + sin θ sin(φ + ψ)∂y + cos θ∂z is the detector’s arm direction (from the first
to the second mirror) expressed in Cartesian coordinates [this can be obtained by rewriting Eq. 268 using
Eqs. 264–267 and taking h → 0].
 43

FIG. 12: Sketch of the geometry of the system leading to Eq. 279, for y = 0.

Let us consider now a laser photon starting from first mirror (time “0”), bouncing on the second mirror
(time “1”) and returning finally to the first mirror (time “2”). From the conservation laws written above
one then obtains
ν0 (1 − cos θ0 ) = ν1 (1 − cos θ1 ) , (275)
   
Q0 Q1
ν0 sin θ0 1 + = ν1 sin θ1 1 + . (276)
2 2

Writing cos θ1 = cos θ0 + δ cos θ + O(h)2 , sin θ1 = sin θ0 + δ sin θ + O(h)2 and ν1 = ν0 + δν + O(h)2 , with
δ cos θ, δ sin θ and δν of order O(h), one can linearize these equations and obtain
ν1 − ν0 δν 1 + cos θ0
= = (Q0 − Q1 ) (277)
ν0 ν0 2
for the change in laser frequency between the first and second mirror. Performing the same calculation
to compute the change in laser frequency when the photon travels in the opposite direction (from the
second mirror to the first), one finds
ν2 − ν1 1 − cos θ0
= (Q1 − Q2 ) . (278)
ν1 2
Note that the different sign of cos θ0 in Eqs. 278 and 277 appears because during the return trip the
photon’s propagation is along −n.
Combining these results, one obtains that the change ∆ν in laser frequency at the first mirror after a
round trip is
 
∆ν ν2 − ν0 1 τ (1 − cos θ) 1
= = (1 + cos θ)Q(t) − cos θ Q t + − (1 − cos θ)Q(t + τ ) (279)
ν ν0 2 2 2
where we have dropped the index of the angle θ since the differences among cos θ0 , cos θ1 and cos θ2
appear at higher order, and we have used the fact that Q0 = Q(t), Q1 = Q(t + τ (1 − cos θ)/2) and
Q2 = Q(t + τ ) with τ the laser round trip time (c.f. Fig. 12). From this expression, one can then get the
change in the laser phase over a round trip,
Z
∆φ = −2π ∆νdt . (280)
 44

The quantity measured at the photodetector is then the difference between the phase changes ∆φ in the
first and second arm of the interferometer.
To compute this phase difference and get some physical insight, let us assume for simplicity θ = π/2,
which yields
ν
∆ν = [h(t) − h(t + τ )] , (281)
2
with h = hTT i j
ij n n . For a monochromatic wave h(t) = h0 sin (2πf t) one then has
Z
h0 ν h  τ i
∆φ = −2π ∆νdt = sin (f πτ ) sin 2πf t + , (282)
f 2
which can then be rewritten more simply as
 τ ν
∆φ = h t + sin (f πτ ). (283)
2 f
Let us interpret this result in terms of a proper distance change δL between the two mirrors, which
produces a change in the laser phase (in a round trip) given by

δL
∆φ = 4πν . (284)
c
Let us think about this distance change as due to an “effective” strain heff :
δL 1
= heff ,
L 2
where the length can be rewritten in terms of the round trip time τ , L = τ c/2. Replacing in Eq. 284,
one then obtains

∆φ = πντ heff . (285)

By comparing to Eq. 282, one gets

heff = h T (f ) , (286)
sin (f πτ )
T (f ) = , (287)
f πτ

where we have introduced the “transfer function” T (f ). This calculation can be generalized to more
generic propagation directions θ 6= π/2.
In summary, we can therefore write the full response of a detector to a Fourier component of frequency
f of a gravitational wave signal as
4πνL
∆φ = [F+ h+ (f ) + F× h× (f )] T (f ) . (288)
c
As can be seen from Eq. 287, at low frequencies T (f ) ≈ 1 and this expression reduces to the low
frequency response of Eq. 258. At high frequencies f  1/τ , the transfer function instead decays as
1/(f τ ), modulated by the oscillations of the numerator of Eq. 287. This explains why the frequency
window of gravitational interferometers scales with their armlength L = cτ .

IX. GRAVITATIONAL WAVE DATA ANALYSIS

As we have seen in the previous section, gravitational wave interferometers measure the interference
pattern (i.e. the phase difference) between laser beams traveling between mirrors. The effect of a gravita-
tional signal, however, is generally so small that it is also crucial to adequately understand the statistical
error (“noise”) affecting this measurement, and to devise statistical techniques to characterize the signal
and its astrophysical source.
 45

FIG. 13: Various contributions to the Advanced LIGO differential armlength (DARM) caused by the noise. This
quantity is proportional to the (square root of the) power spectral density. Taken from Ref. [37].

A. Gaussian noise and power spectral density

A detector’s output can be written as s(t) = h(t) + n(t), where h(t) is the signal and n(t) is the
instrumental noise. Defining the Fourier transform of a time series A(t) as
Z ∞
Ã(f ) = A(t)e2πif t dt , (289)
−∞

the signal is given by h̃(f ) = [F+ h+ (f ) + F× h× (f )] T (f ). As for the noise, a common assumption, which
turns out to be a good approximation for gravitational wave interferometers, is that of stationarity and
Gaussianity.
Stationarity implies that the statistical properties of the noise are time independent. Introducing the
ensemble average h. . .i, i.e. the average over all possible noise realizations 10 , stationarity implies in
particular that hn(t)n(t0 )i should only depend on the time difference τ = t0 − t, i.e. hn(t)n(t0 )i = W (τ )/2.
(Note also that hn(t)i = 0). Let us see what this implies for the quantity hñ∗ (f )ñ(f 0 )i. Using the
definition of Fourier transform, one has
Z Z
0 0
∗ 0
hñ (f )ñ(f )i = dtdt0 hn(t)n(t0 )ie2πi(f t −f t)
Z Z
1 0 0 1
= dtdτ W (τ )e2πi[(f −f )t+f τ ] = S(f )δ(f − f 0 ) , (290)
2 2
where S(f ) is the Fourier transform of W (τ ) and is known as the (single-sided) power spectral density
of the noise. One can also check easily that
Z ∞
2 W (0)
hn(t) i = = S(f )df . (291)
2 0

10 By the ergodic theorem, ensemble averages can be replaced with time averages.
 46

10-15 Acceleration
Laser
10-16
Total

10-17

10-18

10-19

10-20

10-4 0.001 0.010 0.100 1

FIG. 14: Various contributions to the LISA power spectral density. At high frequencies, one can see the degradation and
the oscillations due to the transfer function.

The assumption of Gaussian noise then amounts to saying that each Fourier component has a Gaussian
probability distribution function, with variance given by Eq. 290. In more detail, because detectors only
observe for a finite time T and thus Fourier transforms are implemented as discrete transformations, the
frequencies fi at which measurements are performed are spaced by ∆f = 1/T . In terms of these discrete
frequencies, Eq. 290 becomes
1
hñ∗ (fi )ñ(fj )i = S(fi )δij , (292)
2∆f

and the probability P (n) of having a realization ñ(f ) of the noise is

2|ñ(fi )|2 ∆f 2∆f [(Re ñ(fi ))2 + (Im ñ(fi ))2 ]

Y   Y  
P (n) ∝ exp − = exp −
i
S(fi ) i
S(fi )
 Z ∞
|ñ(f )|2

= exp −2 df . (293)
0 S(f )

Here, the index i spans all positive frequencies in the data [since n(t) is real, ñ∗ (f ) = ñ(−f ) and the
signal at negative frequencies follows from that at positive ones].
As can be seen by changing variables in this equation, while the real and imaginary parts of ñ(f ) are
Gaussian distributed, the phase of ñ(f ) is uniformly distributed, while the norm r(f ) = |ñ(f )| follows
the Rayleigh distribution p[r(fi )] ∝ exp{−2∆f [r(fi )]2 /S(fi )}r(fi ). Introducing the internal product
∞
df Ã∗ (f )B̃(f )
Z
(A | B) ≡ 4 Re , (294)
0 S(f )

between two real functions A(t) and B(t), one can rewrite Eq. 293 as
 
1
P (n) ∝ exp − (n|n) . (295)
2

As we will see, this internal product will also simplify many expressions below.
The behavior of the power spectral density of the noise depends critically on the detector. For ground
based interferometers such as LIGO and Virgo, at low frequencies the noise originates chiefly from the
laser’s radiation pressure, the seismic noise and the thermal noise of the suspensions; at mid-frequencies
 47

the main contributions are the thermal noise from the mirror coatings and the laser noise (radiation
pressure and shot noise); at high frequencies the shot noise dominates (see Fig. 13).
For space-borne detectors like LISA (Fig. 14), the limitations to the interferometer’s sensitivity come
from spurious accelerations of the test masses (due e.g. to cosmic rays, residual gas in the housing,
temperature fluctuations) at low frequencies; from the laser shot noise at mid-frequencies; and from
the antenna transfer function at high frequencies (c.f. the oscillations in Fig. 14, which are due to the
numerator of Eq. 287).

1. Detection in the presence of noise

In order to disentangle the gravitational wave signal from the noise, many techniques have been put
forward, among which one of the most popular is match filtering. The latter essentially amounts to cross
correlating the detector’s noisy output

s(t) = h(t) + n(t) (296)

with a bank of templates h(t, θ). Here, the vector θ denotes the parameters of the source (i.e. for quasi-
circular binaries, the masses, the spins, the distance, the initial phase, the merger time, the inclination,
the sky position and the polarization angle). It is quite intuitive that the cross correlation
Z
dt s(t)h(t, θ) (297)

will “on average” be maximized if the template matches the signal. Taking in fact the average of this
equation, one obtains
Z Z
dt hs(t)h(t, θ)i = dt h(t)h(t, θ) , (298)

and the second integral is small if the signal and template do not match (because in that case the integrand
is highly oscillatory).
Let us try to formalize this statement (see e.g. Ref. [38]) by defining a generic filter
Z
 = dt A(t)K(t) (299)

where A(t) is the time series being filters, and K(t) a filter function. Let us try to define a signal-to-noise
ratio S/N based on this filter. It is natural to take the signal S as the filter of the detector’s output s(t),
averaged over many realizations of the noise:
Z Z Z
S = hŝi = dt hs(t)iK(t) = dt h(t)K(t) = df h̃(f )K̃ ∗ (f ) . (300)

As for the denominator N , since hn̂i = 0, we can define instead

Z Z
1
N 2 = hn̂2 i = dtdt0 K(t)K(t0 )hn(t)n(t0 )i = dtdt0 K(t)K(t0 )W (t0 − t)
2
Z Z
1 0 1
= dtdt0 df K(t)K(t0 )S(f )e−2πif (t −t) = df |K̃(f )|2 S(f ) . (301)
2 2
The ratio S/N is then
R +∞
S −∞
df h̃(f )K̃ ∗ (f )
=h i1/2 . (302)
N R +∞
2
−∞
df (1/2)S(f )|K(f )|

By introducing the internal product of Eq. 294, and using the fact that for a real function A one has
Ã(−f ) = Ã∗ (f ), one can rewrite

S (h|u)
= , (303)
N (u|u)1/2
 48

where we have defined

1
ũ(f ) = S(f )K̃(f ) . (304)
2
Clearly, Eq. 303 is maximized if u and h are parallel, i.e. the filter yielding the optimal signal-to-noise
ratio is

h̃(f )
K̃(f ) ∝ . (305)
S(f )

Therefore, if one is searching for a signal, the optimal filter is the template perfectly matching the signal,
weighted by the noise power spectral density. This is known as Wiener’s optimal filter theorem. Replacing
back into Eq. 303, one obtains that the optimal signal-to-noise ratio is simply given by
2 ∞
|h̃(f )|2
 Z
S
= (h|h) = 4 df . (306)
N 0 S(f )

B. The signal-to-noise ratio for inspiraling binaries

Equation 306 allows one to compute the optimal signal-to-noise ratio, if the frequency domain signal
is known. In the special case of inspiraling binaries, the time domain signal is given, at the lowest PN
order, by Eq. 184. The evolution of the gravitational wave frequency fgw is instead described by Eq. 189,
which yields
 3/8
1 −5/8 5
fgw = M , (307)
π c 256(tc − t)

where we have identified the integration constant tc with the the coalescence time. This identification
is justified because it yields a frequency diverging at t = tc . This divergence is not physical, but simply
signals that the quadrupole formula, used to derive Eq. 189, breaks down. It is natural to assume that this
breakdown happens at the merger, since the assumption of a binary system, implicit in the quadrupole
formula, ceases to apply there.
As we have seen in section VIII, the detector is sensitive to a linear combination of h+ and h× ,
modulated by the transfer function. Neglecting for the moment the transfer function – because its effect
is small for LIGO/Virgo, while for LISA it is usually included in the power spectral density of the noise
(c.f. Fig. 14) – we can write the detector’s response, using Eq. 184, as
5/3
[πfgw (t)]2/3 1 + cos2 ι
 
4Mc
h(t) = F+ h+ + F× h× = F+ cos Φ(t) + cos ιF× sin Φ(t) (308)
D 2
Rt
where Φ(t) = t0 dt0 2πfgw (t0 ).
Note that although we derived Eq. 308 in flat space, it is possible to show (see e.g. Ref. [38]) that it is
valid also in a Robertson-Walker spacetime, provided that the frequency fgw is the one measured at the
detector [where it is redshifted by a factor 1/(1 + z) relative to the frequency at the source], the masses
are redshifted, i.e. Mc → Mc (1 + z), and D is interpreted as the luminosity distance DL .
To perform the Fourier transform analytically, one can use the stationary phase approximation. Con-
sidering first the h+ polarization, let us rewrite Eq. 308 as

A(t) h iΦ(t) i
h(t) = e + e−iΦ(t) (309)
2
by introducing an amplitude A(t), the Fourier transform is
Z
A(t) 2πif t h iΦ(t) i
h̃(f ) = dt e e + e−iΦ(t) . (310)
2

The main contribution to this integral comes from the regions the phase does not change rapidly (the
contribution of the other regions is negligible because the fast oscillations average to zero). Note that
we are only interested in computing h̃(f ) for f > 0 since h̃(−f ) = h̃∗ (f ) [because h(t) is real]. Then,
 49

because Φ̇ = fgw > 0, the first term in square brackets in Eq. 310 has no stationary points, and we can
thus neglect it. As for the second term, let us denote by t∗ the time at which the phase is stationary, i.e.
Φ̇(t∗ ) = 2πf . Taylor expanding near this stationary point, one obtains
Z s Z ∞
1 i[2πf t∗ −Φ(t∗ )] 2
−iΦ̈(t∗ )(t−t∗ ) /2 1 i[2πf t∗ −Φ(t∗ )] 2 2
h̃(f ) ≈ A(t∗ )e dt e = A(t∗ )e dx e−i x
2 2 Φ̈(t∗ ) −∞
s
1 2π
= A(t∗ )ei[2πf t∗ −Φ(t∗ )−π/4] . (311)
2 Φ̈(t∗ )

One can then express t∗ as a function of f by solving f = Φ̇(t∗ )/(2π) = fgw (t∗ ). Using Eq. 307 one then
has
5
tc − t∗ = M −5/3 (πf )−8/3 , (312)
256 c
which can be used into Eq. 311 to obtain the amplitude and phase as functions of frequency. Doing
the same calculation for h× yields the same result, but with an additional factor i (i.e. a π/2 phase
difference). The final result then reads
 1/2
1 5
h̃(f ) = π −2/3 Mc5/6 f −7/6 Qeiψ , (313)
DL 24
ψ = 2πf t∗ − Φ(t∗ ) − π/4 (314)

with t∗ given by Eq. 312 and

1 + cos2 ι
Q = F+ + i cos ιF× . (315)
2
Note that we have expressed this final result in terms of DL , and thus f and Mc must be interpreted as
the detector-frame frequency and the redshifted chirp mass, respectively.
One can then compute the signal-to-noise ratio using Eq. 306 as
2 5/3 fmax
f −7/3
  Z
S 5 1 1 GMc 2
= 2 |Q| df , (316)
N 6 π 4/3 DL c3 0 S(f )

where for clarity we have reinstated G and c. The factor Q depends on the sky position and orientation
of the source (c.f. Eq. 315). For sources with optimal sky position and inclination (ι = θ = 0) one
has Q = 1, while an average over sky position and inclination (assuming isotropic distributions) yields
hQ2 i = 4/25. Signal-to-noise ratios allowing detection are typically S/N & 8.

Exercise 6: Using the stationary phase approximation, show that the signal-to-noise ratio of a
quasi-monochromatic signal
" ! #
√ f˙0 t
h(t) = F+ h+ (t) + F× h× (t) ≡ 2A cos 2π f0 + t + φ0 (317)
2

is given by
2
2A2 T

S
= , (318)
N S(f0 )

where T is the observation time. [Hint: the maximum integration frequency depends on the observation
time.]

C. Parameter estimation

Let us now suppose that a gravitational wave signal has been detected in the data s, and let us examine
whether the parameters θ of the source can be inferred. In principle, the probability distribution function
 50

for the parameters conditional on the data, i.e. the “posterior distribution” P (θ|s), is given by Bayes
theorem in terms of the likelihood P (s|θ) and the “prior distribution” P (θ):

P (θ|s) ∝ P (s|θ)P (θ) . (319)

Here, the prior distribution encodes the knowledge on the parameters before data are taken (e.g. it can be
taken to be uniform if one wants to remain agnostic). As for the likelihood, i.e. the probability to obtain
the data s given the waveform parameters θ, it can be computed from Eq. 295 by replacing n = s − h as
per Eq. 296:
   
1 1
P (s|θ) ∝ exp − (s − hθ |s − hθ ) ∝ exp − (hθ |hθ ) + (hθ |s) , (320)
2 2
where hθ is the template.
Given the large dimensionality of the parameter space (up to 15 parameters for quasi-circular binaries),
the posteriors predicted by Eq. 319 cannot be plotted or sampled by brute force. A more sophisticated
approach based on Markov-Chain-Monte-Carlo or nested sampling techniques is needed, which is outside
the scope of these simple notes. It is possible, however, to expand the posteriors around their maximum
in a Taylor series, obtaining the Fisher matrix of the system. In more detail, let consider the case of
uninformative (uniform) priors, P (θ) = const, and assume that the posterior distribution peaks at θ = θ̄.
Imposing ∂P (s|θ)/∂θi = 0 at the peak, one then gets
! !
∂hθ ∂hθ
(θ̄) s − (θ̄) hθ (θ̄) = 0 . (321)
∂θi ∂θi

Expanding now the posterior distribution at quadratic order near the peak and assuming a high signal-
to-noise ratio [which allows for approximating s ≈ hθ (θ̄)], one obtains
 
1
P (θ|s) ∝ exp − Γij (θi − θ̄i )(θj − θ̄j ) , (322)
2
where we have introduced the Fisher matrix
!
∂hθ ∂hθ
Γij = (θ̄) (θ̄) . (323)
∂θi ∂θj

In the high signal-to-noise ratio limit, the posterior distribution near the peak is therefore described by
a multivariate Gaussian with covariance matrix

Cij = (Γ−1 )ij . (324)

While not useful to analyze real data, the Fisher matrix formalism can be used to obtain predictions
about how well the parameters are measurable by a given detector. This information is encoded in
the covariance matrix. Note that eliminating the row and column corresponding to a parameter in the
covariance matrix amounts to marginalizing over that parameter. Instead, eliminating the row and
column corresponding to a parameter in the Fisher matrix corresponds to assuming that parameter is
known exactly.

Exercise 7: Consider a quasicircular system of two black holes with masses m1,2 of 30 and
35M at 400 Mpc luminosity distance. Compute the signal-to-noise ratio for a detector with power
spectral density

S = 10−49 [x−4.14 − 5x−2 + 111(1 − x2 + x4 /2)/(1 + x2 /2)]Hz−1 ,

with x = f /(215 Hz), tapering the inspiral waveform of Eq. 313 at the ISCO of a particle with mass
µ = m1 m2 /(m1 + m2 ) around a Schwarzschild black hole of mass M = m1 + m2 . Assume optimal
inclination and sky position, and no spins.

Exercise 8: For the quasi-monochromatic signal of Exercise 6, assume f0 = 80 Hz, f˙0 from
Eq. 189 and compute the amplitude A that gives S/N = 20 for an observation time of one year and the
detector of Exercise 7. For such a source, compute the covariance matrix of the parameters (A, f0 and
f˙0 ).
 51

Acknowledgements

I would like to thank a number of friends and colleagues for insightful conversations over the years
on topics related to the physics of gravitational waves, and in particular Stas Babak, Emanuele Berti,
Luc Blanchet, Alessandra Buonanno, Vitor Cardoso, Guillaume Faye, Scott Hughes, Luis Lehner, Eric
Poisson and Luciano Rezzolla. Special thanks should also be extended to Flavio Riccardi and Marika
Giulietti, who transcribed my course at SISSA in May 2020 into a set of notes, the structure of some parts
of which resembles these. This work was supported by the European Union’s H2020 ERC Consolidator
Grant “GRavity from Astrophysical to Microscopic Scales” (Grant No. GRAMS-815673), and by the EU
Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie Grant Agreement
No. 101007855.

[1] B. P. Abbott et al. [LIGO Scientific and Virgo], “Observation of Gravitational Waves from a Binary Black Hole
Merger,” Phys. Rev. Lett. 116, no.6, 061102 (2016) doi:10.1103/PhysRevLett.116.061102 [arXiv:1602.03837
[gr-qc]].
[2] M. J. Pfenning and E. Poisson, Phys. Rev. D 65, 084001 (2002) doi:10.1103/PhysRevD.65.084001 [arXiv:gr-
qc/0012057 [gr-qc]].
[3] E. E. Flanagan and S. A. Hughes, New J. Phys. 7, 204 (2005) doi:10.1088/1367-2630/7/1/204 [arXiv:gr-
qc/0501041 [gr-qc]].
[4] V. F. Mukhanov, H. A. Feldman and R. H. Brandenberger, Phys. Rept. 215, 203-333 (1992) doi:10.1016/0370-
1573(92)90044-Z
[5] J. M. Bardeen, Phys. Rev. D 22, 1882-1905 (1980) doi:10.1103/PhysRevD.22.1882
[6] E. Barausse and L. Lehner, “Post-Newtonian approach to black hole-fluid systems,” Phys. Rev. D 88, no.2,
024029 (2013) doi:10.1103/PhysRevD.88.024029
[7] C. M. Will, “Theory and Experiment in Gravitational Physics,” Cambridge University Press, 2018
[8] E. Poisson and C. M. Will, “Gravity: Newtonian, Post-Newtonian, Relativistic”, Cambridge U. Press, 2014
[9] M. Bonetti and E. Barausse, “Post-Newtonian constraints on Lorentz-violating gravity theories with a
MOND phenomenology,” Phys. Rev. D 91, 084053 (2015) [erratum: Phys. Rev. D 93, 029901 (2016)]
doi:10.1103/PhysRevD.91.084053 [arXiv:1502.05554 [gr-qc]].
[10] M. Bonetti, E. Barausse, G. Faye, F. Haardt and A. Sesana, “About gravitational-wave generation by a three-
body system,” Class. Quant. Grav. 34, no.21, 215004 (2017) doi:10.1088/1361-6382/aa8da5 [arXiv:1707.04902
[gr-qc]].
[11] E. Poisson, “A Relativist’s Toolkit: The Mathematics of Black-Hole Mechanics,” Cambridge University Press,
2009, doi:10.1017/CBO9780511606601
[12] F. de Felice and C. J. S. Clarke, Cambridge University Press, 1992, ISBN 978-0-521-42908-5, 978-0-521-
26639-0
[13] E. Bertschinger, “Stress-Energy Pseudotensors and Gravitational Radiation Power”, http://web.mit.edu/
edbert/GR/gr7.pdf
[14] C. W. Misner, K. S. Thorne and J. A. Wheeler, “Gravitation,” W. H. Freeman, 1973, ISBN 978-0-7167-0344-0,
978-0-691-17779-3
[15] S. L. Shapiro and S. A. Teukolsky, “Black holes, white dwarfs, and neutron stars: The physics of compact
objects,”, 1983, doi:10.1002/9783527617661
[16] R. M. Wald, “General Relativity,” Chicago Univ. Pr., 1984, doi:10.7208/chicago/9780226870373.001.0001
[17] P. K. Townsend, “Black holes: Lecture notes,” [arXiv:gr-qc/9707012 [gr-qc]].
[18] J. H. Taylor, L. A. Fowler and P. M. McCulloch, “Measurements of general relativistic effects in the binary
pulsar PSR 1913+16,” Nature 277, 437-440 (1979) doi:10.1038/277437a0
[19] A. Buonanno and T. Damour, Phys. Rev. D 59, 084006 (1999) doi:10.1103/PhysRevD.59.084006 [arXiv:gr-
qc/9811091 [gr-qc]].
[20] K. Akiyama et al. [Event Horizon Telescope], “First M87 Event Horizon Telescope Results. I. The
Shadow of the Supermassive Black Hole,” Astrophys. J. Lett. 875, L1 (2019) doi:10.3847/2041-8213/ab0ec7
[arXiv:1906.11238 [astro-ph.GA]].
[21] J. M. Bardeen, W. H. Press and S. A. Teukolsky, “Rotating black holes: Locally nonrotating frames, energy
extraction, and scalar synchrotron radiation,” Astrophys. J. 178, 347 (1972) doi:10.1086/151796
[22] S. Babak, A. Taracchini and A. Buonanno, Phys. Rev. D 95, no.2, 024010 (2017)
doi:10.1103/PhysRevD.95.024010 [arXiv:1607.05661 [gr-qc]].
[23] A. Le Tiec, E. Barausse and A. Buonanno, Phys. Rev. Lett. 108, 131103 (2012)
doi:10.1103/PhysRevLett.108.131103 [arXiv:1111.5609 [gr-qc]].
[24] K. S. Thorne, “Disk accretion onto a black hole. 2. Evolution of the hole.,” Astrophys. J. 191, 507-520 (1974)
doi:10.1086/152991
[25] E. Barausse, V. Morozova and L. Rezzolla, “On the mass radiated by coalescing black-hole binaries,” Astro-
phys. J. 758, 63 (2012) doi:10.1088/0004-637X/758/1/63 [arXiv:1206.3803 [gr-qc]].
 52

[26] M. Campanelli, C. O. Lousto and Y. Zlochower, “Spinning-black-hole binaries: The orbital hang up,” Phys.
Rev. D 74, 041501 (2006) doi:10.1103/PhysRevD.74.041501 [arXiv:gr-qc/0604012 [gr-qc]].
[27] E. Berti, V. Cardoso and A. O. Starinets, “Quasinormal modes of black holes and black branes,” Class.
Quant. Grav. 26, 163001 (2009) doi:10.1088/0264-9381/26/16/163001 [arXiv:0905.2975 [gr-qc]].
[28] T. Regge and J. A. Wheeler, “Stability of a Schwarzschild singularity,” Phys. Rev. 108, 1063-1069 (1957)
doi:10.1103/PhysRev.108.1063
[29] F. J. Zerilli, “Effective potential for even parity Regge-Wheeler gravitational perturbation equations,” Phys.
Rev. Lett. 24, 737-738 (1970) doi:10.1103/PhysRevLett.24.737
[30] V. Ferrari and B. Mashhoon, “New approach to the quasinormal modes of a black hole,” Phys. Rev. D 30,
295-304 (1984) doi:10.1103/PhysRevD.30.295
[31] B. F. Schutz and C. M. Will, “Black hole normal modes: a semianalytic approach,” Astrophys. J. Lett. 291,
L33-L36 (1985) doi:10.1086/184453
[32] D. R. Brill, P. L. Chrzanowski, C. Martin Pereira, E. D. Fackerell and J. R. Ipser, “Solution of the
scalar wave equation in a kerr background by separation of variables,” Phys. Rev. D 5, 1913-1915 (1972)
doi:10.1103/PhysRevD.5.1913
[33] S. A. Teukolsky, “Perturbations of a rotating black hole. 1. Fundamental equations for gravitational electro-
magnetic and neutrino field perturbations,” Astrophys. J. 185, 635-647 (1973) doi:10.1086/152444
[34] S. L. Detweiler, “Black holes and gravitational waves. III - The resonant frequencies of rotating holes,”
Astrophys. J. 239, 292-295 (1980) doi:10.1086/158109
[35] J. Weber, “Gravitational Radiation,” Phys. Rev. Lett. 18, no.13, 498-501 (1967)
doi:10.1103/PhysRevLett.18.498
[36] D. M. Eardley, D. L. Lee, A. P. Lightman, R. V. Wagoner and C. M. Will, “Gravitational-wave observations as
a tool for testing relativistic gravity,” Phys. Rev. Lett. 30, 884-886 (1973), doi:10.1103/PhysRevLett.30.884
[37] A. Buikema et al., “Sensitivity and performance of the Advanced LIGO detectors in the third observing run”,
Phys. Rev. D 102, 062003 (2020), doi:10.1103/PhysRevD.102.062003
[38] M. Maggiore, “Gravitational Waves. Vol. 1: Theory and Experiments,” Oxford University Press, 2007, ISBN
978-0-19-857074-5, 978-0-19-852074-0