Session 01 — Solving the detector by matrix inversion
Module: topics/optics/matrix_methods/ · Session: 1 of 3
Spine: Block 1 — the detector is a hierarchy of linear systems
Author: GW101 curriculum · Purpose: Learn one recipe — draw the system → write the
node-coupling (adjacency) matrix M by inspection → solve (I − M)x = x_in — and apply it to a
Fabry–Perot cavity, a SISO feedback loop, and an audio-sideband transfer function.
Requirements: Python 3, numpy, matplotlib, and shared/code/conventions.py (the canonical
sign/phase conventions). Runs top-to-bottom; asserts give private pass/fail feedback.
Overview¶
Goal: see that a cavity and a feedback loop are the same linear-algebra problem, and that scaling up to MIMO / a full interferometer only enlarges
M.Inputs: mirror amplitude coefficients
(r, t), cavity lengthL, optical frequency.Outputs: circulating/reflected fields, buildup vs detuning, closed-loop transfer functions
S,T, and the cavity polef_c— plus four mandatory figures.Prerequisites: the
exp(−iωt)/exp(+ikL)conventions inshared/code/conventions.py. If unfamiliar with complex field amplitudes, start withoptics/ifo_topology.
The notebook fades: a fully worked cavity example → a guided feedback-loop step → an independent transfer-function task. Skip ahead to the open task if it is already easy.
Running in your browser (JupyterLite)? The first run downloads numpy/matplotlib as WebAssembly — give the imports cell ~30–60 s to finish before running the rest, top to bottom.
import pathlib, sys
import numpy as np
import matplotlib.pyplot as plt
# Canonical conventions (single source of truth). In a repo checkout this imports
# shared/code/conventions.py; in the browser (JupyterLite/Pyodide) there is no repo filesystem,
# so fall back to inline definitions of the few symbols this notebook uses.
root = pathlib.Path.cwd().resolve()
while not (root / "shared" / "code" / "conventions.py").exists() and root != root.parent:
root = root.parent
if (root / "shared" / "code" / "conventions.py").exists():
sys.path.insert(0, str(root / "shared" / "code"))
import conventions as cv
else:
import types
cv = types.ModuleType("conventions")
cv.C_SI = 299_792_458.0
cv.TIME_DEPENDENCE = "exp(-i ω t)"
cv.PROPAGATION_PHASE_SIGN = +1
cv.power_from_field = lambda E: np.abs(E) ** 2
np.random.seed(0) # determinism (no randomness is used, but fix it anyway)
plt.rcParams["figure.dpi"] = 110
print("conventions:", cv.TIME_DEPENDENCE, "| propagation exp(+ikL):", cv.PROPAGATION_PHASE_SIGN == 1)The recipe¶
A linear system is a graph: nodes are complex amplitudes, edges multiply by a complex gain. Each node = sum of incoming edges + any external injection:
M_ij is the gain of the edge from node j into node i — read straight off the drawing.
Solving is one np.linalg.solve; the inverse (I − M)^{-1} is the closed-loop response.
Part A — Fabry–Perot cavity (worked)¶
Nodes around the round trip, ITM (r1,t1), ETM (r2,t2), one-way phase φ = kL:
E1leaves ITM toward ETM,E2 = e^{iφ}E1arrives at ETM,E3 = -r2 E2leaves ETM back,E4 = e^{iφ}E3arrives at ITM,ITM closes the loop:
E1 = t1 E_in - r1 E4.
The minus signs are the convention: a field reflecting off the HR/cavity side of a mirror
picks up -r (the field reflecting off the substrate side gets +r). This is what keeps the
cavity’s reflected field physical, |r_cav| <= 1.
Predict¶
Before running: sketch circulating power vs detuning (where is the peak?) and the phase of the reflected field as you scan through resonance. Write your guess down, then run the cells.
# Realistic aLIGO-like arm cavity
L = 4000.0 # m
T1 = 0.014 # ITM power transmissivity
T2 = 5e-6 # ETM power transmissivity (near-perfect reflector)
r1, t1 = np.sqrt(1 - T1), np.sqrt(T1)
r2, t2 = np.sqrt(1 - T2), np.sqrt(T2)
lam = 1064e-9
k = 2 * np.pi / lam
def fp_coupling_matrix(phi, r1, r2):
'''Coupling (adjacency) matrix M for nodes (E1, E2, E3, E4), by inspection.'''
e = np.exp(1j * phi)
# Side-dependent reflection sign (Finesse 'real' convention): the returning field reflects
# off the HR/cavity side of each mirror, so it picks up -r (not +r).
return np.array([[0, 0, 0, -r1],
[e, 0, 0, 0 ],
[0, -r2, 0, 0 ],
[0, 0, e, 0 ]], dtype=complex)
def fp_solve(phi, Ein=1.0):
'''Solve (I - M) E = x_in for the four cavity field nodes; return (E, E_refl).'''
M = fp_coupling_matrix(phi, r1, r2)
x_in = np.array([t1 * Ein, 0, 0, 0], dtype=complex)
E = np.linalg.solve(np.eye(4, dtype=complex) - M, x_in)
E_refl = r1 * Ein + t1 * E[3] # reflection off ITM + leakage of returning field
return E, E_refl
# Sanity: the matrix solution must equal the analytic self-consistency result.
phi_test = 0.013
E, E_refl = fp_solve(phi_test)
E1_analytic = t1 / (1 - r1 * r2 * np.exp(2j * phi_test))
refl_analytic = r1 - t1**2 * r2 * np.exp(2j * phi_test) / (1 - r1 * r2 * np.exp(2j * phi_test))
assert np.isclose(E[0], E1_analytic, rtol=1e-12), "circulating field mismatch"
assert np.isclose(E_refl, refl_analytic, rtol=1e-12), "reflected field mismatch"
print("matrix solution matches analytic self-consistency to machine precision")Model + Measure — buildup and reflection vs detuning¶
Sweep the single-pass phase φ through a resonance. Map it to frequency via f = φ·FSR/π
(round-trip phase advances by 2π over one free spectral range).
FSR = cv.C_SI / (2 * L)
finesse = np.pi * np.sqrt(r1 * r2) / (1 - r1 * r2)
f_c = FSR / (2 * finesse)
print(f"FSR = {FSR:,.0f} Hz | finesse = {finesse:.1f} | cavity pole f_c = {f_c:.2f} Hz")
phi = np.linspace(-0.05, 0.05, 301) # kept modest so the in-browser (Pyodide) kernel keeps up
E_all = np.array([fp_solve(p)[0] for p in phi])
refl = np.array([fp_solve(p)[1] for p in phi])
P_circ = cv.power_from_field(E_all[:, 0])
P_circ_analytic = cv.power_from_field(t1 / (1 - r1 * r2 * np.exp(2j * phi)))
assert np.allclose(P_circ, P_circ_analytic, rtol=1e-10), "Airy curve mismatch"
# Energy conservation: reflected + transmitted-through-ETM must equal the input (lossless).
P_trans_etm = cv.power_from_field(t2 * E_all[:, 1]) # |t2 * E2|^2, E2 = field at the ETM
assert np.all(np.abs(refl) <= 1 + 1e-9), "reflectivity exceeds 1 (non-unitary mirror!)"
assert np.allclose(cv.power_from_field(refl) + P_trans_etm, 1.0, atol=1e-9), "energy not conserved"
print(f"energy conserved across the sweep; on-resonance |r_cav| = {np.abs(refl).min():.4f} "
f"(critical coupling -> 0 when r2 = r1)")
f_axis = phi * FSR / np.pi # Hz
fig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 4))
a1.plot(f_axis, P_circ, lw=2, label="matrix solve")
a1.plot(f_axis, P_circ_analytic, "--", lw=1, label="analytic")
a1.set(xlabel="detuning (Hz)", ylabel=r"circulating $|E_1|^2$ (a.u.)",
title="Fig 1 — Cavity buildup vs detuning")
a1.legend(); a1.grid(alpha=0.3)
a2.plot(f_axis, np.abs(refl), lw=2, color="C3")
a2b = a2.twinx()
a2b.plot(f_axis, np.unwrap(np.angle(refl)), lw=1.5, color="C0")
a2.set(xlabel="detuning (Hz)", ylabel="|reflectivity|", title="Fig 2 — Reflected field")
a2.tick_params(axis="y", labelcolor="C3"); a2b.set_ylabel("phase (rad)", color="C0")
a2b.tick_params(axis="y", labelcolor="C0")
fig.tight_layout(); plt.show()Guided — coupling regimes¶
Vary the ETM reflectivity r2 and read the on-resonance reflected amplitude. The matrix returns
(r1 - r2)/(1 - r1 r2): positive when under-coupled (r2 < r1), zero at critical
coupling (r2 = r1, all light enters the cavity), negative when over-coupled (r2 > r1 —
the LIGO arm, ETM ≈ perfect). Energy stays conserved throughout.
r2_scan = np.linspace(0.970, 0.9995, 161) # all physical (<= 1), brackets r1
r_on_res = []
for rr in r2_scan:
Mc = np.array([[0, 0, 0, -r1],
[1, 0, 0, 0 ],
[0, -rr, 0, 0],
[0, 0, 1, 0 ]], dtype=complex) # phi = 0 (on resonance)
Ec = np.linalg.solve(np.eye(4, dtype=complex) - Mc, np.array([t1, 0, 0, 0], dtype=complex))
r_on_res.append(r1 + t1 * Ec[3])
r_on_res = np.array(r_on_res)
r_analytic = (r1 - r2_scan) / (1 - r1 * r2_scan) # closed form
assert np.allclose(r_on_res, r_analytic, atol=1e-12), "coupling closed-form mismatch"
ic = np.argmin(np.abs(r_on_res))
print(f"critical coupling near r2 = {r2_scan[ic]:.5f} (r1 = {r1:.5f}); "
f"min |r_cav| = {abs(r_on_res[ic]):.2e}")
fig, ax = plt.subplots(figsize=(6.5, 4))
ax.plot(r2_scan, r_on_res.real, lw=2)
ax.axvline(r1, color="C3", ls="--", label="$r_2 = r_1$ (critical)")
ax.axhline(0, color="grey", lw=0.8)
ax.set(xlabel="ETM amplitude reflectivity $r_2$", ylabel="on-resonance $r_\\mathrm{cav}$",
title="Fig 5 — Under / critical / over coupling")
ax.annotate("under-coupled\n($r_2<r_1$)", xy=(0.974, 0.35), fontsize=8, color="C0")
ax.annotate("over-coupled\n($r_2>r_1$)", xy=(0.992, -0.6), fontsize=8, color="C0")
ax.legend(); ax.grid(alpha=0.3)
fig.tight_layout(); plt.show()Part B — SISO feedback loop (guided)¶
The same recipe. Plant P, controller C, nodes (e, u, y):
e = r − y, u = C e, y = P u + d. Edges: y→e (−1), e→u (C), u→y (P).
Your turn: fill in the three non-zero entries of M from the drawing. (The reference is
shown so the notebook runs end-to-end; for live delivery these can be blanked.)
def loop_matrix(C, P):
'''Coupling matrix M for nodes (e, u, y), by inspection of the loop drawing.'''
return np.array([[0, 0, -1], # e gets -1 * y
[C, 0, 0], # u gets C * e
[0, P, 0]], dtype=complex) # y gets P * u
def loop_solve(C, P, r=1.0, d=0.0):
A = np.eye(3, dtype=complex) - loop_matrix(C, P)
return np.linalg.solve(A, np.array([r, 0, d], dtype=complex)) # (e, u, y)
# Closed-loop transfer functions must drop out of the inverse with no extra algebra.
C0, P0 = 4.0 + 1j, 3.0 - 0.5j
e, u, y = loop_solve(C0, P0, r=1.0, d=0.0)
S = 1 / (1 + P0 * C0) # sensitivity
T = P0 * C0 / (1 + P0 * C0) # complementary sensitivity
assert np.isclose(y, T, rtol=1e-12) and np.isclose(e, S, rtol=1e-12), "S/T mismatch"
# Read S and T straight off (I - M)^{-1}.
inv = np.linalg.inv(np.eye(3, dtype=complex) - loop_matrix(C0, P0))
assert np.isclose(inv[2, 0], T) and np.isclose(inv[0, 0], S), "inverse entries mismatch"
print("closed-loop S, T read directly off the matrix inverse ✓")
print(f" the cavity's round-trip gain r1*r2 = {r1*r2:.5f} plays the role of the loop gain PC")Part C — Transfer functions in the audio-sideband approximation (independent)¶
A slow modulation at audio frequency Ω puts sidebands on the carrier at ±Ω. In the linear
regime each sideband sees the same matrix with frequency-dependent edges — so a transfer
function is just (I − M)^{-1} re-evaluated at s = iΩ.
Independent task: the carrier sits on resonance (φ0 = 0). A sideband at offset Ω sees an
extra round-trip phase 2ΩL/c, so the round-trip gain is r1 r2 e^{2iΩL/c}. Compute the
normalized circulating-field response H(Ω) = E(Ω)/E(0), find the −3 dB frequency, and compare to
f_c = FSR/(2·𝓕).
def cavity_response(f):
Omega = 2 * np.pi * np.asarray(f, dtype=float)
g = r1 * r2 * np.exp(2j * Omega * L / cv.C_SI) # carrier on resonance
return 1.0 / (1 - g)
f = np.linspace(0.0, 400.0, 1601)
H = cavity_response(f) / cavity_response(0.0)
mag = np.abs(H)
# numerically locate where |H| crosses 1/sqrt(2) (half power) -> the cavity pole
target = 1 / np.sqrt(2)
i = np.argmax(mag < target)
f_pole_meas = np.interp(target, [mag[i], mag[i-1]], [f[i], f[i-1]])
rel_err = abs(f_pole_meas - f_c) / f_c
print(f"measured pole = {f_pole_meas:.2f} Hz | analytic f_c = {f_c:.2f} Hz | rel.err = {rel_err:.3%}")
assert rel_err < 0.02, "measured cavity pole disagrees with FSR/(2F)"
fig, ax = plt.subplots(figsize=(6.5, 4))
ax.semilogx(f[1:], 20 * np.log10(mag[1:]), lw=2)
ax.axvline(f_c, color="C3", ls="--", label=f"$f_c$ = {f_c:.1f} Hz")
ax.axhline(-3, color="grey", ls=":", lw=1)
ax.set(xlabel="audio frequency Ω/2π (Hz)", ylabel="|H| (dB)",
title="Fig 3 — Cavity frequency response (single pole)")
ax.legend(); ax.grid(which="both", alpha=0.3)
fig.tight_layout(); plt.show()Why near-resonance is numerically delicate¶
As the loop gain r1 r2 e^{2iφ} → 1, the matrix (I − M) becomes nearly singular — the same
reason a high-gain feedback loop is sensitive. The condition number quantifies it.
phi_c = np.linspace(-0.05, 0.05, 301)
cond = np.array([np.linalg.cond(np.eye(4, dtype=complex) - fp_coupling_matrix(p, r1, r2))
for p in phi_c])
fig, ax = plt.subplots(figsize=(6.5, 4))
ax.semilogy(phi_c * FSR / np.pi, cond, lw=2, color="C4")
ax.set(xlabel="detuning (Hz)", ylabel=r"cond$(I-M)$",
title="Fig 4 — Conditioning peaks on resonance")
ax.grid(which="both", alpha=0.3)
fig.tight_layout(); plt.show()Supplementary — the loop’s Bode response (same s = iΩ trick)¶
Give the plant a pole and the controller an integrator, evaluate M(s) at s = iΩ, and read
S,T off the inverse at every frequency. Same recipe, now frequency-resolved.
def PC_of_s(s, g=300.0, wp=2*np.pi*40, ki=2*np.pi*5):
P = g / (1 + s / wp) # single-pole plant
C = 1 + ki / s # proportional + integrator
return C, P
fb = np.logspace(-1, 3, 241)
s = 1j * 2 * np.pi * fb
Sm, Tm = [], []
for sk in s:
Ck, Pk = PC_of_s(sk)
x = np.linalg.solve(np.eye(3, dtype=complex) - loop_matrix(Ck, Pk),
np.array([1, 0, 0], dtype=complex))
Sm.append(x[0]); Tm.append(x[2])
Sm, Tm = np.abs(Sm), np.abs(Tm)
fig, ax = plt.subplots(figsize=(6.5, 4))
ax.loglog(fb, Tm, lw=2, label="|T| (r→y)")
ax.loglog(fb, Sm, lw=2, label="|S| (d→y)")
ax.set(xlabel="frequency (Hz)", ylabel="magnitude", title="Supplementary — loop S and T")
ax.legend(); ax.grid(which="both", alpha=0.3)
fig.tight_layout(); plt.show()Explain (write it down)¶
In 2–3 sentences each:
What is identical between solving the cavity and solving the loop? What is different?
You doubled
T1in Part A’s prediction — did the pole and the buildup move the way you predicted? Why?
Exit prompt¶
If the cavity round-trip gain is the loop gain, what does “resonance” correspond to in feedback language — and why does that make high-finesse cavities hard to control quickly?
Toward MIMO / a full interferometer¶
Replace the scalars P, C with matrices and the optical nodes with vectors: (I − M)^{-1}
becomes a closed-loop response matrix. A Michelson with arm cavities is just a bigger graph →
bigger M, identical recipe. That is session 2 (MIMO) and session 3 (interferometer).