disturbances in the curvature of spacetime

Gravitational waves are disturbances in the curvature of spacetime, generated by accelerated masses, that propagate as waves outward from their source at the speed of…

waves are disturbances in the curvature of spacetime, generated by
accelerated masses, that propagate as waves outward from their source
at the speed of light. They are predicted in General Relativity and
other theories of gravity and since 2017, they have now been

this exercise we will analyse some mock gravitational wave data from
two unknown astrophysical objects merging together and coelescing. We
will use a Monte Carlo Markov Chain (MCMC) to compare a scaled model
that predicts how the wave changes depending on the total mass of the
merging objects and their distance from us to the observed waveform.
This will allow us to determine the nature of the orbiting objects
that merged to form the gravitational wave using MCMC, whether for
instance they could be originating from merging white dwarfs, neutron
stars or black holes.

mock or simulated waveforms measure the strain as two compact, dense
astrophysical objects coalesce. The strain describes the amplitude of
the wave. The system is parameterised by the masses of the merging
objects, π‘€1M1 and π‘€2M2,
and their distance from the observer π·D.

useful parameters are:

mass ratio π‘ž=𝑀2/𝑀1q=M2/M1,
with convention that π‘€1β‰₯𝑀2M1β‰₯M2 and
so π‘žβ‰€1q≀1.

β€œChirp mass”, which is a quantity used in general
relativity, is given by:



all your reasoning for each step. A significant fraction of the
marks are given for explanations and discussion, as they evidence
understanding of the analysis.

of these steps will take a while to run and compile. It’s a good
idea to add in print statements to your code throughout
eg print(β€˜this
step is done’) to
make sure that your bit of code has finished.

the import packages statements in the cell below to the top of your
Jupyter notebook. We will use the pandas package
to read in the data, with eg dataIn=pd.read_csv(β€˜filename.csv’).

may find it useful to look at the following publication from the
LIGO consortium. https://arxiv.org/pdf/1608.01940.pdf

A- Some background [20 marks]

(i) How
do astronomers know that the gravitational waves from the GW150914
event was due to two black holes merging? [5

(ii) Describe
the different parts of the waveform produced due to gravitational
waves from a merging black hole event. [5

(iii) The
file CA3_gravitationalwaveevents.csv contains the properties of
previously observed gravitational waves. Parameters obtained for
these systems are known as posteriors as they are derived using
Bayes Theorem (ie posterior = likelihood Γ—Γ— prior).
The errors in these values are the credible intervals from the
posterior distribution derived for each parameter. Plot
the total mass
of the merging sources against their distance. [10

B – The data [15 marks]

it is your turn to look at observations and measure the mass and
distance of the merging system. We first need to access the
observational data measured with the gravitational wave detectors
(the waveform observed when two compact, dense astrophysical objects
coalesce), and format it correctly.

(i) Read
in the datafile of the observed waveform Dataset1.csv.
These files store the strain as a function of β€œGPS time”
for the merger of two bodies. [5

(ii) The
GPS time of the merger for your waveform is 1205951542.153363. Your
data will need to be shifted so that the merger occurs at time = 0
secs. This is required for when we compare model waveforms with our
data. [4

final data step is that we need to estimate the average noise and
standard deviation in our data. This requires careful thought about
where the noise can be seen in the waveform. [6

C – Using model waveforms to roughly estimate the total mass of two
objects and distance to the system [28 marks]

that we have some data, we ideally want to try and compare predicted
waveforms for a merger with our observed waveform and use this to
derive the mass and distance to the merger.

waveforms are known as templates in
the gravitational wave community. In this part of the question we
will attempt to produce a waveform for different masses and distances
using a reference waveform and scaling it by its change in distance
and mass.

reference waveform/template we will use
This file shows the strain as a function of time (with π‘‘=0t=0 at
the merger) of a system with total mass of the merging objects
as π‘€=40𝑀𝑠𝑒𝑛M=40Msun, π·=1D=1Mpc
and π‘ž=𝑀2/𝑀1=1q=M2/M1=1.
The template has been generated with a lowest frequency
of π‘“low=10flow=10Hz.

we have an equal-mass system (i.e π‘ž=1q=1)
with total mass π‘€=𝑀1+𝑀2M=M1+M2 at
a distance π·D,
then we can scale the strain of a waveform, β„Ž(𝑑,𝑀,𝐷)h(t,M,D) from
a reference waveform with π‘€β€²,𝐷′Mβ€²,Dβ€² as:

β„Ž(𝑑,𝑀,𝐷)=(𝑀𝑀′)(𝐷′𝐷) β„Ž(𝑑′)         [πΈπ‘ž. 1]h(t,M,D)=(MMβ€²)(Dβ€²D) h(tβ€²)         [Eq. 1]


𝑑′=(𝑀′𝑀)𝑑         [πΈπ‘ž. 2]tβ€²=(Mβ€²M)t         [Eq. 2]

and π‘‘=0t=0 is
defined as being the time at which the merger occurred.

that the templates will have different time sampling rates to your
observed data as they are model waveforms.

will need to include the following steps when answering this

(i) Open
the template file using the pandas package.
Write a function in python to scale the time and strain of any
waveform with π‘ž=1q=1,
total mass π‘€M and
distance π·D from
a reference waveform using Equations~1 and 2. In this work we will
use the parameters π‘€=40𝑀𝑠𝑒𝑛M=40Msun and π·=1D=1Mpc
for our reference/template waveform. [10

(ii) Test
your function from part (i) works by subsituting in
parameters π‘€=70𝑀𝑠𝑒𝑛M=70Msun and π·=5D=5Mpc,
and comparing the waveform you produce using your python function
with the model template file generated
for π‘€=70𝑀𝑠𝑒𝑛M=70Msun and π·=5D=5Mpc: reference_Mtot70.0Msun_q1.00_Dist5.0Mpc_f_lower10Hz.csv.
Comment on your result [10

(iii) Use
your function to scale the template waveform
(𝑀=40𝑀𝑠𝑒𝑛M=40Msun, π·=1D=1Mpc
and π‘ž=𝑀2/𝑀1=1q=M2/M1=1)
in mass and distance to make an initial rough estimate β€œby eye”
of the total mass and distance that fit your data (e.g. to within
+/- 5 Msun, +/- 100 Mpc). [8

D – Interpolate the data [13 marks]

now need to make a function that will interpolate over our data
waveform and allow us to produce an observed waveform with the same
time sampling as the template waveform. This will be crucial for
comparing our data (ie our observations) with the templates (ie our
model) at each time step later on. Here we will also consider only
the time range for which we have data.

(i) Our
data waveform starts at some time π‘‘t.
Find out what this is. Next, take your observed data waveform and
only output data for π‘‘t > π‘‘min(data)tmin(data) and π‘‘t < 00 (ie
only keep information for times β‰€0≀0 (before
the merger), or for times where there is data). Verify, by plotting,
that your new observed waveform only has data in this restricted
time range. [5

(ii) OK,
final step before we get onto to estimating the total mass of the
system. To estimate the difference between our observed waveform and
model/template waveforms, we now need to put both observed and
template waveforms on the same time sampling. To do this we need to
interpolate between our observed data so that we can find an
observed strain for each time step and compare it to the predicted
strain for the same timestep from our scaled model. Verify this
works via plotting or checking the length of the original data array
and the new interpolated data array. [8


(i), one can use the following code example index
= np.where((data_old > 5)&(data_old<10))[0].
In simple terms, this type of statement returns a list of indices
where the conditions in the bracket have been met. One can then
use data_old[index] to
pull out the data_old numbers
that satisfy the conditions.

(ii), you should use interpolation, such as that created
by scipy.interpolate.interp1d.
You can use the following code to do this. If x[index] and y[index]
are the observed data from the previous steps, then interp_fn
create an interpolation object of the data and new_data_y
= interp_fn(scaled_template_time) will
take the time values from the scaled template, and output a value
for the observed data for each time step of the template.

E – Estimating the best fit total mass using MCMC [55 marks]

we will derive the best fitting total mass for the two bodies that
collided to produce the wave we see in our data using an MCMC
approach. We will take our β€œby-eye” estimate of the
distance above and then use MCMC to derive the total mass of the
system. We will be comparing the observed waveform with the reference
waveform (40𝑀sunMsun,
D=1Mpc from earlier) (the model) so think carefully about what the
likelihood function will be in this case (see Chapters 5-8).

this question, the following steps will be needed:

(i) Use
MCMC to sample the total mass, π‘€M,
to produce a best-fit value for your data. [25

(ii) Display
the results in an appropriate manner and comment on your findings,
as well as your results from the MCMC. [20

(iii) Report
the median and 90% credible limits on your value of π‘€M. [10

may assume that:

noise is white and described by a Gaussian distribution,

total mass is in the range [20,100] π‘€π‘ π‘’𝑛Msun.


compare the model template with the data, we will need to have both
the reference waveform and the observed data produced for the same
time steps and time range. We can do this by interpolating and
setting the range of time we want to compare to. An example of
interpolating in python is given
by data_interp=interp1d(data_time,data_strain,bounds_error=False)

should work with β€œlog(Likelihood)” to avoid numerical
errors – note this will affect both your posterior and the step in
the MCMC chain where we usually
write π‘proposed/𝑝currentpproposed/pcurrent

very carefully about the form of your likelihood since here we are
comparing observed data with a model when the total mass π‘€M is
changed _

sure that the calculations within each step are as quick as
possible, so perform as much computation as possible β€œoutside”
the MCMC

step size between samples of the MCMC is quite important. A
suggested value is 0.1𝑀𝑠𝑒𝑛0.1Msun

your MCMC on a small number of samples (e.g. 100) before trying it
with a larger number (e.g. 105105 or 106106)

the end, ask yourself if you need to include every sample?

on your step size, this part can take a long time to run. Move all
your plotting routines to a different code cell to save you
re-running everything 10000s of times when you just want to change a
plot command. To find out how long it will take for a Jupyter
notebook to compile the MCMC code cell, add the following snippet to
your code before you go into your MCMC loop (where Nsteps is the
number of steps your MCMC is using):

datetime import datetime,timedelta




duration: {}-{}s’.format(preddur[0],preddur[1]))


Part F –
Estimating mass and distance using MCMC for the merging system [37

your MCMC analysis to now estimate the total mass, π‘€π‘€, and the
distance, π·π·.

(i) Use
MCMC to get π‘€M and π·D. [15

(ii) Display
the results in an appropriate manner and comment on your findings,
as well as your results from the MCMC. [10

(iii) Report
the median and 90% credible limits on your value of π‘€M and
compare your best fit waveform with the data. Comment on your
result. [12

if there are any difficulties completing this component of the
coursework, you can still attempt Part G using your by-eye estimates
for π‘€M and π·D from
Part C, or your π‘€M from
the MCMC from Part E and your by eye estimate for π·D in
Part C.


your theta_current for the MCMC will have 2 values (Mtot and D) as
such you will need to specify 2 initial parameters for
theta_current. Whereas before you may have used something
like theta_current.append([Mtot_guess]),
now you will need something
like theta_current.append(np.array([Mtot_guess,D_guess])) and
similarly your D_theta will
need to have 2 values to tell the MCMC to take a random step in both
Mtot and D values in each step of the MCMC.

careful not to get caught up in messy arrays – if in doubt run the
MCMC for 5-10 steps and print out the theta_current to check your
theta_current[0] ends up looking like it should (ie first element =
initial parameters you set up, next elements either the same as
before or changed to proposed theta values). If in doubt, ask for
help at the weekly exercise classes.

G – Putting it all together [32 marks]

have now measured the total mass for your colliding astrophysical
objects and the distance.

(i) Calculate
the chirp mass for your system and the individual masses of your
merging bodies. Comment on your individual masses. [5

(ii) Plot
your MCMC derived properties alongside the previously discovered
gravitational wave systems. [5

(iii) Assuming
that the objects are not spinning, and that their orbits are
Keplerian and essentially circular before the peak amplitude in the
wave is reached, estimate the orbital separation in km of the two
bodies around peak amplitude. Think carefully about how the orbital
period is related to the period of your gravitational wave. [18

(iv) Comment
on what your analysis suggests are the best astrophysical candidates
for the merging objects? [4

Let’s block ads! (Why?)

Do you need any assistance with this question?
Send us your paper details now
We’ll find the best professional writer for you!


error: Content is protected !!