Michelson response to GW signal and noise

In this notebook we want to model transfer functions and noises in an interferometer so that we can generate our first sensitivity plot.

After this notebook you will be able to:

  • compute the Michelson interferometer response to gravitational wave signal
  • compute the frequency-noise limited sensitivity of a Michelson interferometer

  • compute the quantum-noise limited sensitivity of a Michelson interferometer

  • read and explain a sensitivity curve

IMPORTANT make sure you run this next cell first to import all the code we need to run the notebook.

In [1]:
%run pykat_notebook_defaults.py
pykat.init_pykat_plotting()

import noise

                                              ..-
    PyKat 1.0.19          _                  '(
                          \`.|\.__...-""""-_." )
       ..+-----.._        /  ' `            .-'
   . '            `:      7/* _/._\    \   (
  (        '::;;+;;:      `-"' =" /,`"" `) /
  L.        \`:::a:f            c_/     n_'
  ..`--...___`.  .    ,  
   `^-....____:   +.      www.gwoptics.org/pykat

Imported matplotlib.pyplot as plt
Imported numpy as np

You can now use 'print(kat)' to display the Finesse code of a 'kat' object

How to model a transfer function - Example: power oscillations $\to$ photodiode output

To model a transfer function in Finesse there is a specfic recipe you must follow. This recipe is described in this section.

The steps you need to take when computing a transfer function are:

  • Add an fsig command. e.g. fsig sig1 laser amp 1 0 1, this adds a signal to the simulation called sig1, it takes the laser and modulates the amp (amplitude). It does this at 1Hz, with 0 phase offset, with an amplitude of 1.
  • Add a detector to demodulate this signal, e.g. pd1 tf $fs n1. This measures the signal we added with the fsig command before. $fs is a special variable which is always the same frequency as the fsig command.
  • Scan the signal frequency we are applying, e.g. xaxis sig1 f log 10k 1G 1000. To generate a transfer function we have to calculate how the modulation responds at a range of frequencies. This xaxis command changes the frequency which we modulate the amplitude.
  • Use yaxis abs:deg to get the complex valued transfer function, so we can see amplitude and phase.

In this example we compute how laser power noise affects a photodiode output. We create the power oscillations by modulating the laser's amplitude using the special command fsig. This will modulate the amplitude at particular frequency. We then measure how this modulation propagates to a photodiode output by demodulating a photodiode output at the same frequency as the oscillation. This process calculates the transfer function.

In [2]:
kat = finesse.kat()          # Initialising Finesse
kat.verbose = False          # Tells Finesse to talk less

code = '''
l laser 2 0 n0            # Laser with 0 offset wavelength (thus default 1064nm used)
s s1 1 n0 n1              # Space from laser to mirror (1 m)
fsig sig1 laser amp 1 0 1 # Laser power fluctuations (freq 1 Hz, phase 0, amplitude 1)

# Demodulate our signal to compute a transfer function
# The $fs variable is a special one in Finesse, its value is always the current signal frequency
pd1 tf $fs n1 

# Sweep the frequency of the signal applied from 10kHz to 1GHz in 1000 steps
xaxis sig1 f log 10k 1G 1000
yaxis abs:deg
'''
kat.parseCommands(code)                                  # passing the input text to the Finesse object
out = kat.run()

fig = out.plot(title='Bode diagram',
               ylabel='$|H(\omega)|$ $[\mathrm{W}/\sqrt{\mathrm{W}}]$',
               detectors=["tf"])

This plot shows the magnitude and phase of the transfer function. The units of the transfer function are always Watts per X, where X is the units of thing we are modulating. In this example we are modulating the amplitude, so X = $\sqrt{W}$.

Multiple fsigs

In Finesse you can use multiple fsig commands at once in a simulation. This is useful when you have two effects happening at once. For example, the differential degree of freedom means two mirrors are being moved at the same time but out of phase with each other, as one moves in one direction, the other mirror moves in the opposite direction. This signal can be described with the following two fsig commands:

fsig sig1 MX z 1 0 1
fsig sig2 MY z 1 180 1

Here we shake the MX and MY mirrors but out of phase by 180 degrees.

Options for fsig

This should be used as a reference for the tasks. It tells you what parameters can be modulated with an fsig command.

**Component** **Target** **Unit** **Description**
Mirror/Beamsplitter `phase` `W/deg` Mirror tuning modulation
Mirror/Beamsplitter `z` `W/m` Mirror tuning modulation
Laser `freq` `W/Hz` Modulates laser frequency
Laser `phase` `W/deg` Modulates laser phase
Laser `amp` `W/sqrt{W}` Modulates laser amplitude
Space `phase` `W/h` Modulates the optical phase for a
gravitational wave strain of `h`

Michelson model

Here we provide the Finesse code for a Michelson, use this for all the tasks below. It is based on this image:

You should look at the commands that have been used and understand how the components are connected together.

#---------------------------------------------------------------------------
# The core optical layout, without any detectors or simulation instructions.
#---------------------------------------------------------------------------

base = finesse.kat()          # Initialising Finesse
base.verbose = False          # Tells Finesse to talk less

base.parseCommands( """
## Parameters ##
const Pin 1M          # Laser power [W] (1M = 1 MW)
const LX 4000         # Length of X arm [m]
const LY 4000.01      # Length of Y arm [m]

## Laser & Beam splitter ##
l laser $Pin 0 n0                          # Laser (P=Pin, 0 wavelength offset from 1064nm)
s s1 1 n0 nBSb                             # Space from laser to beam splitter (1 m)
bs BS 0.5 0.5 0 45 nBSb nBSy nBSx nBSd     # Central beam splitter (0.5 refl, 0.5 trans, 0 deg. tuning,
                                           #                        45 degrees angle of incidence)

## X arm ##
s LX $LX nBSx nMX1          # Space between beam slitter and mirror MX (length = LX)
m MX 1 0 -45.003 nMX1 nMX2        # Mirror MX (1 refl, 0 trans, 0 deg. tuning)

## Y arm ##
s LY $LY nBSy nMY1          # Space between beam splitter and mirror MY
m MY 1 0 45.003 nMY1 nMY2       # Mirror MY (1 refl, 0 trans, 90 deg. tuning)

## Output ##
s sout 1 nBSd nout          # Space from BS to the photodiode (1 m)
""")

Sidenote on Michelson model (not needed to completing the task)

LIGO's laser can output 200W which is one of the most powerful lasers made. This Michelson model uses a laser that has an output power of 1MW! This is very high and also slightly unrealisitic. You will see in the lectures on Thursday how LIGO uses optical cavities to make a 200W laser into nearly megawatts of power. You should note that the end mirrors, MX and MY, have a differential tuning of $\pm 45.003$ degrees. This sets the opearting point of our Michelson to be at the dark fringe, so very little light is leaving the output port. The output port in this model is at the node nout. This is where we will add detectors to measure the optical fields to readout gravitational wave signals and noise. You should also see that the X and Y arm lengths are not exactly the same. This is a design decision that is too complicated to describe now. For all your models, the Y-arm should be 1cm longer than the X-arm.

Task 1: The gravitational wave signal

We now get to the stage we can model how a gravitational wave signal would appear in our detector. The gravitational wave will modulate the phase of the light travelling along a space. We can apply this type of signal to a space component in Finesse using: fsig sig1 space phase 1 0 1

(a) Use the code for the Michelson model given above:

  • Add two fsig sig1 space_name phase 1 0 1 commands to modulate the X-arm and Y-arm space at the same time, one of the fsigs phase should be 180 degrees different to the other, as the gravitational wave will make one arm bigger whilst the other gets smaller
  • Add in a photodiode to measure the transfer function at the Michelsons output port
  • Scan the signal from 1Hz to 100kHz in 1000 steps logarithmically
  • Plot the magnitude and phase of this transfer function

The units of this transfer function will be $W/h$ where $h$ is the gravitational wave strain. The strain value will typically be something very small, $h \sim 10^{-23}$.

(b) Is the Michelson more susceptible to low or high frequency gravitational waves?

(c) The magnitude of the transfer function goes to 0 at some frequencies, why is this?

In [3]:
base = finesse.kat()          # Initialising Finesse
base.verbose = False          # Tells Finesse to talk less

base.parseCommands( """
## Parameters ##
const Pin 1M          # Laser power [W] (1M = 1 MW)
const LX 4000         # Length of X arm [m]
const LY 4000.1       # Length of Y arm [m]

## Laser & Beam splitter ##
l laser $Pin 0 n0                          # Laser (P=Pin, 0 wavelength offset from 1064nm)
s s1 1 n0 nBSb                             # Space from laser to beam splitter (1 m)
bs BS 0.5 0.5 0 45 nBSb nBSy nBSx nBSd     # Central beam splitter (0.5 refl, 0.5 trans, 0 deg. tuning,
                                           #                        45 degrees angle of incidence)

## X arm ##
s LX $LX nBSx nMX1          # Space between beam slitter and mirror MX (length = LX)
m MX 1 0 -45.003 nMX1 nMX2        # Mirror MX (1 refl, 0 trans, 0 deg. tuning)

## Y arm ##
s LY $LY nBSy nMY1          # Space between beam splitter and mirror MY
m MY 1 0 45.003 nMY1 nMY2       # Mirror MY (1 refl, 0 trans, 90 deg. tuning)

## Output ##
s sout 1 nBSd nout          # Space from BS to the photodiode (1 m)
""")
In [4]:
kat = base.deepcopy()
kat.parseCommands("""
fsig sig1 LX 1 0 1
fsig sig1 LY 1 180 1

pd1 output $fs nout
xaxis sig1 f log 1 100k 1000
yaxis abs:deg
""")

out = kat.run()
fig = out.plot(ylabel="W/h")

Task 2: Frequency noise response

In this task we will compute how laser frequency noise couples to the Michelson output.

(a) Using the basic Michelson model from Task 1:

  • Add one fsig command to modulate both the frequency of the input laser. Use the reference table above to get the correct command
  • Add in the photodiode to measure a transfer function to the Michelson output port
  • Scan the signal from 1Hz to 100kHz in 100 steps logarithmically
  • Plot the magnitude and phase of this transfer function
In [5]:
kat = base.deepcopy()
kat.parseCommands("""
fsig sig1 laser freq 1 0 1

pd1 output $fs nout
xaxis sig1 f log 1 100k 1000
yaxis abs:deg
""")

out = kat.run()
fig = out.plot()

Task 3: Frequency noise limited sensitivity

In this task we will produce our first sensitivity plot. This will show how the laser frequency noise transfers to the output port where we measure for a gravitational wave, and then limits our sensitivity to measure GWs.

To compute the frequency noise limited sensitivity you must calculate:

$S_f(f) = \frac{N_f(f) |H_n(f)|}{|H_s(f)|}$

$N_f$ is the amplitude spectral density of the frequency noise of the laser, $H_n$ the transfer function from laser frequency to the output port, and $H_s$ being the transfer function of the signal to the output port. The units of this work out as:

$\frac{[Hz/\sqrt{Hz}] [W/Hz]}{[W/h]} = \frac{h}{\sqrt{Hz}}$

Which is the equivalent strain amplitude spectral density from frequency noise.

Here we generate some frequency noise that would come out of a laser that we have just bought. We will use the following code to generate this:

f = out.x # you should use the same frequency values that Finesse
          # outputs on the xaxis when you scan the frequency.
N_f_asd = noise.frequency_noise_ASD(f)

plt.loglog(f, N_f_asd)
plt.xlabel("Frequency [Hz]")
plt.ylabel("ASD [$Hz/\sqrt{Hz}$]")

image.png

(a) Use your results from tasks 1 and 2 for the transfer functions, calculate a fake $N_f$ as shown above, and then calculate $S_f(f)$ from above. Plot the $S_f(f)$ from 1Hz to 100kHz

In [6]:
import noise

kat = base.deepcopy()
kat.parseCommands("""
fsig sig1 LX 1 0 1
fsig sig1 LY 1 180 1

pd1 output $fs nout
xaxis sig1 f log 1 100k 1000
yaxis abs:deg
""")

out = kat.run()
Hs = abs(out['output'])

kat = base.deepcopy()
kat.parseCommands("""
fsig sig1 laser freq 1 0 1

pd1 output $fs nout
xaxis sig1 f log 1 100k 1000
yaxis abs:deg
""")

out = kat.run()
Hn = abs(out['output'])

f = out.x               

N_f_asd = noise.frequency_noise_ASD(f)

# Here we compute the frequency noise limited sensitivity
S_f_asd = N_f_asd * Hn / Hs            

plt.loglog(f, S_f_asd)
plt.xlabel("Frequency [Hz]")
plt.ylabel("Sensitivity [$h/\sqrt{Hz}$]")
Out[6]:
<matplotlib.text.Text at 0x11a7fe518>

Task 3: Mirror motion to output port transfer function - seismic noise

In this tasks we will model how the motion of a mirror appears in the output detector. This is an important transfer function, as many sources of noise will shake the mirrors individually. For example, seismic noise from the Earth will shake the mirrors in the detectors. However, this shaking is not coherent between all the mirrors because the mirrors are so far away from each other.

$S_s(f) = \frac{N_s(f) |H_m(f)|}{|H_s(f)|}$

where $S_s(f)$ is the equivalent strain amplitude spectral density from seismic noise, $N_s(f)$ is the motion in meters of the mirror due to seismic noise, $H_m(f)$ is the transfer function from the mirror to the output detector, and $H_s(f)$ is the transfer function of the gravitational wave signal.

You can get the seismic noise using this function:

f = out.x # you should use the same frequency values that Finesse
          # outputs on the xaxis when you scan the frequency.
N_s_asd = noise.seismic_noise_ASD(f)

(a) Use the table above to find the fsig command to modulate the mirror position in units of [W/m].

(b) The seimic noise will shake the two end mirrors separately. Compute and plot the transfer function in units of [W/m] for the:

  • MX mirror to the output port
  • MY mirror to the output port

(c) Compute and plot the equivalent strain amplitude spectral density from seismic noise. For this you will need to add the seismic noise from the MX and MY mirror.

In [7]:
kat = base.deepcopy()
kat.parseCommands("""
fsig sig1 MX z 1 0 1

pd1 output $fs nout
xaxis sig1 f log 1 100k 1000
yaxis abs:deg
""")

out = kat.run()
Hmx = abs(out['output']) # Here we save the transfer function to use later too
fig = out.plot(ylabel="MX [W/m]")

kat = base.deepcopy()
kat.parseCommands("""
fsig sig1 MY z 1 0 1

pd1 output $fs nout
xaxis sig1 f log 1 100k 1000
yaxis abs:deg
""")

out = kat.run()
Hmy = abs(out['output']) # Here we save the transfer function to use later too
fig = out.plot(ylabel="MY [W/m]")
In [8]:
f = out.x               

N_s_asd = noise.seismic_noise_ASD(f)

# Here we compute the frequency noise limited sensitivity
S_f_asd = np.sqrt((N_s_asd * Hmx / Hs)**2 + (N_s_asd * Hmy / Hs)**2)

plt.loglog(f, S_f_asd)
plt.xlabel("Frequency [Hz]")
plt.ylabel("Sensitivity [$h/\sqrt{Hz}$]")
Out[8]:
<matplotlib.text.Text at 0x11add9080>

Task 4: The sensitivity plot

Now we are ready to make a sensitivity plot for our Michelson! We will add the seismic noise, frequency noise, and finally the quantum noise and compute the total sensitivity. Remember that these are not all the noise sources! This is just generating an example plot for you to experiment with.

You will hear more about quantum noise in the lectures tomorrow, do not worry about understanding it! Just run the code below to compute the quantum noise limited sensitivity that you can use in your plots. The qnoised detector in this code just measures quantum noise amplitude spectral density at the output port.

kat = base.deepcopy()
kat.parseCommands("""
fsig noise 1
qnoised qn 1 $fs nout
xaxis noise f log 1 100k 1000
""")

out = kat.run()
f = out.x

# Here we compute the Quantum noise limited sensitivity
# The variable Hs is the signal transfer function for a gravitational wave
# NOTE! Hs is not computed here, you must replace Hs with the transfer function computed in Task 1!
S_q_asd = out['qn'] / Hs

plt.loglog(f, S_q_asd, label="Quantum")
plt.xlabel("Frequency [Hz]")
plt.ylabel("Sensitivity [$h/\sqrt{Hz}$]")
plt.legend()
plt.tight_layout()

NOTE the variable Hs you should replace with the GW signal transfer function from task 1! You should then be able to make a plot like this:

image.png

(a) Your task is to plot on the same graph:

  • The quantum noise limited sensitivity
  • The frequency noise limited sensitivity
  • The seismic noise limited sensitivity
  • The sum of these three noises to get the total sensitivity

Tip 1: Remember, when adding noise sources $A$, $B$, and $C$ together, $\sqrt{A^2 + B^2 + C^2}$.

Tip 2: When plotting this curve it might be useful to limit the y-axis using the function plt.ylim(lower, upper). For example, plt.ylim(1e-24, 1e-20) to limit the sensitivity between 1e-24 and 1e-20.

(b) How much do we need to reduce the frequency noise by to measure a gravitational wave at 100Hz with a strength of $h\sim10^{-23}$?

In [65]:
kat = base.deepcopy()
kat.parseCommands("""
fsig noise 1
qnoised qn 1 $fs nout
xaxis noise f log 1 100k 1000
""")

out = kat.run()
f = out.x   # Here is the frequency vector we will use for the sensitivity

# The quantum noise limited strain sensitivity
S_q_asd = out['qn'] / Hs 



# The frequency noise limited strain sensitivity
freq_factor = 200000 # The frequency factor we need to improve our frequency noise by
S_f_asd_1 = S_f_asd/freq_factor 



# The seismic noise limited strain sensitivity
N_s_asd = seismic_noise_ASD(f)
# Here we compute the frequency noise limited sensitivity
S_s_asd = np.sqrt((N_s_asd * Hmx / Hs)**2 + (N_s_asd * Hmy / Hs)**2)


plt.loglog(f, S_q_asd, label="Quantum")
plt.loglog(f, S_f_asd_1, label="Frequency")
plt.loglog(f, S_s_asd, label="Seismic")

total = np.sqrt(S_q_asd**2 + S_f_asd_1**2 + S_s_asd**2)
plt.loglog(f, total, label="Total")

plt.ylim(1e-24, 1e-20)
plt.legend(loc=0)
plt.xlabel("Frequency [Hz]")
plt.ylabel("Sensitivity [$h/\sqrt{Hz}$]")
plt.margins(0, 0)
In [ ]: