University of Amsterdam
MSc Astronomy & Astrophysics
Gravitation AstroParticle Physics Amsterdam
Master Thesis
Reliability of the Parameterised Test of General Relativity
on GW150914 and GW151226
An assessment of model waveforms and robustness of the test’s results
by
Janna Goldstein
10186832
August 2016
60 ECTS
September 2015 - July 2016
Supervisor:
Dr. Chris van den Broeck
Examiners:
Dr. Chris van den Broeck &
Dr. Dorothea Samtleben
Summary
The parameterised test of general relativity (GR) is a generic test of GR that uses a Bayesian approach to constraining parameterised deviations. It uses gravitational waves from coalescing compact binaries such as the two signals observed by the LIGO detectors in the past year, which consist of an inspiral stage, merger and ringdown of the final black hole. The analysis depends on the accuracy of the phenomenological gravitational waves model IMRPhenomPv2 that produces a full inspiral-merger-ringdown waveform. The model can only be directly verified using the limited number of numerical relativity waveforms that have been produced. In this project, the accuracy of IMRPhenomPv2 is studied more extensively by comparing it with the different, quasi-analytical model SEOBNRv3, which can be done for any choice of binary parameters. This comparison is done by calculating the faithfulness F , a measure of the overlap of the two waveforms. Around the parameters of both events GW150914 and GW151226, it was found that the models are mostly in good agreement, with F ≥ 0.96, and that lower matches can often be explained by a lack of calibration numerical relativity waveforms used in constructing the models. It was also found that generally, the stronger the precession effects in the spins of the black holes, the lower the faithfulness of the two models.
Additionally, the results of the parameterised test of GR on GW150914 are investigated in two tests. Seemingly, deviations from GR of 2–2.5σ are found for the inspiral region, which is suspected to be due to fluctuations in the noise. This is investigated by running the parameterised test on fifteen stretches of detector noise, in which a numerical relativity waveform is injected. The results are deemed inconclusive as the recovered signal to noise ratio of most runs is much higher than that of GW150914, however offsets are found at 1–2σ. A second test uses injections of model waveforms (IMRPhenomPv2) with built-in deviations from GR, to probe the sensitivity of the parameterised test. It was found that deviations in the inspiral regime are recovered somewhat less significantly than the injected deviations, which might be due to the previously mentioned noise effect. Some deviations in the merger-ringdown regime are measured very strongly, while others are not, indicating a possibly varying sensitivity of the parameterised test.
Contents
Summary 1
Acknowledgements 3
Introduction 4
1 Gravitational waves 5
1.1 Black hole binaries . . . 6
1.2 Detections of gravitational waves . . . 7
2 Testing general relativity 9 2.1 The strong-field regime . . . 9
2.2 The parameterised test of GR . . . 9
2.2.1 The nested sampling algorithm . . . 11
2.2.2 Results for GW150914 and GW151226 . . . 12
3 Waveforms 15 3.1 Waveform approximations . . . 15
3.1.1 Effective One Body . . . 16
3.1.2 Phenomenological models . . . 18
3.2 Overlaps . . . 21
3.2.1 Faithfulness . . . 21
3.2.2 Probing parameter space . . . 23
3.2.3 Results . . . 25
3.2.4 Discussion . . . 25
4 Reliability and sensitivity of the parameterised test 39 4.1 The effect of noise . . . 39
4.1.1 Results . . . 39
4.1.2 Discussion . . . 40
4.2 Sensitivity of the parameterised test . . . 40
4.2.1 Results . . . 43
4.2.2 Discussion . . . 43
Conclusions 47
A Results for the noise investigation 50
Acknowledgements
I would like to sincerely thank Dr. Chris van den Broeck, my supervisor, for his guidance, help and comments on my work during a busy time for gravitational waves. My gratitude extends as well to Dr. Dorothea Samtleben, for filling the role of examiner. Special thanks to Jeroen Meidam for all his help and patience. I have greatly appreciated being part of the gravitational waves group of Nikhef; thank you all. Finally, I wish to express my gratitude towards my friends and family for their support.
Introduction
The detection of gravitational waves from the merger of a binary black hole made by the LIGO Virgo Collaboration on the 14th of September 2015, GW150914 (Abbott et al., 2016b), has heralded a new era of gravitational wave astronomy. Together with the second detection GW151226 (Abbott et al., 2016a), also of two merging black holes, this shows that it is now possible to observe events that are hardly, or perhaps not at all, visible with electromagnetic telescopes. Given the distance to which the most energetic events can be seen, this provides an exciting prospect for both astrophysics and cosmology. From a physics perspective, gravitational waves are interesting because they probe strong, dynamical gravity. Thus, they provide the first possibility to test general relativity (GR) beyond the weak field regime.
Such a test is implemented in the parameterised test of GR, discussed in Chapter 2, which has been de-veloped in the gravitational waves group of Nikhef (National Institute of Subatomic Physics). In this work, several aspects of the reliability of the test’s methods and results are studied, focusing on its application to GW150914 and GW151226. The parameterised test is a search for generic deviation from GR in signals from coalescing binaries, consisting of neutron stars and/or black holes. In order to do so, deviation parameters are introduced in analytical waveforms that are approximations of gravitational waves as predicted by GR. These approximants are needed as it takes too long to calculate a large number of waveforms precisely, since this can only be done numerically. Naturally, it is necessary to study the reliability of these approximants, which is addressed in Chapter 3 of this thesis. Chapter 4 is an assessment of the robustness and sensitivity of the parameterised test; how much do random noise fluctuations affect the result? And how strong must a deviation from GR be to be picked up? Leading up to all this, the first chapter will provide an introduction on gravitational waves, specifically those from binary black hole mergers, and the detections GW150914 and GW151226.
Chapter 1
Gravitational waves
Gravitational waves were first predicted by Einstein (1916, 1918) as a vacuum solution to the equations of general relativity (GR): Rµν− 1 2gµνR = 8πG c4 Tµν, (1.1)
where gµν is the metric, Rµν is the Ricci curvature tensor, R the Ricci scalar, G the Newtonian gravitational
constant, c the speed of light and Tµν the stress-energy tensor. These equations are highly non-linear and
can not be solved exactly, save for a few special cases (for example, a single black hole). It is possible to simulate GR spacetime, which is called numerical relativity. However, far away from any strong curvature, for example at the site of a gravitational waves detector, the equations can be linearised. This is done by introducing a small perturbation h on a background Minkowski metric:
gµν = ηµν+ hµν, |h| 1. (1.2)
Furthermore, the detector site is considered a vacuum, so Tµν = 0. With this the Einstein equations can be
rewritten, using some gauge freedom as well, to obtain:
hµν = 0, (1.3)
or with the d’Alambertian written out: ∂2 ∂x2 + ∂2 ∂y2 + ∂2 ∂z2 hµν = 1 c2 ∂2 ∂t2hµν, (1.4)
which is a wave equation. The speed of these gravitational waves (GW) is c, the speed of light. The solutions are plane waves, and can be written (in the transverse-traceless gauge), as:
hT T = h× hm+ 0 h+ −hm× 0 0 0 0 cos(ω(t − z/c)). (1.5) Here, the wave is traveling in the z-direction, and has an angular frequency ω. h+ and h× are the two
polarisations of the wave. The effect of a wave with only non-zero h+is to expand and compress space in the
horizontal and vertical directions of the plane, perpendicular to the propagation direction. For a non-zero h×, the same thing happens under a rotation of 45◦.
The shape of the polarisations depends on the source of the gravitational waves, which can be anything with a mass quadrupole moment, or higher multipole moment (there is no gravitational mono- or dipole radiation). Therefore, even though the waves near a detector are described in linearised gravity, they do contain full non-linear GR from the source. This is what can be measured with GW detectors. Since gravity is such a weak force, these waves are generally incredibly tiny. Any detectable waves must come from a source with very strong spacetime curvature, and thus from very massive and compact objects with a rapidly varying quadrupole (or higher multipole) moment.
Figure 1.1: Model (template) waveform and numerical relativity waveform constructed using the estimated parameters of GW150914. On top, a schematic picture of the binary black hole coalescence is plotted, each stage corresponding to the region of the waveform below. Picture taken from Abbott et al. (2016b).
1.1
Black hole binaries
Compact binaries, consisting of two neutron stars and/or black holes, are among the most massive and compact sources of gravitational waves (GW) possible, and GWs from merging binary black holes are the only ones detected so far (Abbott et al., 2016b). In a compact binary, the orbit of the stars creates a mass quadrupole moment, such that GWs carry away energy. Thus, the orbit decays into a closer orbit with a more rapidly changing quadrupole, such that stronger GWs are emitted and so on, until eventually the objects merge. The decrease in the orbital period was famously observed in the Hulse-Taylor binary (Hulse and Taylor, 1975; Taylor and Weisberg, 1982), consisting of a neutron star and a pulsar. It is, however, far from the point of merger. On the other hand, binaries that are observed with gravitational waves are in the last stages of the inspiral, where the waves’ amplitude and frequency are high enough to be detectable by laser interferometers such as LIGO. The waves can be observed through the merger of the binary, after which there is an additional stage called the ringdown. This ringdown is caused by the damped oscillations of the final black hole that was created during merger. A schematic of the inspiral-merger-ringdown and the corresponding waveform is shown in Fig. 1.1, as given by Abbott et al. (2016b).
There is no exact solution for the waves from compact binaries, but an excellent approximation for the inspiral stage is made by expanding the equations in powers of v/c, where v is the characteristic velocity of the binary and c the speed of light. With Keplers third law, this velocity is given by:
v = (πGM fGW)1/3, (1.6)
where fGW is the gravitational wave frequency, which is twice the orbital frequency. Such an expansion is
called a post-Newtonian approximation, or PN for short. Due to a historical mishap, the PN orders are indicated with half-integers, referring to the power in v/c by twice as much; for example, 1.5PN would refer to the coefficients multiplying (v/c)3beyond leading order.
is false only for the last stages of the inspiral. In this adiabatic approximation, the system can be described with two differential equations:
˙ φ(t) = v 3 GM, (1.7a) ˙v = −F (v) E0(v) . (1.7b)
The first equation gives the time derivative of the orbital phase ˙φ, which is obtained by setting it equal to the orbital angular frequency ωorbit= πfGW, and rewriting fGW using the expression for v in Eq. (1.6). The
second equation simply equates the rate of binding energy E lost with the energy flux of the gravitational waves F , using that dE
dt = ∂E
∂v ∂v ∂t = E
0˙v. Both the flux F and the energy E can be written as a
post-Newtonian expansion in v/c, up to 3PN order for the energy (which only has integer orders) and 3.5PN order for the flux (which has both integer and half-integer orders). The fraction F (v)/E0(v) can be re-expanded in a number of different ways, which leads to different solutions due to the finite order of the expansions. As such, different PN approximations to the gravitational wave phase are obtained, which differ slightly (see for example Buonanno et al. (2009)). The general form is:
Φ(v) = n=7 X n=0 h φn+ φ(l)n ln v c iv c n−5 . (1.8)
Since the energy and flux are known up to 3.5PN order, the expansion goes up to n = 7, although in principle there are an infinite number of coefficients. These coefficients φnare functions of the binary parameters: two
masses and two spin vectors. The φ(l)n multiply an additional log term at some orders, for which the first
non-zero one is φ(l)5 at 2.5PN. The amplitude of the wave can also be calculated, but higher order corrections to the amplitude are more difficult to measure than for the phase, which makes it of lesser interest for data analysis purposes.
The approximation works well as long as v c and the adiabatic approximations holds true, which fails for the late stages of the inspiral of a coalescing binary. For example, the peak velocity in the binary black hole merger GW150914 was v ∼ 0.6c. Furthermore, the PN expansion does not describe the merger and ringdown stages of the wave. However, models have been developed with input from numerical simulations, that are able to describe the full inspiral-merger-ringdown (IMR) stages of the gravitational waves. One of these is the IMRPhenomPv2 model, for which the phase is consists of three parts for the inspiral, intermediate (late inspiral) and merger-ringdown regions of the model (Khan et al., 2016):
φ(f ) = φP N(M f ) +1η σ0+ σ1f +34σ2f4/3+35σ3f5/3+12σ4f2 for M f ≤ 0.018 1 η β0+ β1f + β2log(f ) − 1 3β3f −3 for 0.018 < M f ≤ 0.5 fRD 1 η α0+ α1f − α2f −1+4 3α3f 3/4+ α 4tan−1 f − αf 5fRD damp for M f > 0.5 fRD , (1.9) where G = c = 1 for brevity, and fRD is the ringdown frequency of the final balck hole. φP N is a
post-Newtonian approximant as a function of GW frequency instead of velocity using Eq. (1.6). The other coefficients are determined with fits to numerical relativity waveforms, {σi} for the inspiral, {βi} for the
intermediate and {αi} for the merger-ringdown regime. η = m1m2/M2 is the symmetric mass ratio. This
model, along with a different IMR model based on the Effective One Body formalism, will be further discussed in Chapter 3.
1.2
Detections of gravitational waves
In the first run (O1, from September 2015 till January 2016) of the Advanced LIGO (aLIGO) detectors, two detections of gravitational waves from merging binary black holes have been made with significance > 5.3σ, GW150914 and GW152126 (Abbott et al., 2016b,a). aLIGO consists of two laser interferometers, one in
Figure 1.2: Signals of binary black hole mergers in aLIGO’s O1 run. The plots show model reconstructions of the detected waveforms. The two detections are GW150914 with a (matched filter) SNR of 23.7 and GW151226 with an SNR of 13.0. LVT151012 is below the detection threshold with an SNR of 9.7. Figure taken from Abbott et al. (2016a).
Hanford, Washington and one in Livingston, Louisiana. Having two detectors in different locations allows for a better distinction between an incoherent glitch, such as a seismic disturbance at one of the sites, and a coherent signal of a gravitational wave passing through both detectors within a light travel time. The interferometers are built of two perpendicular 4 km arms with mirrors at the end. A laser beam is split in two and sent up and down both arms, tuned such that the reflected beams are exactly out of phase. A passing gravitational wave will lengthen one arm while shortening the other (the exact effect depending on the strength of the polarisations of the wave and the angle of incidence on the detector), thus changing the length of the path traveled by the laser. This changes the phase of the laser beams such that they interfere differently, which is measured with a photo detector.
Besides the two detections in O1, candidate signal LVT151012 was found to fall below the detection threshold with a significance of 1.7σ. (It was nevertheless investigated since it was observed coherently with both LIGO detectors and no specific disturbances in the data during the event were found.) All three signals match a model waveform of a merging binary black hole, for which parameters have been estimated (Abbott et al., 2016a). GW150914 was measured with a signal to noise ratio (SNR) of 23.7 and component masses m1 = 36.2+5.2−3.8 and m2 = 29.1+3.7−4.4, given with 90% confidence intervals. GW151226 is a quieter event with
an SNR of 13.0, and component masses m1= 14.2+8.3−3.7and m2= 7.5+2.3−2.3. The length of the measured signal
depends on its frequency range compared to the sensitivity band of the detector, which is roughly from 30 Hz to ∼ 104 Hz. This is determined by the total mass of the system, which is inversely proportional to the
frequency at merger. Hence, a much longer part of the inspiral is measured for GW151226, which has a lower mass and higher merger frequency, than GW150914. This difference can be seen in Fig. 1.2, which shows the model reconstructions of the three signals as given by Abbott et al. (2016a).
Chapter 2
Testing general relativity
2.1
The strong-field regime
In the century since Einstein’s publication of the theory of General Relativity (GR), GR has been successfully tested with many different observations (for a review, see Will (2014)). Since gravity is such a weak force, in almost any circumstance, only the linearised version of GR is necessary to describe physics. Therefore, the dynamical, strong-field regime of GR can not be tested with these observations. This holds even for the most relativistic binary neutron star that has been observed electromagnetically, the binary pulsar PSR J0737-3039. Gravitational wave observations of merging compact binaries, such as GW150914, do allow the study of the dynamical, strong-field regime. Compare the orbital velocity of PSR J0737-3039 of v ∼ 2 × 10−3c to that of GW150916, which is v ∼ 0.6c at its peak. The gravitational field strengths can be expressed with the dimensionless quantity GM/(c2R), where M is the binary’s total mass and R its separation, being 4.4 × 10−6
for PSR J0737-3039 and ∼ 0.1 for GW150914.
If reality is governed by some alternative theory to GR that manifests itself in the strong-field regime, this should leave a specific mark on gravitational waves emitted by binary coalescences. It would be relatively straightforward to devise a test that discriminates between said alternative theory and GR. However, a new test would be needed for each alternative theory, and even then possible unknown alternatives have not been accounted for. It is more desirable to devise a generic test, that can detect any possible deviation from GR. The TIGER method (Li et al., 2012) is an implementation of such a generic test. It uses Bayesian inference to calculate the ratio between the evidence for GR and the evidence for “not GR”, called the Bayes factor (the probability ratio of the two hypotheses with the ratio of prior probabilities taken out). It is a computationally expensive calculation, but part of its infrastructure can be used on its own to constrain deviations from GR. This part, which will be referred to as the parameterised test of GR, has already been performed using the data from GW150914 and GW151226 (Abbott et al., 2016d,a). Results of these tests can be found in section 2.2.2.
2.2
The parameterised test of GR
In the TIGER method, the hypotheses “gravity follows GR” and “gravity deviates from GR” are formalised using the coefficients of a waveform. For an inspiral gravitational wave signal, a post-Newtonian expansion could be used, which has convenient coefficients for each order in v/c (see Eq. (1.8)). For the inspiral-merger-ringdown (IMR) waveforms from coalescing binaries, the IMRPhenomPv2 model is used (see Eq. (1.9) and section 3.1.2), which has besides the PN coefficients {φi}, additional coefficients {σi}, {βi} and {αi}. Say
all the coefficients of the chosen model are the set {φn}. Then, the GR hypothesis in the TIGER method
is “All of the coefficients φn are as predicted by GR”. Every other possibility is captured by the non-GR
The phrase “one or more of the coefficients” can be described equivalently as “any of the non-empty subsets of {φn}” (the empty subset corresponds to the GR hypothesis). Each subset S can be associated
to a sub-hypothesis HS, “the coefficients in the set S deviate from the GR value, and other coefficients
do not”. For example, in a waveform model with two coefficients {φ1, φ2}, the non-emtpy subsets are
{{φ1}, {φ2}, {φ1, φ2}}, corresponding to the sub-hypotheses H1, H2, and H1,2. In order to parameterise the
deviations from GR, each coefficients φn is multiplied by a factor (1 + δ ˆφn) in the waveform model, such
that all δ ˆφn = 0 represents the GR waveform. The set {δ ˆφn} is called the deviation parameters.
In the full TIGER calculation, the evidence for each of the sub-hypotheses is calculated separately, after which they are combined to form the evidence for the non-GR hypothesis. (The non-GR hypothesis is the logical OR of the sub-hypotheses. Since they are logically disjoint, the evidences can be averaged over, see also Van Den Broeck (2014).) However, before calculating the combined Bayes factor, valuable results can be obtained from the evidence calculations on the sub-hypotheses, especially those corresponding to a single deviation parameter. The evidence for the sub-hypothesis Hk, corresponding to δ ˆφk 6= 0, is formally given
by the probability of measuring the data d, given that Hk and any background information I are true:
P (d|Hk, I) =
Z
d~θ p(d|~θ, Hk, I) p(~θ|Hk, I). (2.1)
The integral over ~θ represents a marginalisation over the parameter space of the model, which includes all the parameters from GR, such as the masses and spins of the binary, as well as the deviation parameter δ ˆφk.
The probability density p(~θ|Hk, I) is the prior probability distribution of the parameters, and p(d|~θ, Hk, I)
is the likelihood function. The likelihood is calculated from the data, by matching it with a modelled signal. This integral is evaluated using the nested sampling algorithm (Skilling, 2004), as explained in the following section 2.2.1.
Using the calculation of the evidence integral with nested sampling, it is trivial to extract the posterior probability distribution on any parameter as well. For the posterior of the deviation parameter δ ˆφk, with
Bayes’ theorem, we have:
p(δ ˆφk|d, Hk, I) =
p(d|δ ˆφk, Hk, I)p(δ ˆφk|Hk, I)
P (d|Hk, I)
. (2.2)
The denominator is given by Eq (2.1), and the numerator can be found by integrating the same equation over only the parameters found in GR, leaving out the integration over δ ˆφk:
p(d|δ ˆφk, Hk, I)p(δ ˆφk|Hk, I) =
Z
d~θGR p(d|~θ, Hk, I) p(~θ|Hk, I). (2.3)
If the data d is determined by GR, the distribution (2.2) should be peaked around δ ˆφk = 0, or possibly
with a small offset from zero due to noise in the data. Given the width of the distribution, one can constrain the value of δ ˆφk to a certain interval (around zero in the case of GR), and thus exclude any values outside
the support of the distribution. This then, is the parameterised test of GR; constraining each deviation coefficient by its posterior probability distribution.
The data d represents a signal from a single binary coalescence event, but data from multiple events can easily be combined. Given data from two distinct events, dA and dB, the combined posterior probability
distrubtion of the deviation parameter δ ˆφi can found using Bayes’ theorem. It is also assumed that the data
are independent, such that p(dA|dB) = p(dA) and vice versa, which yields:
p(δ ˆφk|dA, dB, I) =
p(δ ˆφk|dA, I)p(δ ˆφk|dB, I)
p(δ ˆφk, I)
. (2.4)
Equivalently, any number of posterior distributions can be combined by multiplying them, and dividing by the prior of δ ˆφk for any event after the first.
2.2.1
The nested sampling algorithm
The nested sampling algorithm can be used to calculate the integrated evidence from Eq. (2.1) by producing samples of the liklihood p(d|~θ, Hk, I) and the prior p(~θ|Hk, I) over parameter space. Conveniently, these
samples can also be used to calculate the posterior probability distribution of any of the parameters in ~θ, for example those for the parameterised test of GR as in Eq. (2.2), for which the total evidence is a normalisation factor. In the following, the algorithm will be explained by evaluating the evidence first and from there the posterior distribution, although the latter is of greater interest for the parameterised test; it is the same algorithm either way.
In evaluating the evidence integral, first, the likelihood and prior need to be quantified. The total prior is simply a multiplication of a chosen prior for each of the model parameters. Since there is no a priori reason to prefer any parameter values, the priors are mostly uniform distributions on an astrophysically or computationally sensible interval (Abbott et al., 2016c). The two angles determining the sky position (right ascension and declination), and the inclination of the source are naturally limited on different intervals ([0, π], [−π/2, π/2] and [0, π/2] respectively). The distance is the only non-flat prior, since it is chosen such that the distribution of sources is homogeneous over a volume. A maximum distance can be chosen such that even the brightest sources would not be detectable from further away, and so an interval of [1, 2000] Mpc is used. The spin magnitudes are naturally limited on [0, 1], and the spin orientations on the sphere; flat priors are chosen for both. For the initial searches, the component masses are limited on [1, 100] M , but
for later analyses this interval can be tightened to decrease computational cost, depending on the estimated parameters. For GW150914, a mass prior of [10, 80] M is used, and for GW151226 this is [3.022, 54.398]
M . Lastly, the prior for each deviation parameter is chosen wide enough so that it has no effect on the
measurement.
The likelihood function expresses how well the data fit the model waveform hk(~θ) predicted by Hk.
Assuming the data d consist of a signal and stationary, gaussian noise, the likelihood is defined as follows (Veitch and Vecchio, 2010; Van Den Broeck, 2014):
p(d|~θ, Hk, I) = N e−(d−hk(~θ)|d−hk(~θ))/2. (2.5)
Here, N is some normalisation factor. This holds for data from a single detector, but data from multiple detectors can be combined by multiplying the likelihood functions. The inner product that is used is given by: (h|g) = 4Re Z ∞ 0 ˜ h∗(f ) ˜g(f ) Sn(f ) df, (2.6)
where ˜h and ˜g are the Fourier transformed functions, and ˜h∗ represents the complex conjugate. Sn(f ) is
the noise power spectral density, i.e. the power of the noise per frequency bin df , as a function of frequency. Its form depends on the characteristics of the detector, and can be measured in the absence of a signal. In practice, either an analytical estimate or a fit based on noise measurements is used. A plot of a typical noise curve for the Advanced LIGO detector is given in figure 3.2.1. Since the noise gets loud for either low or high frequencies, the integral in Eq. (2.6) can be cut off on both sides. Furthermore, since the data are not continuous but a discrete set of points, the integral is replaced by a summation, and the continuous Fourier transform is replaced by the discrete version.
Having quantified the evidence integral (2.1), it is not trivial to calculate over a highly dimensional parameter space as the one for gravitational wave models. Veitch and Vecchio (2010) have described a method to make this computationally feasible using the nested sampling algorithm (Skilling, 2004). As the integral can not be analytically computed, it needs to be approximated by a sum over N points in parameter space: Z d~θ p(d|~θ, Hk, I) p(~θ|Hk, I) ≈ N X i=1 Liwi, (2.7)
where Liis the likelihood at the point ~θi(as in Eq. (2.5)), and wi is the fraction of the prior associated to it.
A simple way of distributing N values is to lay a regular grid over parameter space, since then each weight is given by wi= 1/N (for a flat prior). However, only a small fraction of the points contribute the most to
the total sum, since the likelihood is only high in a small section of parameter space, namely that where the model waveform is close to the signal. Thus, it would require many points to get an accurate approximation of the integral (see (Veitch and Vecchio, 2010) for a quantification). In the nested sampling algorithm, instead, only relatively few points are used at any time to sample the likelihood, and the results are used to evolve the collection of points to regions of higher likelihood iteratively. At each step, the contribution to the integral of the lowest likelihood point is calculated, after which it is replaced by randomly drawing a new point, on the condition that is has a higher likelihood than the original point.
A remaining difficulty is assigning the weights wi to the sample points. Imagine arranging the points by
increasing likelihood, and each point laying on an equal likelihood surface enclosing a volume in parameter space. The total volume of parameter space is assigned a prior probability of one, and so each smaller volume represents a fraction of the prior. Moving along the collection of points towards higher likelihood values, the enclosed fraction of the prior decreases with each step. The weight wi assigned to a point ~θi, then, is the
difference between its enclosed fraction and that of the previous point. These weights can not be determined exactly, but a probability distribution is found for the ratio between consecutive prior fractions, such that the average ratio can be calculated. This number is used to approximate the consecutive fractions, and thereby the weights wi, also using that the outermost possible sample point would enclose the total prior of
one. (For details, see Veitch and Vecchio (2010).)
The algorithm is terminated on the basis of the highest likelihood point that has been found so far. An estimate of the remaining evidence is obtained by assuming that all points are at this highest likelihood value. If this is less than a small fraction (e−5) of the total evidence accumulated so far, the iterations are stopped and the estimated remaining evidence is added to the total. This describes the complete outline of the algorithm, although it can be run multiple times (simultaniously) to obtain a more accurate result. In that case, each point ~θi is saved along with its likelihood value Li, after which the combined set of points is
rearranged and new weights wiare assigned.
As mentioned before, the calculation of the total evidence can also be used to find the posterior probability density of the deviation parameter(s). This is done by re-using the samples Liwi, and the total evidence as a
normalisation factor. Given equations (2.2) and (2.3), the postetior distribution of the deviation parameter δ ˆφk is given by:
p(δ ˆφk|d, Hk, I) =
R d~θGR p(d|~θ, Hk, I) p(~θ|Hk, I)
P (d|Hk, I)
(2.8) As in equation (2.7), the integral can by approximated by a sum over the posterior samples Liwi, only
here the samples are first divided over bins in δ ˆφk before they are summed. However, the weights need to
be corrected for the increased density of samples towards higher likelihood. For sample Liwi, where the
likelihood contour at Li encloses a fraction Xi of the prior, the sample density is boosted by a factor 1/Xi.
To compensate, corrected weights wiXi are used. Note that, outside of the parameterised test of GR, other
model parameters are estimated in the same way as a deviation parameter.
2.2.2
Results for GW150914 and GW151226
For the parameterised test of GR on data from GW150914 and GW151226, the waveform model IMRPhenomPv2 (see sections 3.1.2 and 3.1.2) was used. This phenomenological model is characterised by a set of coefficients for the phase and amplitude of the wave, both divided into three parts for the inspiral, intermediate and merger-ringdown regimes of the model. Only the phase coefficients are used in the parameterised test, since it is more accurately measured than the wave’s amplitude. Some coefficients are excluded because they have the same effect as a constant shift in the wave’s initial phase φ0 or starting time t0, and
there-fore can not be measured. Furthermore, the set of phenomenological coefficients that provide corrections to the post-Newtonian description of the inspiral phase are excluded as well. They are only determined with high uncertainties in the calibration of the model. As their GR predicted value is not very precise, it
would be difficult to find a significant deviation in any signal with these coefficients. For each remaining parameter in the testing set, a deviation coefficient is built in the model as explained in section 2.2. They are {δ ˆχ0, δ ˆχ1, δ ˆχ2, δ ˆχ3, δ ˆχ4, δ ˆχ5l, δ ˆχ6, δ ˆχ6l, δ ˆχ7} for the inspiral regime (see the inspiral parameters {φn} in
appendix A in Khan et al. (2016), only note that the logarithmic parts of φ5 and φ6 are treated as
sepa-rate coefficients indicated with a subscript l), {δ ˆβ2, δ ˆβ3} for the intermediate, and {δ ˆα2, δ ˆα3, δ ˆα4} for the
merger-ringdown regime (see Table V in Khan et al. (2016)).
The results of the parameterised test on GW150914 and GW151226 published in Abbott et al. (2016d,a), are given in Fig. 2.1. The combined posteriors, as in Eq. (2.4), are found in the bottom panel of the same figure. The posteriors of the deviation coefficients on the post-Newtonian parameters of GW150914 are not well centered around zero. This may be due to the short length of the inspiral part of the waveform that was measured, in combination with effects of the detector noise. For GW151226, for which a much longer part of the inspiral was in band for aLIGO, the posterior are much more centered on the GR value. An investigation on the effects of noise is discussed in section 4.1.
Figure 2.1: Results of the parameterised test of GR on GW150914 and GW151226, as published in Abbott et al. (2016a). The posteriors are given as violin plots, with bars indicating the 90% confidence interval. The top two rows are results from the individual events, while the bottom row shows the combined posteriors. Coefficients from the inspiral regime of the model IMRPhenomPv2 are indicated by their post-Newtonian order. The last panel contains the tested coefficients from the intermediate (βn) and merger-ringdown (αn)
Chapter 3
Waveforms
In the analysis of a gravitational wave (GW) signal from a binary coalescence, many different model wave-forms need to be generated. For example, to calculate posterior distributions for the binary’s parameters, template waveforms are needed across the possible parameter space to be compared to the measured wave-form. Exact solutions to the Einstein equations only exist for some very specific mass-energy configurations, such as a single black hole or a homogeneous and isotropic universe. For other cases they can be solved numerically, however, numerical relativity (NR) does not suffice for these GW analyses. It takes NR codes roughly a month to generate a single waveform. Moreover, in the parameterised test, a non-GR parameter is included in the search. Template non-GR waveforms can not be produced using NR, since it simulates GR spacetime.
Gravitational waves from coalescing binaries will have to be modelled with approximations. Two of these have already been discussed in Section 1.1, the first being the post-Newtonian (PN) expansion in powers of v/c. This is a very good approximation as long as the orbital velocity is not very high, which is true during the early inspiral of the binary, but not when the objects get closer to each other and merge. For signals from binary neutron stars, waveforms using the PN approximation are sufficient, because the merger happens at such a high frequency that it is not well measured by the detector anyway. On the other hand, for binaries with one or two black holes, the late inspiral and the merger can be in band of the detector. Additionally, the dynamics are complicated by the spins of the black hole(s) and the presence of a typical ringdown signal of the new black hole that was created during merger. As these more complicated parts of the signal contain a lot of interesting information - they are created in a stronger, more dynamical gravitational field - a waveform is required that can model this complexity.
Different approximation models exist, such as the IMRPhenomPv2 model that was briefly discussed in Section 1.1. Another important method used the Effective One Body (EOB) formalism, and both are explained in the first section 3.1 of this chapter. In both methods, the model is calibrated against a set of NR waveforms. It can be easily verified that the models are within reasonable error from both the NR waveforms used for calibration and others within the calibration parameter range. However, a different test is needed to assess the accuracy of the models at any point in parameter space, where no NR waveform exists. Section 3.2 contains the methods and results of a test that has been done based on comparing different approximant models to each other.
3.1
Waveform approximations
Both methods of approximating GR waveforms discussed in this section take the post-Newtonian expansion of GR as a starting point, and manage its accuracy by calibrating the waveform to NR results. As the PN model can not (accurately) describe the late inspiral, merger and ringdown of a GR waveform, it has to be extended in some way. An additional difficulty in modelling these waveforms are possible precession effects of the black hole spins. To leading order, these are due to spin-orbit coupling ~S · ~L, where ~S is one of the
black hole spins and ~L the orbital angular momentum of the binary. For spins aligned with ~L, there is no precession.
3.1.1
Effective One Body
The EOB formalism
The EOB formalism is a method that maps the dynamics of a two-body system (for example, a black hole binary), to the dynamics of one (extended) body and a test particle (for a more detailed explanation, on which the following is based, see Damour (2013), or Buonanno et al. (2009)). As a starting point, it uses the Hamiltonian of a two-body system in GR, that has been calculated using a post-Newtonian expansion up to 3PN order. The first, non-GR term (0PN) describes a ’test particle’ of mass µ that moves in the gravitational field of an ’external mass’ M . Here, M is the total mass of the two bodies and µ the reduced mass µ = m1m2/M . The idea behind EOB is to extend this description to GR by finding an effective
geometry for the test particle to move in, such that its dynamics are equivalent to the actual system. An ansatz is made for the effective geometry (Buonanno et al., 2009):
ds2ef f = −A(r)dt2+D(r) A(r)dr
2+ r2(dθ2+ sin(θ)2dφ2) (3.1)
Where A(r) and D(r) are written as post-Newtonian expansions: Ak(r) = i=k+1 X i=0 ai(η) ri (3.2a) Dk(r) = i=k X i=0 di(η) ri (3.2b)
ai(η) and di(η) are all functions of the symmetric mass ratio η = µ/M that need to be determined by the
mapping between the actual and effective dynamics. To what order k this can be done, depends on the PN order of the dynamics of the original problem; for 3PN, k = 4. The geometry should reduce to the Schwarzschil geometry for η → 0, as that describes the case of an actual test particle of zero mass around a single mass M . This requirement sets some of the coefficients to a fixed value.
The mapping between the two descriptions is inspired by quantum mechanics: the energy levels are com-puted as though they were quantised (only here the quantum numbers n and l are very large). The energy levels for the 3PN dynamics can be calculated from the Hamiltionian. Similarly, the energy levels of the effective description can be calculated from the geometry, with quantum numbers nef f and lef f. A unique
mapping between the energy levels of the actual and the effective dynamics can be found by identifying n = nef f, l = lef f. This solution then completes the description of the effective dynamics.
Having a description of the dynamics of the EOB system, the gravitational radiation needs to be intro-duced. This is done by adding a radiation reaction force in the equation of motion of the momentum in the azimuthal direction, assuming that the binary is circularised. After resumming a number of quantities, explicit expressions for the radiation reaction and the emitted waveform are obtained.
The waveform produced by the EOB method encompasses only the inspiral part of the binary coalescence. In order to complete it, a superposition of quasi normal modes (QNM) is pasted after it, which represents the ringdown of the Kerr black hole that was formed. The cut is made at the peak of the orbital frequency, which is identified as the time of effective merger tm. The coefficients of the QNM’s are adjusted so that the
two parts of the waveform (and its time derivative) match each other at a number of points chosen around tm. The complete waveform now consists of the EOB waveform up until tm, and the matched ringdown
EOBNR waveforms
The complete EOB waveform is constructed off of a post-Newtonian expansion of the dynamics of the binary. It is therefore not surprising that these waveforms are not quite accurate compared to numerical relativ-ity (NR) waveforms. In order to fix this, input from the NR waveforms is used. Higher order PN terms are added in the dynamics, although the coefficients are unknown. The complete EOB waveform is calcu-lated with the undetermined coefficients. Additionally, unknown non-quasi circular parameters are added to the waveform, to account for a departure from an adiabatic quasi-circular inspiral as described by the PN expansions. Then, the unknowns are determined by fitting the EOB waveform to NR data. The result is cal-ibrated EOBNR waveforms, that are much more accurate than the original EOB waveforms. (Damour, 2013) The EOBNR method has been extended to include spin, and the newest waveform model including spins aligned with the orbital angular momentum (i.e. non-precessing spins) is SEOBNRv2 (Taracchini et al., 2014). It has been calibrated using 27 of a total of 38 NR waveforms listed in Table 3.1.1. In the calibration set, mass ratios range from 1.0 to 8.0. Spins range from -0.95 to 0.98, but only for q = 1, whereas the largest spin value for unequal mass binaries is 0.5.
nr. q η χ1 χ2 1 1.0 0.25 0 0 2 1.5 0.24 0 0 3 2.0 0.22 0 0 4 3.0 0.19 0 0 5 4.0 0.16 0 0 6 5.0 0.14 0 0 7 6.0 0.12 0 0 8 8.0 0.099 0 0 9 1.0 0.25 0.98 0.98 10 1.0 0.25 0.97 0.97 11 1.0 0.25 0.95 0.95 12 1.0 0.25 -0.95 -0.95 13 1.0 0.25 0.9 0.9 14 1.0 0.25 -0.9 -0.9 15 1.0 0.25 0.85 0.85 16 1.0 0.25 0.8 0.8 17 1.0 0.25 -0.8 -0.8 18 1.0 0.25 0.6 0.6 19 1.0 0.25 -0.6 -0.6 nr. q η χ1 χ2 20 1.0 0.25 0.44 0.44 21 1.0 0.25 -0.44 -0.44 22 1.0 0.25 0.2 0.2 23 1.0 0.25 -0.2 -0.2 24 3.0 0.19 0.5 0.5 25 3.0 0.19 -0.5 -0.5 26 1.0 0.25 0.5 0.0 27 1.0 0.25 -0.5 0.0 28 1.5 0.24 0.5 0.0 29 1.5 0.24 -0.5 0.0 30 3.0 0.19 0.5 0.0 31 3.0 0.19 -0.5 0.0 32 5.0 0.14 0.5 0.0 33 5.0 0.14 -0.5 0.0 34 8.0 0.099 0.5 0.0 35 8.0 0.099 -0.5 0.0 36 1.5 0.24 0.5 -0.5 37 1.5 0.24 -0.5 0.5 38 2.0 0.22 0.6 0.0
Table 3.1: NR waveforms used in the calibration of the SEOBNRv2model (Taracchini et al., 2014). 1-8 are non-spinning waveforms, and 9-25 have equal spins (χ1= χ2).
Furthermore, a waveform model with precessing spins has been implemented under the name SEOBNRv3 (Pan et al., 2014). Its strategy is to compute the precession effects separately form the rest of the waveform, by constructing a frame that precesses along with the orbital angular momentum ~L. In this frame, the inspiral EOB waveform is produced using the spin components aligned with ~L only. Then, it is transformed to the frame aligned with the spin of the final black hole, in which the ringdown waveform (the superposition of QNMs) is added and matched around the merger time. Finally, another transformation yields the waveform in the binary’s original orientation. The non-precessing inspiral waveform is produced using the SEOBNRv2 model, and no further (precessing) calibration NR waveforms are used in SEOBNRv3.
Reduced order modelling
Although an EOBNR waveform is much faster generated than a full NR waveform, it is not a closed-form analytical formula as it involves numerically solving differential equations. For example, it takes SEOBNRv3 about 2.4 seconds to make a waveform for a binary black hole with total mass 60M (equal masses and zero
spins). This makes is still too slow to be of any use when calculating posteriors with the parameterised test. It is still useful for creating mock data by injecting a model waveform in detector noise, since then only one waveform needs to be calculated.
A method exists to turn a slower model, such as SEOBNRv2, into a faster one. Reduced Order Mod-elling (ROM) (P¨urrer, 2014) uses a basis of many different pre-calculated waveforms across parameter space and interpolates for any given point between them. The basis is chosen such that it captures all the most important features of the model. Different versions of a ROM model for SEOBNRv2 exist, of which SEOBNRv2ROMDoubleSpinis the most recent version that can take any two values for the aligned spins of the two compact objects. To compare with SEOBNRv3, this model produces a waveform with M = 60M
in about 0.036 seconds. It has been shown that the error of the ROM model relative to the original model is small compared to the error in the original model to NR waveforms (P¨urrer, 2014), hence it is just as accurate to use the ROM model.
Unfortunately, no ROM model exists for SEOBNRv3. Since this approximant uses precessing spins, there are four additional parameters compared to SEOBNRv2 (the x-components and y-components of both spins, although their in-plane angle is usually not very relevant). A ROM model would therefore require a basis of calculated waveforms over a parameter space with four additional dimensions, which is a non-trivial upgrade.
3.1.2
Phenomenological models
A different approach to making a fast waveform approximant is building a phenomenological model, for which the procedure (Husa et al., 2016; Khan et al., 2016) is outlined here. First, example waveforms are made by combining an approximate waveform for the inspiral - for example from a PN expansion - and numerical relativity (NR) data for the late inspiral, merger and ringdown. These hybrid waveforms are constructed in the frequency domain, as is the eventual phenomenological model, which is convenient for use in most analysis methods. For this purpose, the NR data is converted with a discrete Fourier transform. The parameters of the inspiral waveform are adjusted so that it matches the NR waveform well, over a chosen frequency window in which both are reliable approximations.
Then, these hybrid waveforms are used to calibrate a phenomenological ansatz. The ansatz is a closed-form analytical expression with, before the calibration, unknown phenomenological parameters. These are expressed as functions of the masses (m1, m2) and spins (χ1, χ2) of the binary. With these functions, any
point in parameter space can be used to obtain a waveform. Generally, the approximation is only valid in the region of parameter space that was covered by the hybrid waveforms used for calibration.
PhenomD
IMRPhenomD(Husa et al., 2016; Khan et al., 2016) is the latest in a series of phenomenological models with non-precessing spins. It uses hybrids of uncalibrated SEOBNRv2 and 19 different NR waveforms to fit the model to. The range of NR waveforms includes mass ratios up to 16, and non-precessing spins from -0.95 to 0.98 for equal mass binaries, and roughly in between ±0.8 for non-equal masses (see Table 3.1.2 or Table 1 in Khan et al. (2016)). Configuration with spin magnitudes greater than 0.85 and non-equal masses for example, are not supported by this calibration range.
The ansatz is split into two different models for the wave’s amplitude and phase, which are in turn split into three parts for the inspiral, intermediate and merger-ringdown regions of the model (Husa et al., 2016;
Khan et al., 2016). The complete IMR phase is given by: φ(f ) = φP N(M f ) +η1 σ0+ σ1f +34σ2f4/3+35σ3f5/3+12σ4f2 for M f ≤ 0.018 1 η β0+ β1f + β2log(f ) − 1 3β3f −3 for 0.018 < M f ≤ 0.5fRD 1 η α0+ α1f − α2f−1+43α3f3/4+ α4tan−1 f − αf 5fRD damp for M f > 0.5fRD , (3.3) where η = m1m2/M2 is the symmetric mass ratio, and G = c = 1. For the inspiral region, both the
phase and amplitude are based on the TaylorF2 model, a known post-Newtonian expansion in the frequency domain (see, for example, Buonanno et al. (2009)), for which the phase is denoted φP N. Some higher order
terms are included with extra phenomenological parameters {σi} as prefactors. At M f = 0.018, the example
waveforms consist of purely NR data, at which point the intermediate regime starts. For this region as well as the merger-ringdown, purely phenomenological guesses are used based on the evolution of the amplitude and phase, with coefficients {βi} for the intermediate and {αi} for the merger ringdown regime. Some of
these, namely β0, β1, α0 and α1, are used solely to ensure the continuity of the phase and its derivative.
The boundary between the two regions is at M f = 0.5fRD, where fRD is the ringdown frequency of the
final black hole which is calculated from its mass and spin. fdamp, used in modelling the ringdown, is the
damping frequency of the final black hole, which can also be calculated from the mass and spin parameters. The inspiral coefficients {σi} are expressed as polynomials in the physical parameters m1, m2, χ1and χ2.
Since this is a non-precessing model, χ1and χ2are considered as the spin components parallel to the orbital
angular momentum. In the fit of the ansatz to the hybrid waveforms, the coefficients of the polynomials are determined. (The exact expressions can be found in appendix B of Khan et al. (2016).) This allows for the coefficients, and thus the complete model, to be constructed with any choice of physical parameters. The intermediate and merger-ringdown coefficients are also written as polynomials, but only using the physical parameters η (the symmetric mass ratio) and χP N, a combination of the spins of the two objects. χP N is
chosen as it determines the leading-order spin effect in the phase of the gravitational wave. Using it instead of the two spins is a good approximation for most binaries, but may not hold for high mass ratios or high spins. Its value is given by (see also Eq. (2) and (3) in Khan et al. (2016)):
χe= m1χ1+ m2χ2 M , (3.4) χP N = χe− 38η 113(χ1+ χ2). (3.5) Effective precession
The IMRPhenomD model is only applicable to binaries with non-precessing spins, i.e. spins aligned with the orbital angular momentum ~L. A full spin phenomenological model would require using four additional physical parameters (two additional spin components for each object) in the ansatz, and many more NR waveforms to cover the extended parameter space. However, a method has been developed by Hannam et al. (2014) to turn a non-precessing into a precessing waveform model; no additional calibration with precessing waveforms is used. The resulting model for precessing spin configuration is called IMRPhenomPv2.
In the method by Hannam et al. (2014), the complete inspiral-merger-ringdown waveform is calculated first using only the spin components aligned with the orbital angular momentum, a valid approximation if the in-plane spin components do not have a significant effect on the inspiral. This is true as long as the direction of the total angular momentum ~J does not change (transitional precession), which is expected to happen only in rare cases (Schmidt et al., 2015). Then, the precession of ~L around ~J is modelled using the in-plane spin components, and used to find a rotation that transforms the non-precessing waveform into a precessing one.
In general, this rotation is defined by two angles that depend on all four in-plane spin components and the binary’s mass ratio. They can be calculated using leading-order (or higher) PN equations for the precession. For a generic spin configuration, there are no analytic solutions to these equations, but for the special case
nr. q η χ1 χ2 χˆ A1 1.0 0.25 -0.95 -0.95 -0.95 A2 1.0 0.25 -0.6 -0.6 -0.6 A3 1.0 0.25 0.0 0.0 0.0 A4 1.0 0.25 0.6 0.6 0.6 A5 1.0 0.25 0.98 0.98 0.98 A6 4.0 0.16 -0.75 -0.75 -0.75 A7 4.0 0.16 -0.5 -0.5 -0.5 A8 4.0 0.16 0.0 0.0 0.0 A9 4.0 0.16 0.5 0.5 0.5 A10 4.0 0.16 0.75 0.75 0.75 A11 8.0 0.099 -0.85 -0.85 -0.85 A12 8.0 0.099 -0.5 0.0 -0.458 A13 8.0 0.099 0.0 0.0 0.0 A14 8.0 0.099 0.5 0.0 0.458 A15 8.0 0.099 0.85 0.85 0.85 A16 18.0 0.05 -0.8 0.0 -0.77 A17 18.0 0.05 -0.4 0.0 -0.385 A18 18.0 0.05 0.0 0.0 0.0 A19 18.0 0.05 0.4 0.0 0.385
Table 3.2: NR waveforms used in the calibration of the PhenomD model. Adapted from Khan et al. (2016), using the same identification numbers. Given are the mass ratio q, the symmetric mass ratio η, both dimensionless spins χ1and χ2, and ˆχ, which is χP N normalised between -1 and 1 (see Eq. (4) in Khan et al.
(2016)).
where one of the black holes has zero in-plane spin, there are. In the model by Hannam et al. (2014), the four spin components are reduced to one effective precession spin parameter χp. A single spin configuration
is obtained by assigning χp as the spin of the larger black hole, which serves as an approximation of the
full configuration. The motivation for χp is presented here, following Schmidt et al. (2015). It starts by
considering the leading-order precession effect, namely spin-orbit coupling at 1.5PN order: ˙
~ L = L
r3A1S~1⊥+ A2S~2⊥ × ˆL. (3.6)
Here, r is the separation and A1= 2+2q3, A2= 2+3q2. (Although Schmidt et al. (2015) define q ≡ m2/m1> 0,
the equations are adapted to the convention of q ≡ m1/m2> 0 used throughout this work.)
As the spins both precess at, in general, different rates around ~L, their relative angle changes continuously. Therefore, the term in square brackets in Eq. (3.6) varies between its minimum when the two spins points in opposite directions and its maximum when the two spins are aligned. An approximation is made for the average of the spin-term, given by the average of the minimum and maximum values:
Sp≡
1
2(A1S1⊥+ A2S2⊥+ |A1S1⊥− A2S2⊥|) = max(A1S1⊥, A2S2⊥). (3.7) The spin Sp can either be determined by the spin of the larger or the smaller black hole. Usually, S1⊥
dominates the precession, so in this approximation all spin is assigned to the larger black hole. This then defines the effective spin parameter:
χp≡
Sp
A1m21
. (3.8)
In some cases, the representation of two in-plane spin components with χp is not a good approximation.
spins precess at different rates around ~L, such that the weighted sum of the two in-plane spin components could be approximated by Sp as in Eq. (3.7). The rate of precession depends on the mass of the black hole,
thus the rates are (almost) the same for (almost) equal mass systems. In such a case, the relative angle between the two in-plane spins will be constant. This determines the size of A1S~1⊥+ A2S~2⊥, and therefore
the strength of the spin-orbit coupling. For example, for equal masses and an in-plane angle φ = 0, the spin components add directly and the precession effect will be underestimated by χp. If on the other hand φ = π,
the spin components cancel out and there is no spin-orbit coupling at all. Schmidt et al. (2015) found that this becomes a problem for mass ratios q . 1.2, although the effect is comparable to that of neglecting the higher order spin-spin coupling.
The second assumption was that the larger black hole will dominate the spin-orbit coupling, which is why χp is assigned as its spin. However, this is not the case if A2S2⊥ > A1S1⊥. Using χi ≡ Si/m2i, this
condition can be rewritten as:
χ1⊥<
A2
A1q2
χ2⊥. (3.9)
For q ' 1.0, this is simply when the in-plane spin of the smaller black hole is larger than that of the bigger black hole. The higher the mass ratio, the greater the difference between the spins must be for the con-dition to be met. (For large q, the factor A2/(A1q2) is proportional to 1/q.) Schmidt et al. (2015) note
that although one of the precession angles - the one between ~L and ~J - is affected by this problem, only a small percentage of the resulting waveforms deviate significantly from waveforms constructed with the full spin configuration. The problems occurs mostly for configurations where the in-plane spin component of the larger black hole is small, χ1≤ 0.08.
3.2
Overlaps
In order for a waveform model to do its job, it needs to be both fast and approximate the actual GR solution accurately. The parameterised test of GR requires the calculation of so many instances of the waveform, that only IMRPhenomPv2 and SEOBNRv2ROM meet the speed requirement. As the latter admits no precessing spins, IMRPhenomPv2 is the usual approximant of choice for these kinds of searches. It is therefore important that its accuracy is verified.
The most direct way of checking a GR approximation is comparing the waveform to NR results. But NR waveforms only exist for a select set of binary parameters. A test that allows any choice of parameters is comparing two waveform models to each other. If they both accurately represent the GR solution at that point, they should be very similar. Such a test is discussed in this section, in which the overlap between IMRPhenomPv2and SEOBNRv2ROM or SEOBNRv3 is examined. Note that a similar test for IMRPhenomD and SEOBNRv2 has already been done by Khan et al. (2016), for a number of different masses and mass ratios with equal spin configurations (see also the discussion of the results 3.2.4). Here, the test will be extended to non-equal spins and precessing spins for the waveforms IMRPhenomPv2 and SEOBNRv3.
3.2.1
Faithfulness
In order to compare two waveforms, a quantitative measure for their similarity is required. As is used in the calculation of the evidence integral (see section 2.2.1 for a definition of the symbols), the inner product between to waveforms is defined as:
(h|g) = 4Re Z ∞ 0 ˜ h∗(f ) ˜g(f ) Sn(f ) df. (3.10)
Sn(f ) is the noise power spectral density (PSD), i.e. the power of the noise per frequency bin df , as a function
or a fit based on measured noise. In this work, an estimate of the noise curve for the early stages of Advanced LIGO (Abbott et al., 2013) is interpolated and the integral is cut off from 20 Hz to the lowest maximum frequency of the two waveforms that are used. A plot of the noise PSD is given in Fig. 3.2.1.
For the product to be of practical use as a similarity measure, the waveforms are normalised such that (h|h) = 1, ensuring that the product of two different waveforms always lies in the range of minus one to one, with negative values for anti-correlated waveforms. Each waveform can be given an arbitrary start time t0 and an arbitrary phase at this time, φ0. It makes sense to maximise over these two parameters when
matching two waveforms (either a model and a signal or two models). Applying this to Eq. (3.10) yields the faithfulness:
F (h, g) = max
t0,φ0
(h(λ)|g(λ))
p(h(λ)|h(λ)p(g(λ)|g(λ). (3.11) Here, the faithfulness F is evaluated for the physical parameters λ. When matching a waveform model to a signal with unknown parameters λ0, the match is maximised over λ as well. The overlap measure that is obtained in this way is called the effectualness:
E(h, g) = max
t0,φ0,λ
(h(λ0)|g(λ))
p(h(λ0)|h(λ0)p(g(λ)|g(λ). (3.12)
A waveform model that is effectual does not have to be faithful as well; it could be that the maximum match is obtained for model parameters λ that deviate significantly from the signal parameters λ0. It would be sufficient for the detection of a signal, but not suitable for parameter estimation. As we are concerned with the accuracy of waveform models for given parameters, the measure of choice is the faithfulness.
Figure 3.1: Predicted noise curve for the early stages of the Advanced LIGO detectors. The red dashed vertical line indicates the lower frequency cut-off used in calculating Faithfulness values.
The faithfulness calculation is done numerically, and therefore, the integral (3.10) is replaced with a sum. The waveforms and noise PSD are discrete data sets of equal length, which is managed by cutting them off at
the same minimum frequency of 20 Hz and the same maximum frequency, determined by whichever one of the waveforms has the lowest. The waveforms are generated with code built into LALSuite, and only one of the polarisation, h+, is used. The h× polarisation is the same as h+ with possibly a different overall amplitude
(depending on the inclination of the binary), which would be lost in normalising the waveform. It therefore does not contain any additional information for the overlap. Since the SEOBNRv3 code produces a time domain waveform, its output is transformed with a fast Fourier transform (FFT). Before the FFT, a tapered window is added to the waveform at frequencies below 20 Hz, to avoid sharp cut-offs. The SEOBNRv2ROM and IMRPhenomPv2 waveforms are directly generated in the frequency domain.
For the maximisation of the product, both waveforms are generated once with the given parameters. Afterwards, the external parameters φ0 and t0 are changed by multiplying the waveform with a factor
e2iφ0 and an array [e−2iπf t0] for each frequency f , respectively. For the optimisation, the python routine
optimize.fminfrom the Scipy module is used (a minus sign is introduced in the faithfulness function so that it can be minimised instead of maximised).
3.2.2
Probing parameter space
For the parameterised test of GR on GW150914 and GW151226, the waveform IMRPhenomPv2 has been used. In order to back up the GR tests, IMRPhenomPv2 would only need to be checked around the binary parameters of both events. The total masses and mass ratios are quite well-measured (M = 64.8+4.6−3.9 M ,
q = 1.22+0.42−0.20and M = 21.8+5.9−1.7M , q = 1.92+2.43−0.83 respectively, (Abbott et al., 2016c)), but this is not true
for most other parameters. Besides this limitation, a test that covers a larger part of parameter space is also more desirable than a more restricted test. It would provide extended information on the waveforms themselves, an also possibly be of use when new gravitational waves are measured. In the most general case, the intrinsic parameter space of coalescing binaries is eight dimensional: two masses and two spins vectors with three components each. The position and the orientation of the system relative to the observer add another four parameters. It is impossible to vary all of these parameters at the same time, so a feasible test will have to be restricted in some way. It is therefore useful to investigate which parameters have any, or the most significant, effect on the waveform.
For the extrinsic parameters, only the inclination of the binary system can affect the shape of the waveform. The sky position of the source determines the angle at which the wave passes through the detector but this only affects the relative strength of the two waveform polarisations, not their shape. Only the latter is considered in this work, hence a single polarisation is used and the sky position is irrelevant. This is also true for the distance, which is responsible for an overall factor in the amplitude. This information is lost in the faithfulness calculation, since the waveforms are normalised. The inclination only has an effect on waves from binaries with precessing spins; the amplitude and phase modulations caused by precession are weakest for face-on and strongest for edge-on systems Schmidt et al. (2015).
On the other hand, most of the intrinsic parameters of the binary are important for the shape of the waveform. For non-precessing spins, there are four: m1, m2, χ1 and χ2, where both spins are of course
aligned with the orbital angular momentum. An equivalent description of the two component masses is a combination of the total mass M and the mass ratio q, which is used in this work. Precessing spins need three parameters each, which are chosen to be the magnitude χi, the angle with the orbital angular momentum θi
and the angle of the in-plane projection of the spin with the x-axis φi (for an overview of the angles involved,
see Fig. 3.2). The first two are clearly relevant parameters; their combination determines the size of the spin components aligned with ~L and perpendicular to ~L. In the phenomenological models, these are reduced to two parameters χefrom the aligned and χp from the perpendicular components, and appear to be the final
measurable parameters (Khan et al., 2016; Schmidt et al., 2015). However, these are approximations and in general all four spin components are relevant.
On the other hand, both in-plane angles φi may not be relevant at all. The black holes orbit each other,
and there is no fixed start time (it is maximised over in the faithfulness calculation), so the direction of the x-axis can be chosen arbitrarily. Then, both angles φ1 and φ2 reduce to their difference, φ = |φ1− φ2|.
Figure 3.2: Overview of the angles involved with the spin and angular momentum vectors. ~L is the orbital angular momentum of the binary, ~S1 the spin of the bigger black hole and ~S2 that of the smaller. Vector
sizes are chosen arbitrarily in this figure.
that this is what the χp approximation is built upon (see section 3.1.2). If the spins rotate at the same
rate however - which is the case for equal mass equal spins systems - the initial value for φ is maintained. Only in those configurations, φ is relevant as is determines the strength of the main precession effect, namely spin-orbit coupling at 1.5PN order (see Eq. (3.6)).
Identifying all the possibly relevant parameters yields the following list: M , q, χ1, χ2, θ1, θ2, φ and ι
(see also Table 3.2.2). Even for non-precessing spins, which require only the first four parameters, this is too much to vary over a large range at the same time, let alone visualise in a two-dimensional plot. Focusing on the parameters of GW150914 or GW151226 allows for a choice of the masses and mass ratios. Respecting the uncertainty in these estimated parameters, the values at the upper and lower 90% confidence intervals (CI) are used as well (Abbott et al., 2016a). For GW150914, an extra mass ratio value of 1.43 is added to fill the gap between the estimated value of 1.22 and the upper 90% CI of 1.64. A similar thing is done for GW151226, but for both the upper and the lower 90% CI since the uncertainties are rather large.
Having a few fixed values for the masses, the spin magnitudes remain as the only relevant variables for the non-precessing case. The dimensionless spins χi= | ~Si|/m2i have a physically limited range from zero to
one, so they are easy to use as variables. To use both directions of the spin (one with a component aligned with ~L and one with a component anti-aligned with ~L), the range is extended from minus one to one. χ1 is
chosen as the x-variable and χ2as the y-variable, and both ranges of [-1, 1] are probed at 50 points, resulting
in a total of 50 × 50 = 2500 faithfulness calculations per test. These points can be visualised as a colour plot with χ1 and χ2 on the x- and y-axes, respectively.
The same is done for the precessing spin configuration, where a choice needs to be made for the remaining parameters θ1, θ2, φ and ι. The faithfulness is a somewhat computationally expensive calculation, since both
waveforms need to be calculated for each set of parameters, and the value needs to be optimised over t0
and φ0. Besides, the number of resulting plots will multiply for the number of choices of each remaining
parameter, so too many choices will clutter the results. Therefore, the choices for θ1, θ2 are restricted to
be equal, and vary over the values π/6, π/3, π/2, such that the precession effects are of different strength for each choice. Likewise, two values for the inclination are chosen, 0 and π/4, where the precession effects should be more visible for the second than for the first choice. As the φ parameter will only be relevant (at least to leading order) for equal mass, equal spin systems, it is less interesting to vary this parameter and so only the value φ = π is used, meaning that the two in-plane components of the spins point in opposite direction at the start time t0. All discussed choices for parameters are summarised in Table 3.2.2.
parameter non-precessing spins precessing spins values for GW150914 test values for GW151226 test M - total mass (M ) yes yes 60.9, 64.8, 69.4 20.1, 21.8, 27.7
q - mass ratio (m1/m2 > 0) yes yes 1.0, 1.22, 1.43, 1.64 1.09, 1.51, 1.92,
3.14, 4.35
χ1 - magnitude spin 1 yes yes x-variable: [-1, 1] x-variable: [-1, 1]
χ2 - magnitude spin 2 yes yes y-variable: [-1, 1] y-variable: [-1, 1]
θ1 - angle with ~L of spin 1 no yes 0, π/6, π/3, π/2 0,π/6, π/3, π/2
θ2 - angle with ~L of spin 2 no yes 0, π/6, π/3, π/2 0,π/6, π/3, π/2
φ - in-plane spin angle no in some cases π π ι - inclination no yes 0, = π/4 0,π/4
Table 3.3: Parameters used in the overlaps test of IMRPhenomPv2. The second and third columns indicate whether the parameter is relevant for respectively a configuration with non-precessing or precessing spins. The fourth and fifth columns list the values for each parameter that have been used in the test. For the total mass and mass ratio, the estimated median values for GW150914 and GW151226 are in bold font, while the other values are based on the 90% confidence intervals (Abbott et al., 2016c,a). Since no final results for the precessing spin configurations for GW151226 were obtained, these values are shown in gray.
3.2.3
Results
Non-precessing spins
Both the SEOBNRv2ROM and the SEOBNRv3 model can produce waveforms from non-precessing spin con-figurations, so both are used in a comparison with IMRPhenomPv2 (the ROM version of SEOBNRv2 is used since it is faster). First, the results of the studies with GW150914 parameters are shown in Fig. 3.3 and 3.4, followed by the ones for GW151226 in Fig. 3.5 and 3.6. In each study, all combinations of choices for M and q (given in Table 3.2.2) are used, resulting in twelve plots for GW150914 and fifteen for GW151226. The plots are ordered by increasing M from left to right and increasing q from bottom to top.
Precessing spins
Overlap studies with precessing spin configuration were performed using the waveforms IMRPhenomPv2and SEOBNRv3. No final results were obtained for the studies with mass parameters of GW151226, so only results for GW150914 are presented here. Studies with an orbital inclination of zero are shown first, with angles between the spin components and ~L of π/6 in Fig. 3.7, π/3 in Fig. 3.8 and π/2 in Fig. 3.9. Results for configurations with the same precession angles but with an inclination of π/4 are shown in Fig. 3.10, 3.11 and 3.12. The individual plots with different values for M and q are ordered in the same way as for the non-precessing spin studies.
3.2.4
Discussion
Meaning of the Faithfulness scale
In the results in general, faithfulness values vary drastically from close to 1.0 to low values of around 0.7. In order to interpret these results, it is useful to provide some context. Concerning the accuracy of the parameterised test (or other data analysis purposes), what finally matters is by how much the inaccuracy of the waveform affects the result, compared to for example the error due to noise. Abbott et al. (2016d) report, besides the parameterised test of GR on GW150914, properties of the residuals of the data after subtracting the highest likelihood model waveform. An upper bound on the signal to noise ratio (SNR) in the residuals of SNRres≤ 7.3 (at 95% confidence) was found, consistent with SNR found in noise at times close to the signal.
