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

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.

Add

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

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.

(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

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.

Comment on your result [10

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))[0].

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[0],preddur[1]))

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

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

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]

Let’s block ads! (Why?)

Send us your paper details now

We’ll find the best professional writer for you!

This template supports the sidebar's widgets. Add one or use Full Width layout.

This website uses cookies to improve your experience. We'll assume you're ok with this, but you can opt-out if you wish.Accept Read More

Privacy & Cookies Policy

error: Content is protected !!