Supplementary code for:
Korobko, M., Ma, Y., Chen, Y. and Schnabel, R.
Quantum expander for gravitational - wave observatories.
Light Sci Appl 8, 118 (2019) 
DOI: 10.1038/s41377-019-0230-2

Mikhail Korobko
Email: mkorobko@physnet.uni-hamburg.de
Affiliation: Institut für Laserphysik und Zentrum für Optische Quantentechnologien, Universität Hamburg, 
Luruper Chaussee 149, 22761 Hamburg, Germany

Abstract

This notebook implements the equations in the Supplementary Materials for the paper and shows how to reproduce Fig.1 in the main text of the paper.  It can further be used to reproduce other plots in the paper. 

The parameters are chosen to be not specific to any existing or planned GW observatory, but the code presented here is general and can be used to simulate those as well. It can further be used for  simulating other internal squeezing proposals, such as [39] and [42].

[39] Korobko, M. et al. Beating the standard sensitivity-bandwidth limit of cavity-enhanced interferometers with internal squeezed-light generation. Phys. Rev. Lett. 118, 143601 (2017).
[42] Korobko, M., Khalili, F. Y. & Schnabel, R. Engineering the optical spring via intra-cavity optical-parametric amplification. Phys. Lett. A 382, 2238–2244 (2018).

Parameters

Main;
L= 20000; (*Arm cavity length*)
ls =   56;(*SE cavity length*)
Ti= (Sqrt[ 0.014 5   ]); (*Input mirror amplitude transmission*)
Ts=  (Sqrt[0.35]); (*SE mirror amplitude transmission*)
m= 200 ; (*Mirror mass*)
Ic=  2 4000 (10^3); (*Intra-cavity power, factor of 2 for scaling law*) 

Losses;
Note that for simplicity the internal loss is assumed zero instead the detection efficiency represents total loss in the system. The complete treatment as in Fig.S1  can be achieved by setting other values ;
η = 0.995; (*Detection efficiency*)
λsr = 0 1500 (10^(-6));(*Intra SE cavity loss*)
Tend=0 (Sqrt[100 (10^(-6))]);(*End mirror amplitude transmission - intra arm cavity loss*)

Internal squeezing;
qth = -(1/2)(Log[Rs]) ;(*Approximate parametric threshold value*)
θ=π/2 (*Internal squeezing phase*);

External squeezing;
ExtSqzdB = 0 (*External squeezing, dB*);
qext = (1/2)Log[(10^(-ExtSqzdB/10))](*External squeezing, factor*);
ϕext = 0π/2(*External squeezing angle*);

Tunings;
ϕ1= 0 (*First phase of SE cavity*);
ϕ2= π/2(*Second phase of SE cavity*);
ζ=0(*Homodyne detection phase*);
γf = 0 (*Filter cavity for variational readout linewidth*);
δa=0(1/τ)(π/2) (*Detuning of the arm cavity*);

Other parameters;
ℏ=1.054571628535353535353*(10^(-34)) (*Plank constant*);
c =  299792458 (*Speed of light*);
λ= 1550 (10^(-9))(*Central wavelength*);
kp=2 (π/λ) (*Wavenumber*);
ωm= 0 (*Resonanse frequency of the suspension - assume free mass*);
γm= 0 (*Mechanical loss of the suspension*);
Ri=(Sqrt[1-(Ti^2)]) (*Input mirror amplitude reflectivity*);
Rs=(Sqrt[1-(Ts^2)])(*SE mirror amplitude transmission*);
Rend=(Sqrt[1-(Tend^2)])(*End mirror amplitude transmission*);
τ=L/c(*Arm cavity travel time*);
Ω = 2 π f(*Angular frequency*);
fend = 7400; (*Change the end plotting frequency*)

Matrices

𝕆[ϕ_]:=(({{Cos[ϕ], -Sin[ϕ]}, {Sin[ϕ], Cos[ϕ]}}))(*Eq. 45*);
𝕊 = (({{(E^q), 0}, {0, (E^(-q))}}))(*Eq. 47*);
𝕊ext = (({{Cosh[qext]+Cos[2ϕext]Sinh[qext], Sin[2ϕext]Sinh[qext]}, {Sin[2ϕext]Sinh[qext], Cosh[qext]-Cos[2ϕext]Sinh[qext]}}))(*Eq. 93*);
Ee = (Sqrt[2 (Ic/(ℏ kp c))])(({{1}, {0}}))(*Eq. 86*);
𝕄12 = 𝕆[ϕ1]. 𝕆[θ].𝕊. 𝕆[θ]\[HermitianConjugate]. 𝕆[ϕ2](E^(I Ω (ls/c)))//ComplexExpand//Simplify(*Eq. 58, ℳ[ϕ,ψ]*);
𝕄21 = 𝕆[ϕ2]. 𝕆[θ].𝕊. 𝕆[θ]\[HermitianConjugate]. 𝕆[ϕ1](E^(I Ω (ls/c)))//ComplexExpand//Simplify(*Eq. 58, ℳ[ψ,ψ]*);
𝕀 = (({{1, 0}, {0, 1}}));
𝕐 = 𝕆[π/2](*Eq. 46*);
𝔻b = Inverse[𝕀 + Ri Rs (1-λsr)𝕄21.𝕄12]//Simplify(*Eq. 59*);
𝔻d = Inverse[𝕀 + Ri Rs(1-λsr) 𝕄12.𝕄21]//Simplify(*Eq. 60*);
ℝb = Rs 𝕀 + (1-λsr)Ri (Ts^2)𝕄12.𝔻b.𝕄21//Simplify(*Eq. 63*);
𝕋b = (Sqrt[(1-λsr)])Ti Ts 𝕄12.𝔻b//Simplify(*Eq. 65*);
ℝd = Ri 𝕀+ (1-λsr)Rs (Ti^2)𝕄21.𝔻d.𝕄12//Simplify(*Eq. 64*);
𝕋d = (Sqrt[(1-λsr)])Ti Ts 𝕄21.𝔻d//Simplify(*Eq. 66*);
𝕃b1 = -Ts Ri (Sqrt[1-λsr])λsr 𝕄21.𝔻b.𝕄12(*Eq. 67*);
𝕃b2 = Ts Ri Rs (Sqrt[λsr])(1-λsr) 𝕄21.𝔻b.𝕄12 - (Sqrt[λsr])(*Eq. 68*);
𝕃d1 = -Ri Rs Ti (Sqrt[λsr])(1-λsr) 𝕄12.𝔻d.𝕄21 + (Sqrt[λsr])𝕄12(*Eq. 69*);
𝕃d2 = -(Sqrt[λsr (1-λsr)])Ti Rs 𝕄12.𝔻d(*Eq. 70*);
𝔻c = Inverse[𝕀 - Rend 𝕆[δa τ].𝕆[δa τ].ℝd (E^(2I Ω τ))]//Simplify(*Eq. 73*);
𝔻e = Inverse[𝕀 - Rend 𝕆[δa τ].ℝd.𝕆[δa τ] (E^(2I Ω τ))]//Simplify(*Eq. 74*);
ℝ= -ℝb +Rend 𝕋b.𝔻c.𝕆[δa τ].𝕆[δa τ].𝕋d (E^(2I Ω τ))(*Eq. 76*);
𝕋 = Tend 𝕋b.𝔻c.𝕆[δa τ](E^(I Ω τ))(*Eq. 77*);
ℤ= 2kp Rend 𝕋b.𝔻c.𝕆[δa τ].𝕐.Ee (E^(I Ω τ))(*Eq.78*);
𝕃1 = Rend 𝕋b.𝔻c.𝕆[δa τ].𝕆[δa τ].𝕃d1 (E^(2I Ω τ)) + 𝕃b1(*Eq.79*);
𝕃2 = Rend 𝕋b.𝔻c.𝕆[δa τ].𝕆[δa τ].𝕃d2 (E^(2I Ω τ))+ 𝕃b2(*Eq.80*);
𝕃v = (1+(Rend^2))𝔻e.𝕆[δa τ].ℝd.𝕆[δa τ](E^(I Ω τ)) + Rend(*Eq. 84*);

𝒦 =- 4 ℏ (kp^2) Ee\[Transpose].𝔻e.𝕆[δa τ].ℝd.𝕆[δa τ].𝕐.Ee (E^(2I Ω τ))(*Eq. 85*);
χeff = (1/(m (((ωm)^2) - (Ω^2)) + 2I γm Ω+𝒦))(*Eq. 88*);
ℍ = (({{Cos[ζ-ArcTan[2(((γf)^2)/(Ω^2))]]}, {Sin[ζ-ArcTan[2(((γf)^2)/(Ω^2))]]}}))(*Eq. 90, 99*);

Spectral densities

Sxx = ((ℍ\[Transpose].(ℝ.𝕊ext.𝕊ext\[HermitianConjugate].ℝ\[HermitianConjugate] + 𝕋.𝕋\[HermitianConjugate]+𝕃1.𝕃1\[HermitianConjugate]+𝕃2.𝕃2\[HermitianConjugate]).ℍ)/((Abs[ℍ\[Transpose].ℤ])^2))+((1-η)/η)(1/((Abs[ℍ\[Transpose].ℤ])^2))(*Eq. 95*);
SFF = (ℏ^2)(kp^2)(((1+(Rend^2)))^2)Ee\[Transpose].𝔻e.𝕆[δa τ].(𝕋d.𝕊ext.𝕊ext\[HermitianConjugate].𝕋d\[HermitianConjugate] + 𝕃d1.𝕃d1\[HermitianConjugate]+𝕃d2.𝕃d2\[HermitianConjugate]).𝕆[δa τ]\[HermitianConjugate].𝔻e\[HermitianConjugate].Ee + (ℏ^2)(kp^2) (Tend^2) Ee\[Transpose].𝕃v.𝕃v\[HermitianConjugate].Ee(*Eq. 96*);
SxF = ℏ (kp/(ℍ\[Transpose].ℤ))((1+(Rend^2))ℍ\[Transpose].(-ℝ.𝕊ext.𝕊ext\[HermitianConjugate].𝕋d\[HermitianConjugate] + 𝕃1.𝕃d1\[HermitianConjugate]+𝕃2.𝕃d2\[HermitianConjugate]).𝕆[δa τ]\[HermitianConjugate].𝔻e\[HermitianConjugate].Ee (E^(-I Ω τ)) + Tend ℍ\[Transpose].𝕋.𝕃v\[HermitianConjugate].Ee)(*Eq. 97*);

Sf = (Sxx/((Abs[χeff])^2)) + SFF+ 2 Re[((χeff)^(-1)) SxF];
Sx = Sxx + ((Abs[χeff])^2) SFF+2 Re[χeff\[Conjugate] SxF](*Eq. 94*);
Sh = Sx (4/((m^2)(L^2) (Ω^4)((Abs[χeff])^2)))(((((Sin[Ω τ])/(Ω τ))))^(-2))(*Eq. 98*);
SQLh = ((8ℏ Abs[m (Ω^2)])/( (L^2) ((Abs[m (Ω^2)])^2)));

Plot spectral densities

ShBase = (Sqrt[Sh]) (10^24)/.q-> 0;
ShQE = (Sqrt[Sh]) (10^24)/.q->qth;

LogLogPlot[{(Sqrt[SQLh]) (10^24),Re[ShBase],Re[ShQE]},{f,10,7000}, PlotRange->{{0.2, 5}},PlotLegends->{"SQL","Baseline observatory", "Quantum Expander, 0.5% loss"},Frame-> True,FrameLabel->{"Frequency, Hz", "Strain sensitivity, (10^(-24))/(Sqrt[Hz])"},PlotStyle->{Black, Darker[Blue],Darker[Red]}]