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 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.
%run pykat_notebook_defaults.py
pykat.init_pykat_plotting()
import noise
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:
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 1
Hz, with 0
phase offset, with an amplitude of 1
.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.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.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.
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}$.
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.
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` |
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)
""")
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.
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:
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 fsig
s phase should be 180 degrees different to the other, as the gravitational wave will make one arm bigger whilst the other gets smallerThe 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?
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)
""")
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")
In this task we will compute how laser frequency noise couples to the Michelson output.
(a) Using the basic Michelson model from Task 1:
fsig
command to modulate both the frequency of the input laser. Use the reference table above to get the correct commandkat = 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()
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}$]")
(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
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}$]")
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 portMY
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.
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]")
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}$]")
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:
(a) Your task is to plot on the same graph:
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}$?
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)