Gravitational Wave Divisions

Miles Oleksak & Alvaro Feito Boirac

Unknown Mergers

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!

Disclaimer: We find deplorable the historical idea of joy divisions and what the album title Unknown Pleasures may suggest.

What does the Graph show?

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.


The Laser Interferometer Gravitational-Wave Observatory (LIGO)

Hanford LIGO Interferometer

Image taken from LIGO's website.

Detecting Gravitational Waves

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.

Simplified operation of a gravitational wave observatory

Diagram taken from LIGO's website.

In Figure 1 of the diagram, a beamsplitter (green line) splits coherent light emitted from source (white box) into two beams, which then reflect off the mirrors (cyan rectangles) at the end of each arm. The reflected beams recombine (after approximately 280 round trips down the 4 km length). This "recombination" results in constructive or destructive interference depending on the path difference. An interference pattern is detected at the circle. In Figure 2, a gravitational wave (yellow) passes over the left arm and so changes its length slightly, resulting in a measurably different interference pattern detected at the circle.

To summarize, LIGO aims to detect gravitational waves by comparing the time of propagation of light in mutually orthogonal paths in the distorted space between freely suspended test masses separated by 4 km using the technique of laser interferometry.

Y-axis squiggles

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

Graphs lifted from Abbott et al.

Note that the y-axis on these plots is called strain. Strain is amount of deformation experienced by a body in the direction of force applied, divided by initial dimensions of the body. Or:

$$ \mathrm{Strain} = \frac{\Delta L }{L} $$

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.

Hanford, Livingston, and Virgo

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.

LIGO gravitational wave plots

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 Graph

The tools

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 data

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.

The process

Plotting unfiltered data

The gwpy library has an object called TimeSeries which allow us to read the data files (.hdf5 files) and contains methods (functions) to read, process & do simple plots of the data.

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

Raw data from Ligo

Filtering data

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:

A wave before applying filters

A wave after applying filters

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?

Step 1: Defining functions for whitening, bandpassing, and plotting

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.

Step 2: Selecting time and frequency cut-offs for each event

LIGO frequency-time plot

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.

Step 3: Correcting shift

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.

Preparing data to plot

Step 4: Subsampling

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.

Step 5: Normalising Virgo events

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

LIGO with spike from Virgo

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

Step 6: Shape-defining function

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. A wave before applying the function

A wave after applying the function

The final plot

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.

Acknowledgments:

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$

Diving Deeper

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.