40m Quantum Limited Sensitivity - Simple Michelson

In this example, we check that the most simple version of the model we could possibly make, matches our analytical expressions and expectations.

We will build the 40m model, with no substrates, a simple beamsplitter, no auxillary modulations and simple suspensions.

We start with our required imports:

import finesse
import numpy as np
from scipy.constants import speed_of_light, pi, hbar
from finesse_40m.factory import FortyMeterFactory
import finesse.analysis.actions as fac
from finesse.analysis.actions.axes import Noxaxis as Noxaxis
from finesse.analysis.actions.axes import Xaxis as Xaxis
import finesse.components as fc
from finesse.plotting import bode
import matplotlib.pyplot as plt
finesse.init_plotting()
import pprint

We then disable all the “realism”

factory = FortyMeterFactory()
factory.reset()
factory.options.add_11MHz=False
factory.options.add_33MHz=False
factory.options.add_55MHz=False
factory.options.add_165MHz=False
factory.options.BS_type='fake'
factory.options.LSC.add_DOFs = False
factory.options.add_detectors = "power"

# Print out build options
pprint.PrettyPrinter(indent=4, sort_dicts=True).pprint(
    factory.options.toDict()
)

# Build model
model = factory.make()
model.modes("off")  # planewave

{   'ASC': {   'add': False,
               'add_AC_loops': True,
               'add_DOFs': True,
               'add_locks': True,
               'add_output_detectors': True,
               'add_readouts': True,
               'close_AC_loops': False},
    'BS_trans_arm': 'X',
    'BS_type': 'fake',
    'ETMAR': False,
    'INPUT': {'add_IMC_and_IM1': False, 'set_IMC_mode': True},
    'LSC': {   'add_AC_loops': False,
               'add_DOFs': False,
               'add_locks': False,
               'add_output_detectors': False,
               'add_readouts': True,
               'close_AC_loops': False},
    'OUTPUT': {'readout': None},
    'PR2_PR3_substrates': True,
    'QUAD_suspension_kwargs': {},
    'QUAD_suspension_model': None,
    'add_118MHz': False,
    'add_11MHz': False,
    'add_165MHz': False,
    'add_33MHz': False,
    'add_45MHz': False,
    'add_55MHz': False,
    'add_9MHz': False,
    'add_detectors': 'power',
    'add_transmon': False,
    'fake_prc_gouy': False,
    'fake_src_gouy': False,
    'materials': {   'test_mass_substrate': <finesse.materials.Material object at 0x7f34d5743b70>},
    'test_mass_wedge': False}

To see interesting radiation pressure, we are going to need a lot of power. 150 kW arm power will be enough for interesting effects.

We then disable all of the extra mirrors and check that we are at the operating point.

dc_offset = 125e-6 # degrees
model.L0.P=300e3
model.PRM.set_RTL(R=0,T=1)
model.PRMAR.set_RTL(R=0,T=1)
model.ITMX.set_RTL(R=0,T=1)
model.ITMY.set_RTL(R=0,T=1)
model.ITMXAR.set_RTL(R=0,T=1)
model.ITMYAR.set_RTL(R=0,T=1)
model.SRM.set_RTL(R=0,T=1)
model.SRMAR.set_RTL(R=0,T=1)
model.ETMY.phi=90 - dc_offset

# Output powers
sim = model.run(fac.Noxaxis())
for detector in sim.detectors:
    if detector[0] == 'P':
        print(f"{detector} = {sim[detector]:.1f} W")

Px = 149835.6 W
Pinx = 149841.4 W
Py = 149820.7 W
Piny = 149826.5 W
Pprc = 300000.0 W
Psrc = 0.0 W
Pin = 300000.0 W
PreflPRM = 299325.2 W
Prefl = 299325.2 W
Ppop = 291.2 W

We can now parse some additional katscript to suspend our mirrors and set up our simulation.

model2 = model.deepcopy()
model2.parse("""
# Differentially modulate the arm lengths
# These are toy numbers only, for a simple example
pendulum susX ETMX.mech mass=1 fz=0.8 Qz=100
pendulum susY ETMY.mech mass=1 fz=1.2 Qz=100

# Differentially modulate the arm lengths
fsig(1)
sgen darmx LX.h
sgen darmy LY.h phase=180

# Output the full quantum noise limited sensitivity
qnoised NSR_with_RP SRM.p2.o nsr=True
# Output just the shot noise limited sensitivity
qshot NSR_without_RP SRM.p2.o nsr=True

# We could also display the quantum noise and the signal
# separately by uncommenting these two lines.
# qnoised noise srm.p2.o
# pd1 signal srm.p2.o f=fsig
""")

Now, we can encode our analytical expression for the Noise-to-Signal-Ratio (NSR) of the Michelson response to a GW fluctuation.

def michelson_response(omega_gw, model, out):
    """ Equation 6.21 in Bond2017

    Bond2017 = Interferometer Techniques for Gravitational Wave Detection
    Bond, Brown, Freise & Strain
    """
    P0 = out['Pprc']
    omega0 = 2*pi*speed_of_light / model.lambda0
    common_arm = 0.5*(model.LX.L + model.LY.L ) # + model.lx1.L + model.ly1.L
    delta_L = (model.ETMY.phi - model.ETMX.phi)

    term1 = np.sqrt(2*hbar / (P0*omega0))
    term2 = omega_gw / np.sin(omega_gw * common_arm / speed_of_light)
    return np.abs(term1*term2)

Between 100 Hz and 1 MHz we expect a flat response. We don’t need to plot that, so we can optimize our simulation by splitting it into two runs and using named axes.

hf_ax = Xaxis(model2.darmx.f, 'log', 400e3, 20e6, 1000, name= 'HF')
lf_ax = Xaxis(model2.darmx.f, 'log', 0.1, 200, 1000, name='LF')

out = model2.run(fac.Series(lf_ax,hf_ax))

Now strain is defined as \(\Delta L / L\). In our simulation we apply the strain twice, once on the X-ARM and once on the Y-ARM, with 180 degrees phase. This means the total signal is a factor 2 larger than we expect. As a result, our Noise-to-Signal-Ratio is a factor 2 smaller than we expect. Therefore, we multiply our simulated experimental data by 2.

Plotting our solution, we find that at high frequencies we have agreement between the shot noise limited sensitivity of the Finesse simulation and the analytics.

At low frequencies, we can see that Finesse is modelling an additional effect - radiation pressure!

fig, (ax1,ax2) = plt.subplots(ncols=2,sharey=True,figsize=(7,3))
ax1.loglog(out['LF'].x0,2*out['LF']['NSR_without_RP'], c='b', ls='-',
    label='NSR without RP'
)
ax2.loglog(out['HF'].x0,2*out['HF']['NSR_without_RP'], c='b', ls='-')

ax1.loglog(out['LF'].x0,2*out['LF']['NSR_with_RP'], c='r', ls='--',
    label='NSR with RP'
)
ax2.loglog(out['HF'].x0,2*out['HF']['NSR_with_RP'], c='r', ls='--')

analytic = {}
analytic['HF'] = michelson_response(2*pi*out['HF'].x0, model2, out['HF'])
analytic['LF'] = michelson_response(2*pi*out['LF'].x0, model2, out['LF'])

ax1.loglog(out['LF'].x0,analytic['LF'], c='g', ls=':', label='Analytical')
ax2.loglog(out['HF'].x0,analytic['HF'], c='g', ls=':')

# hide the spines between ax and ax2
ax1.spines['right'].set_visible(False)
ax2.spines['left'].set_visible(False)
ax1.yaxis.tick_left()
#ax2.yaxis.set_ticklabels([])
ax1.tick_params(labelright='off')
ax2.tick_params(axis='y',which='both',width=0)

ax1.legend()

ax1.set_ylabel('Sensitivity [1/sqrt(Hz)]')

Text(0, 0.5, 'Sensitivity [1/sqrt(Hz)]')

Last of all, we add a check to ensure that the documentation build will fail if this validation test fails.

assert np.all(
    np.isclose(2*out['LF']['NSR_without_RP'], analytic['LF'], rtol=1e-2, atol=1e-24)
)
assert np.all(
    np.isclose(2*out['HF']['NSR_without_RP'], analytic['HF'], rtol=1e-2, atol=1e-24)
)
assert np.all(
    np.isclose(
        2*out['HF']['NSR_without_RP'], 2*out['HF']['NSR_with_RP'],
        rtol=1e-2, atol=1e-24)
)