Miles Oleksak & Alvaro Feito Boirac

In 1970, Physicist Harold D. Craft, Jr. published his PhD thesis. It included a novel visualization of signals arriving from neutron stars.

The central image was the inspiration for designer Peter Saville's now iconic album cover:

Inspired by the art and science, we decided to give an update to that iconic image. The pulsars (or neutron stars) that star in that album cover were first measured half a century ago by Jocelyn Bell. Could we display more recent discoveries or data in a similar fashion? One of the most astonishing experiments in astrophysics is the detection of gravitational waves. They were first observed 5 years ago and were so ground-breaking that the team immediately got the Nobel prize the following year.

The following image shows data measured in the LIGO observatory between 2015 and 2019. In this article we will try to explain the meaning of this graph and the Physics behind it.

ðŸ‘• *If you would like to support initiatives like these, order a Wave Divisions T-shirt!*

The image above shows data coming from *Gravitational Waves*.
But what are they? This discussion will necessarily ommit
most of the details. Our intention is to give a flavour of the type of Physics
which describe gravitational waves.

Image taken from LIGO's website.

Gravitational waves affect distances. When they pass through the Earth,
they shrink some distances and lengthen others. Many words in that sentence
contain subtleties (*affect*, *distance*, *shrink*), but let us keep it simple for now.

A way to measure distance is to send light down a tunnel, reflect it on a mirror, and receive it back. In the image above there are two perpendicular tunnels. Most of the time, light takes the same time on either path. However, when gravitational waves stretch one arm, light arrives slightly later than on the other one.

The *shrinking* and the delay, however, are minuscule. The tunnels are 4 km long (2.5 miles),
so the return trip takes
$t = \frac{d}{c} \approx \frac{2 \times 4 \times 10^3}{3 \times 10^8} \approx 0.00003s $.

A rough estimate tells us that when the *short* tunnel measures 4000m, the *long*
or *stretched* tunnel measures 4000.000000000000000001m. Or just slightly longer.
Where *slightly* is probably the biggest understatement **ever**!

As you can imagine, a stopwatch won't tell the difference between 0.00003s and 0.0000299999999999999999.

The measurement is carried out by *interferometry*,
the process of purposefully superimposing (overlaying, or combining)
waves to determine information about them.

In the above image of LIGO's Hanford Interferometer, you can see an aerial view of its ultra high vacuum tunnels. Each 'arm' of the L-shape experiment has a perfectly straight, single colour and coherent light beam.

The function of each arm in the overall interferometer is illustrated below.

Diagram taken from LIGO's website.

In

To summarize, LIGO aims to detect gravitational waves by

The following plots show gravitational wave data from LIGO (the same data we used to make *Wave Divisions*).

Graphs lifted from Abbott et al.

Note that the y-axis on these plots is called strain.

With $L$ the initial length, and $\Delta L$ the change in said length. In our case, it's the relative change in one arm of the LIGO intereferometer compared to the other. Thus, strain is unitless.

Typical LIGO strains are in the range $10^{-19} - 20^{-23}$.
This means that if you had a $100$km metal rod ($10^5$m), you would need to measure differences
in length which are $10^{-20}$. In this example $L = 10^5$, so $\Delta L \propto 10^{-15}$m.
However, atomic nucleai are roughly $10^{-15}$ long. In other words, if we had a solid rod from New York to Philadelphia, a
gravitational wave would change its length by one or two nuclei (not atoms, **nuclei!**).
We probably don't need to explain that such a rod is impossible to build, and that thermal expansion,
atomic vibrations, etc. would completely defeat the purpose.
It is for this reason LIGO doesn't measure distance directly, instead opting for laser
interferometry as
described above.

LIGO operates two gravitational wave observatories: one in Livingston, Louisiana
(signified by **L1** in our code) and one in Hanford, in Washington (signified by **H1**).
These sites are 3,002 kilometers (1,865 miles) apart when measured straight through the earth,
and just over this value when measured over the surface. Since gravitational waves are expected
to travel at the speed of light, this distance means we would expect a difference in the arrival
times of any incident waves of up to 10 milliseconds. Physicists working to observe gravitational waves
also consider data from a third, similar instrument located in Europe, Virgo (**V1**).

Through the use of trilateration, the difference
in arrival times can help determine the source of the wave. In the first observed gravitational wave
(event **GW150914** in the code below), this source was identified as
the merging of two distant black holes: hence **Unknown Mergers**.

The above simulation by SXS gives a visual idea of what such a merging may look like.

The following videos provide some very helpful background information and further visualisations of LIGO's methodology - we encourage the curious reader to give them a watch before reading on!

The code for all stages of generating the final plot - downloading, naming, whitening, plotting - was written in Python, using jupyter notebooks (and Deepnote for collaboration).

`gwpy`

for Gravitational Wave analysis`pandas`

&`numpy`

to process the data`matplotlib`

for all graphs`signal`

from`scipy`

for subsampling.

The Neutron Star and Black Hole collision data is available from LIGO and Virgo's online
Event List.
To that list we appended a couple of *off catalogue* but notable events from the O3 Discovery Papers.

In [ ]:

```
fig, ax = plt.subplots(6, figsize=(15,20))
fig.tight_layout(pad=1.0)
ind = 0
for ev in full_event_list[:6]:
fn = "./ligodata/{}-{}-{}.hdf5".format(ev[1], ev[2], ev[3])
strain = TimeSeries.read(fn,format='hdf5.losc' )
center = int(ev[2])
strain = strain.crop(center-15, center+15)
ax[ind].plot(strain)
ax[ind].title.set_text(ev[1] + ' - ' + ev[3]) # title is event & lab
ind = ind + 1
```

Each of LIGO's events contained $32$ seconds of strain data at $4$KHz. This means that for each of the graphs below, the number of points totals $32 \times 4,000 = 128,000$.

Filtering LIGO's data was undoubtedly necessary - the raw data and filtered data look nothing alike.

The first graph is the raw time series, and the second graph is the filtered one:

So how on earth did we get from the first image to the second image? Notice that the second graph has an amplitude of $10^{-22}$ and the first one $10^{-19}$. This means that the real signal is hiding inside the first curve, in the tiny wiggles shown by the red arrow.

The instruments used to collect the data are *ultra* sensitive; someone sneezing in the lab
could have an effect on the plots. In the first image, the large, wavy patterns
that aren't present in the second one are very likely caused by noise.

According to LIGO,

from cars driving on nearby roads to waves crashing on distant ocean shores, quantum fluctuations within the laser itself, nanometer-scale changes in the shapes of optics, even molecules crossing the path of the laser

can add noise to the signal. Hidden under these bigger
patterns are the 'true' ones: those caused by the gravitational waves, travelling through Earth from far,
far out in space. This is why it is critical to have 3 measurements on 3 separate locations on Earth.

If a signal is only present in one of them, it probably did not come from space, but was caused by local noise.
If the signal is present in the 3 locations with the delay expected by a signal travelling at the speed of light ...
you are probably onto something. But how did we eliminate the false patterns?

As we experimented with different data cleaning strategies, our 'de-noisifying' function took several forms.

In our final strategy, we combined a **bandpass filter** and a few **notch filters**.

The idea is that any signal can be expressed
as a sum of cosines.

Our own LIGO signal can be decomposed into a sum of cosines, each with an amplitude and a frequency.
The frequency inside the cosine can suggest where it comes from.

For example, the electric grid in the US
has a frequency of 60Hz which is likely to add noise to our data.

Our signal, decomposed as a sum of cosines will have a term that shows this.
If we remove the cosine with the 60Hz frequency, and add the cosines again,
our signal will look *cleaner*. Equally, we can *knock out* frequencies
which are too high or too low for gravitational waves.

To plot the de-noisified data to check it for cleanness, we also defined `plot_strain()`

.
This both returned plots of the filtered data and extracted their time-axis cut-off values; we would
later use these to select the interesting sections of each event.

Both functions were developed following these instructions
by Duncan Macleod (author of the GWpy package). Here is our `clean_it_up()`

function:

In [ ]:

```
def clean_it_up(event_index, full_event_list, low_bandpass, high_bandpass):
'''
:event_index int: chose event_index (a number between 0 and full_event_list)
:full_event_list list: defined further up this notebook
:low_bandpass int: lower cut for frequencies (typically 30Hz - 50Hz)
:high_bandpass int: higher cut for frequencies (typically 200Hz - 300Hz)
:grid_frequency int: frequency of electrical grid (which induces noise in signal), 50Hz Europe, 60Hz US
:return: filt_strain (a filtered TimeSeries object)
'''
ev = full_event_list[event_index]
# chose file
fn = "ligodata/{}-{}-{}.hdf5".format(ev[1], ev[2], ev[3])
# get strain from file
strain = TimeSeries.read(fn,format='hdf5.losc' )
# bandpass
bp = filter_design.bandpass(low_bandpass, high_bandpass, strain.sample_rate)
# remove electrical noise 60Hz
lab = ev[3]
if lab == 'H1' or lab == 'L1': # US lab
notches = [filter_design.notch(line, strain.sample_rate) for \
line in (59, 59.5, 60, 60.5, 61, 120, 180, 240)]
if lab == 'V1': # EU Lab
notches = [filter_design.notch(line, strain.sample_rate) for \
line in (49, 49.5, 50, 50.5, 51, 98, 99, 100,101,102, 148, 149, 150, 151, 152, 200)]
zpk = filter_design.concatenate_zpks(bp, *notches)
filt_strain = strain.filter(zpk, filtfilt=True)
filt_strain = filt_strain.crop(*filt_strain.span.contract(1))
return filt_strain
```

However, because each event may have had different high/low frequency filters - and may not have been exactly centered around the "center" time - we had to manually investigate and filter each of them, one by one, to obtain the best looking possible plot each time.

Inspecting plots like the above helped us determine frequencies of interest and the approximate timing of the event.

The brighter-coloured region represents the frequencies with highest amplitudes in the data (post-filtering),
so by drawing lines across
the graph to find the approximate ranges of such regions, we could safely exclude all frequencies using the bandpass
filter. This was done, quite simply, by plugging the limits of each range into `clean_it_up()`

.

We could also gauge approximate timings by repeating the same process on the x-axis, which would yield the times the data's most important frequencies occured at.

From these approximate timings, we could manually experiment with and select the exact intervals
where each graph's event occured. Once the event had been located and positioned approximately
in the center of its section, we could apply a crop to isolate the chosen interval and discard the noisy data
on either side. The sections were then appended to the new time series `the_interesting_data`

.

Several features in the data sent us on some wild goose chases. For example, `TimeSeries.plot()`

and matplotlib's `plt.plot()`

would crop data differently, even when given the same time-interval.
What is more, they would crop slightly differently for different events. Eventually we had to
inspect each one by hand and shift accordingly.

Because some of the interesting signals were longer in time than others, not all of our sections of interest had the same number of points. To allow us to plot all of the events on the same X-axis, some points needed to be cut from the longest sections of interest. We used a method of **subsampling** to achieve this.
The lowest number of data points in any one of the events was 1400. Graphs with higher numbers of points (up to ~5000) were downsampled.

We observed that some strains were several orders of magnitude larger. For example, this wave from the Virgo observatory:

It just dwarfs all others, so we scaled those events to make them comparable to the ones from Hanford and Livingstson.

(You may notice that in the above image of the un-normalised Virgo events, you can see that the lines of data get shorter the further up they are positioned. This was intentional; it gave the image a perspective effect we found appealing. Unfortunately, as neat as this visual effect was, we ultimately ommited it. We wanted to make our final plot as recognisable as possible as a recreation of the original, which plotted its lines with equal lengths.)

Lastly, we applied a shape-defining function to each event to better reflect the expectations of numerical relativity. For example, no waves are expected after the black-hole merger, but a lot of noise remains in the data. To reflect this we attenuated the signal. This marked a departure from data processing as we entered into a more artistic endeavour.

We defined separate functions for the rising and falling phases of the plot. The rising phase was attenuated with a function proportional to $(t+k)^4$, and the post-merger with an $e^{-a(t-t_0)}$ function, $t_0$ being the approximate time of the merger. More specifically:

In [ ]:

```
ramp_up = ramp_up_height + 1/ramp_up_factor*(X+2.1)**4
ramp_up = np.where(ramp_up > 1, 1, ramp_up)
ramp_down = np.exp(-ramp_down_speed*(X-ramp_down_center/nr_points+1/2)) # 1/(100*X)
ramp_down = np.where(ramp_down > 1, 1, ramp_down)
```

The exact parameters of the function were changed depending on the individual
shapes of the original plots. The images below show
how the functions above attenuate the final versions of event **GW150914 - H1**.

Finally, we arrived at the *Wave Division*...

On the right side of each signal are the technical names of each event and the observatories from which they came.

We would like to thank Connectech Bermuda for kindly letting us use their facilities during the pandemic.

This plot and article is an example of what learning Physics in High School can be like. If you would like to support initiatives like these, order a Wave Divisions T-shirt ðŸ‘•! $\Longleftarrow$

The full code used to create the images is available on Github.

For those interested in the Physics and Math,
a companion article is available at alvarofeito.com/articles/ligo/.

A full treatment would require laying a deep foundation which involves calculus, differential equations, algebra, tensors, modern mechanics and Einstein's relativity. Instead we will try to use math and concepts available to a curious and motivated high school student.

- General Relativity for Laypeople - a primer
- The Equivalence principle & the deflection of light
- time dilation in GR
- What is space - nautil.us

- A no nonsense introduction to General Relativity - Sean Caroll
- Lecture notes on General Relativity - Sean Caroll
- Introduction to general relativity - H. Blöte
- Metric, curved space & General relavivity - Syksy RÃ¤sÃ¤nen
- Edbert's GR 1

- Studying the universe with gravitational waves
- Relativistic Astrophysics lecture - Wheeler (intuition + partial derivation)
- Piled Higher and Deeper Explanation

- https://www.quora.com/What-is-the-wave-equation-of-gravitational-waves
- https://www.pnas.org/content/pnas/113/42/11662.full.pdf

- https://arxiv.org/pdf/1005.4735.pdf
- https://www.tat.physik.uni-tuebingen.de/~kokkotas/Teaching/NS.BH.GW_files/GW_Physics.pdf
- https://www.ams.org/publications/journals/notices/201707/rnoti-p684.pdf