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…

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 light. They are predicted in General Relativity and
other theories of gravity and since 2017, they have now been
observed!

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

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

Other
useful parameters are:

The
mass ratio 𝑞=𝑀2/𝑀1q=M2/M1,
with convention that 𝑀1≥𝑀2M1≥M2 and
so 𝑞≤1q≤1.

The
“Chirp mass”, which is a quantity used in general
relativity, is given by:

𝑀𝑐ℎ=(𝑀1𝑀2)3/5(𝑀1+𝑀2)1/5Mch=(M1M2)3/5(M1+M2)1/5

General
tips:

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

Some
of these steps will take a while to run and compile. It’s a good
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

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

Part
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
marks]

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

(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
marks]

Part
B – The data [15 marks]

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

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
marks]

(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
marks]

(iii)The
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
marks]

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

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

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

The
reference waveform/template we will use
isreference_Mtot40.0Msun_q1.00_Dist1.0Mpc_f_lower10Hz.csv.
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.

If
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]

where:

𝑡′=(𝑀′𝑀)𝑡         [𝐸𝑞. 2]t′=(M′M)t         [Eq. 2]

and 𝑡=0t=0 is
defined as being the time at which the merger occurred.

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

You
will need to include the following steps when answering this
question:

(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
marks]

(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.
marks]

(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
marks]

Part
D – Interpolate the data [13 marks]

We
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
marks]

(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
marks]

Hints:

For
(i), one can use the following code example index
= np.where((data_old > 5)&(data_old<10)).
In simple terms, this type of statement returns a list of indices
(index)
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.

For
(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
=interp1d(x[index],y[index],bounds_error=False)will
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.

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

Here
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).

For
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
marks]

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

(iii) Report
the median and 90% credible limits on your value of 𝑀M. [10
marks]

You
may assume that:

the
noise is white and described by a Gaussian distribution,

the
total mass is in the range [20,100] 𝑀𝑠𝑢𝑛Msun.

Hints:

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

You
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

_Think
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 _

make
sure that the calculations within each step are as quick as
possible, so perform as much computation as possible “outside”
the MCMC

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

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

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

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

from
datetime import datetime,timedelta

tstart=datetime.now()

print(‘start
time:’,tstart)

preddur=[Nsteps*0.01,Nsteps*0.02]

print(‘predicted
duration: {}-{}s’.format(preddur,preddur))

predtend=[tstart+timedelta(seconds=preddur),tstart+timedelta(seconds=preddur)]

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

Repeat
your MCMC analysis to now estimate the total mass, 𝑀𝑀, and the
distance, 𝐷𝐷.

(i) Use
MCMC to get 𝑀M and 𝐷D. [15
marks]

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

(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
marks]

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

Hints:

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

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

Part
G – Putting it all together [32 marks]

You
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
marks]

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

(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
marks]

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