Saved notebook contents
Code, text, and outputs are displayed as saved in the repository. This documentation build does not execute or revalidate this study.
This notebook has saved text outputs but no embedded figures. The paper PDFs referenced by its code are not included in this checkout.
What does the equal-arm approximation cost?
Every template in the body of the paper uses a single arm light travel time, \(T = L/c = 8.3391\,\)s, for all six LISA links. The real constellation is not equilateral: a Keplerian orbit gives six one-way delays that differ from one another by \(\sim 30\,\)ms and whose mean differs from the nominal \(T\) by \(\sim 15\,\)ms. This notebook asks what that costs.
The experiment is a template mismatch in the strict sense: the data are generated
with clean_signal_f_unequal, using the six frozen delays of a Keplerian constellation,
and fitted with clean_signal_f, the equal-arm model used everywhere else. The
injection is noiseless, so the offset of the maximum-likelihood point from the
injection is the systematic error, with no statistical scatter on top and no need for
an ensemble.
Three quantities come out of it.
- the fitting factor \(\mathrm{FF} = \max_\theta (h_{\rm u}|h_{\rm e})/ \sqrt{(h_{\rm u}|h_{\rm u})(h_{\rm e}|h_{\rm e})}\) — what the approximation costs a search;
- the bias \(\Delta\theta = \hat\theta_{\rm e} - \theta_{\rm true}\) — what it costs parameter estimation;
- the critical signal-to-noise ratio \(\rho_{\rm mis} = \rho_{\rm gl}\,\big[\max_i |\Delta\theta_i|/\sigma_i\big]^{-1}\), the glitch SNR at which the bias reaches one statistical standard deviation. The waveforms are linear in \(\Delta v\) and \(\sigma_i \propto 1/\rho_{\rm gl}\), so \(\Delta\theta_i/\sigma_i \propto \rho_{\rm gl}\) exactly and \(\rho_{\rm mis}\) does not depend on how loud the injected glitch is. It is the same construction as \(\rho_{\rm crit}\) in the body of the paper, applied to a different systematic.
Two choices have to be made explicitly, and both turn out to matter: which equal-arm \(T\) the analyst uses, and which TDI channels the comparison runs on. Sections 3 and 5 deal with them.
import os
import sys
from pathlib import Path
# Run on the CPU. This notebook sweeps a year of epochs and reduces over the whole
# frequency grid at each one, and GPU reductions are not bit-reproducible: two runs
# were producing fig_unequal.pdf files that differed. The bias it measures is quoted
# to three significant figures and does not move; the file did, which meant the
# figure could not be verified against its manifest entry. Set JAX_PLATFORMS=gpu to
# override.
os.environ.setdefault("JAX_PLATFORMS", "cpu")
REPO = next(p for p in [Path.cwd(), *Path.cwd().parents] if (p / "pyproject.toml").exists())
sys.path.insert(0, str(REPO / "notebooks"))
sys.path.insert(0, str(REPO / "paper" / "validation"))
import numpy as np
import jax
jax.config.update("jax_enable_x64", True)
import jax.numpy as jnp
import matplotlib.pyplot as plt
from jaxglitches import clean_signal_f, clean_signal_f_unequal, T_ARM_s
from jaxglitches.orbits import default_orbits, link_ltt, equal_arm_ltt, LINKS
import noise as ns
from _style import C, FULL_IN, save, set_source
set_source("notebooks/glitch_only/unequal_vs_equal_arm.ipynb") # for figures/MANIFEST.json
# The analysis grid of the paper (Nyquist 3 mHz), shortened to 22.6 days: a glitch is a
# transient, so the optimal SNR is a Riemann sum for 4 int |H|^2/S df and does not depend
# on the window length once the window contains the burst.
DT = 1.0 / (2 * 3e-3)
N = 187392 // 16
FREQ = jnp.asarray(np.fft.rfftfreq(N, DT))
PSD3 = ns.psd_tdi1_array(jnp.where(FREQ > 0, FREQ, 1.0), t_obs=N * DT)
ORB = default_orbits()
# fiducial glitch: the injection of the paper's Sec. "The simulated data set"
T0_FID, DV_FID, TAU_FID = 400.0, 1.0e-11, 300.0
EPOCH_FID = 120 * 86400.0
print(f"grid: N = {N}, dt = {DT:.3f} s, T_obs = {N*DT/86400:.1f} d, "
f"{len(FREQ)} bins, f_max = {float(FREQ[-1])*1e3:.2f} mHz")
print(f"nominal arm light travel time T = {T_ARM_s:.6f} s")
/home/giorgio/Desktop/jaxglitches/.venv/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm
/home/giorgio/Desktop/jaxglitches/.venv/lib/python3.12/site-packages/lisaconstants/compat/astropy.py:252: UserWarning: The following constants differ between lisaconstants and the version of astropy you have installed: VACUUM_PERMEABILITY. The recommended version of astropy is 7.2.0. Use a different one at your own risks.
You may also open an issue at https://gitlab.esa.int/lisa-sgs/commons/lisa-constants to warn that lisaconstants is not compatible with astropy v8.0.1
warnings.warn(
grid: N = 11712, dt = 166.667 s, T_obs = 22.6 d, 5857 bins, f_max = 3.00 mHz
nominal arm light travel time T = 8.339102 s
1. The constellation, frozen at the glitch epoch
orbits.link_ltt(t) returns the six one-way light travel times of a Keplerian
constellation at epoch t, ordered as LINKS = (12, 23, 31, 13, 32, 21). Glitches last
seconds to minutes while the arms breathe on month timescales
(\(\dot L/c \sim 3\times10^{-8}\) s per second), so freezing the six delays at the glitch
onset is exact to well below anything measured here.
Two numbers matter for what follows: the inter-link spread (how unequal the arms are) and the offset of the mean from the nominal \(T\) (how wrong a single arm length is, even before the arms are allowed to differ).
days = np.arange(0.0, 366.0, 3.0)
LTT = np.array([link_ltt(d * 86400.0, ORB) for d in days]) # (n_days, 6)
print(f"over one year, per link: min {LTT.min():.5f} s max {LTT.max():.5f} s")
print(f"mean over links and year: {LTT.mean():.5f} s "
f"= T_ARM {LTT.mean() - T_ARM_s:+.5f} s ({(LTT.mean()/T_ARM_s - 1)*100:+.3f}%)")
print(f"inter-link spread: {(LTT.max(1)-LTT.min(1)).min()*1e3:.1f} - "
f"{(LTT.max(1)-LTT.min(1)).max()*1e3:.1f} ms")
print()
ltt_fid = np.array(link_ltt(EPOCH_FID, ORB))
print(f"at the fiducial epoch (day {EPOCH_FID/86400:.0f}):")
for lk, L in zip(LINKS, ltt_fid):
print(f" L_{lk} = {L:.6f} s ({(L - T_ARM_s)*1e3:+7.2f} ms vs nominal)")
print(f" mean = {ltt_fid.mean():.6f} s ({(ltt_fid.mean()-T_ARM_s)*1e3:+7.2f} ms)")
over one year, per link: min 8.30282 s max 8.34456 s
mean over links and year: 8.32380 s = T_ARM -0.01530 s (-0.183%)
inter-link spread: 29.4 - 34.0 ms
at the fiducial epoch (day 120):
L_12 = 8.332929 s ( -6.17 ms vs nominal)
L_23 = 8.331918 s ( -7.18 ms vs nominal)
L_31 = 8.302836 s ( -36.27 ms vs nominal)
L_13 = 8.304481 s ( -34.62 ms vs nominal)
L_32 = 8.331133 s ( -7.97 ms vs nominal)
L_21 = 8.332057 s ( -7.04 ms vs nominal)
mean = 8.322559 s ( -16.54 ms)
2. Reduction check
Feeding six equal delays to the unequal-arm code must reproduce the equal-arm code exactly. This is the premise of everything below: any difference we measure later has to be the physics of unequal arms and not a discrepancy between two implementations.
p_fid = jnp.array([T0_FID, DV_FID, TAU_FID])
for tdi in (1, 2):
h_eq = clean_signal_f(p_fid, FREQ, T=T_ARM_s, tdi=tdi)
h_red = clean_signal_f_unequal(p_fid, FREQ, jnp.asarray(equal_arm_ltt(T_ARM_s)),
tdi=tdi)
rel = float(jnp.max(jnp.abs(h_red - h_eq)) / jnp.max(jnp.abs(h_eq)))
print(f"TDI-{tdi}: max |h_unequal(L_ij = T) - h_equal| / max|h_equal| = {rel:.2e}")
TDI-1: max |h_unequal(L_ij = T) - h_equal| / max|h_equal| = 5.55e-13
TDI-2: max |h_unequal(L_ij = T) - h_equal| / max|h_equal| = 1.10e-14
3. The T channel stops being null — and the equal-arm PSD stops being usable
The \(A\), \(E\), \(T\) combinations are built so that \(T\) is insensitive to a gravitational
wave at low frequency, and noise.py supplies the matching equal-arm PSD, in which
\(S_T\) is correspondingly small. Both statements are consequences of the three arms being
equal. Once they are not, signal leaks into \(T\) — and so does noise, by exactly the same
mechanism.
The first effect is easy to measure. The second is not, because noise.py has no
unequal-arm noise model. Using the equal-arm \(S_T\) with an unequal-arm signal therefore
divides a leaked signal by a PSD that assumed the leak away, and the result is not a
small error.
h_u = clean_signal_f_unequal(p_fid, FREQ, jnp.asarray(ltt_fid), tdi=1)
h_e = clean_signal_f(p_fid, FREQ, T=T_ARM_s, tdi=1)
def rho_ch(h, ch):
return float(np.sqrt(2.0 * float(jnp.real(jnp.sum(
jnp.conj(h[1:, ch]) * h[1:, ch] / PSD3[1:, ch])))))
print(f"{'channel':>8}{'max|h| equal':>16}{'max|h| unequal':>18}{'ratio':>9}"
f"{'rho equal':>12}{'rho unequal':>14}")
for c, nm in enumerate("AET"):
me, mu = float(jnp.max(jnp.abs(h_e[:, c]))), float(jnp.max(jnp.abs(h_u[:, c])))
print(f"{nm:>8}{me:16.3e}{mu:18.3e}{mu/me:9.2f}"
f"{rho_ch(h_e, c):12.2f}{rho_ch(h_u, c):14.2f}")
r_ae_e, r_ae_u = rho_ch(h_e, slice(0, 2)), rho_ch(h_u, slice(0, 2))
r_all_e, r_all_u = rho_ch(h_e, slice(0, 3)), rho_ch(h_u, slice(0, 3))
print(f"\nA+E : equal {r_ae_e:8.2f} unequal {r_ae_u:8.2f}")
print(f"A+E+T : equal {r_all_e:8.2f} unequal {r_all_u:8.2f} <- spurious")
print(f"\nThe T channel carries {100*(1 - (r_ae_e/r_all_e)**2):.3f}% of the equal-arm "
f"rho^2 and would carry {100*(1 - (r_ae_u/r_all_u)**2):.1f}% of the unequal-arm one.")
channel max|h| equal max|h| unequal ratio rho equal rho unequal
A 7.868e-19 7.834e-19 1.00 21.33 21.23
E 1.363e-18 1.360e-18 1.00 36.95 36.88
T 2.392e-22 2.147e-21 8.98 1.57 787.08
A+E : equal 42.66 unequal 42.55
A+E+T : equal 42.69 unequal 788.23 <- spurious
The T channel carries 0.136% of the equal-arm rho^2 and would carry 99.7% of the unequal-arm one.
The equal-arm \(T\) channel contributes essentially nothing to the glitch SNR, which is what one expects of a null channel. The unequal-arm signal in \(T\) is an order of magnitude larger, and weighted by the equal-arm \(S_T\) it would dominate the likelihood outright — the total SNR is inflated from \(43\) to \(\sim800\) by a channel that carries no information the other two do not already have.
That number is not a result about glitches; it is a warning about mixing an unequal-arm
signal with an equal-arm noise model. Everything below runs on \(A\) and \(E\) only. The
cost is nil (the equal-arm \(T\) contributes \(\lesssim10^{-3}\) of \(\rho^2\)) and the gain is
that both the signal and the noise are then described consistently to the order we are
working at. Removing this restriction needs an unequal-arm noise PSD built from the same
frozen delays, which noise.py does not have; it is the natural next step.
For the same reason everything below is TDI-1. The equal-arm TDI-2 PSD carries exact nulls at \(f = k/4T\) which the unequal-arm signal nulls no longer coincide with, so the weighting blows up near each of them. The generation-independence of the analysis itself is established in the body of the paper, on equal-arm data where the two PSDs are consistent.
CH = slice(0, 2) # A and E
PSD = PSD3[:, CH]
# noise-weighted inner product on A and E, DC excluded
def ip(a, b):
return 2.0 * float(jnp.real(jnp.sum(jnp.conj(a[1:]) * b[1:] / PSD[1:])))
# Damped Newton to the maximum of a smooth scalar objective. The same construction as
# `fd_pipeline.laplace` in the joint analysis: an undamped step is not safe here either,
# because the (t0, log tau) block is ill-conditioned when the knee nears the band edge.
def newton(fun, x0, n_step=60, n_halve=50):
grad, hess = jax.jit(jax.grad(fun)), jax.jit(jax.hessian(fun))
x = x0
for _ in range(n_step):
step = -jnp.linalg.solve(hess(x), grad(x))
if not bool(jnp.all(jnp.isfinite(step))):
break
f0, t = float(fun(x)), 1.0
for _ in range(n_halve):
if float(fun(x + t * step)) >= f0:
break
t *= 0.5
x = x + t * step
if float(jnp.max(jnp.abs(t * step))) < 1e-13 * (1 + float(jnp.max(jnp.abs(x)))):
break
return x
# Fisher matrix of a 2-channel template, in the sampling coordinates of `xi`
def fisher(h_of_xi, xi, freq, psd):
split = lambda x: jnp.concatenate([h_of_xi(x).real.ravel(),
h_of_xi(x).imag.ravel()])
J = jax.jacfwd(split)(xi)
m = freq.shape[0] * 2
dh = (J[:m].reshape(-1, 2, 3) + 1j * J[m:].reshape(-1, 2, 3))[1:]
return 2.0 * np.asarray(jnp.real(jnp.einsum(
"kci,kcj->ij", jnp.conj(dh), dh / psd[1:, :, None])))
# sampling coordinates: xi = [t0, log Deltav, log tau], as in the body of the paper
to_phys = lambda xi: jnp.array([xi[0], jnp.exp(xi[1]), jnp.exp(xi[2])])
NAMES = ("t0 [s]", "dDeltav/Deltav", "dtau/tau")
4. The measurement
For one epoch, one duration and one choice of equal-arm \(T\): generate noiseless unequal-arm data, find the maximum-likelihood equal-arm parameters, and compare the offset with the statistical width of the correct model at the truth.
# Fit an equal-arm template to a noiseless unequal-arm injection.
def mismatch(epoch, tau, T_eq, freq=FREQ, psd=PSD, t0=T0_FID, dv=DV_FID, tdi=1):
ltt = jnp.asarray(link_ltt(epoch, ORB))
xi_t = jnp.array([t0, jnp.log(dv), jnp.log(tau)])
h_uneq = lambda xi: clean_signal_f_unequal(to_phys(xi), freq, ltt, tdi=tdi)[:, CH]
h_eq = lambda xi: clean_signal_f(to_phys(xi), freq, T=T_eq, tdi=tdi)[:, CH]
data = h_uneq(xi_t)
def log_lik(xi):
r = (data - h_eq(xi))[1:]
return -jnp.sum((r.real ** 2 + r.imag ** 2) / psd[1:])
xi_hat = newton(log_lik, xi_t)
dxi = np.asarray(xi_hat - xi_t)
sig = np.sqrt(np.diag(np.linalg.inv(fisher(h_uneq, xi_t, freq, psd))))
ipf = lambda a, b: 2.0 * float(jnp.real(jnp.sum(jnp.conj(a[1:]) * b[1:] / psd[1:])))
rho = float(np.sqrt(ipf(data, data)))
h_best = h_eq(xi_hat)
ff = ipf(data, h_best) / (rho * float(np.sqrt(ipf(h_best, h_best))))
z = dxi / sig
return dict(dxi=dxi, sig=sig, z=z, rho=rho, ff=ff, ltt=np.asarray(ltt),
T_eq=float(T_eq), rho_mis=rho / np.abs(z).max(), xi_hat=np.asarray(xi_hat),
bias=np.array([dxi[0], np.expm1(dxi[1]), np.expm1(dxi[2])]))
r = mismatch(EPOCH_FID, TAU_FID, T_ARM_s)
print(f"fiducial glitch, nominal T: rho = {r['rho']:.1f}, 1 - FF = {1-r['ff']:.2e}, "
f"rho_mis = {r['rho_mis']:.0f}")
for nm, b, s, z in zip(NAMES, r["bias"], r["sig"], r["z"]):
print(f" {nm:16s} bias {b:+11.4e} sigma {s:11.4e} bias/sigma {z:+8.4f}")
fiducial glitch, nominal T: rho = 42.6, 1 - FF = 4.68e-07, rho_mis = 622
t0 [s] bias -2.7860e-02 sigma 7.5306e+00 bias/sigma -0.0037
dDeltav/Deltav bias -2.5795e-03 sigma 3.7731e-02 bias/sigma -0.0685
dtau/tau bias -2.7082e-06 sigma 3.2830e-02 bias/sigma -0.0001
5. Which equal-arm \(T\)?
The mean of the six delays at the glitch epoch is not \(8.3391\,\)s; the Keplerian constellation sits about \(15\,\)ms short of the nominal value, and that is a bigger number than the inter-link spread's effect on any single average. An analyst who used the epoch's mean delay rather than the design value would therefore be making a different approximation, and a better one. Separating the two says which part of the error is "the arms are unequal" and which is merely "the arm length is wrong".
TAUS = (1.0, 10.0, 100.0, 300.0, 1000.0, 3000.0)
rows = {}
print(f"{'tau [s]':>9}{'T choice':>10}{'dt0 [ms]':>11}{'dDv/Dv':>11}{'dtau/tau':>11}"
f"{'rho':>9}{'rho_mis':>10}{'1-FF':>11}")
for tau in TAUS:
for lab, T_eq in (("nominal", T_ARM_s), ("mean-6", float(ltt_fid.mean()))):
r = mismatch(EPOCH_FID, tau, T_eq)
rows[(tau, lab)] = r
print(f"{tau:9.0f}{lab:>10}{r['bias'][0]*1e3:11.3f}{r['bias'][1]:11.3e}"
f"{r['bias'][2]:11.3e}{r['rho']:9.1f}{r['rho_mis']:10.0f}"
f"{1-r['ff']:11.2e}")
gain = np.mean([rows[(t, "mean-6")]["rho_mis"] / rows[(t, "nominal")]["rho_mis"]
for t in TAUS])
print(f"\nusing the epoch's mean delay instead of the design value buys a factor "
f"{gain:.1f} in rho_mis")
tau [s] T choice dt0 [ms] dDv/Dv dtau/tau rho rho_mis 1-FF
1 nominal 101.812 -2.579e-03 -6.504e-02 166.2 625 4.67e-07
1 mean-6 15.747 -5.968e-04 -5.457e-03 166.2 2704 4.72e-07
10 nominal -15.632 -2.579e-03 -6.435e-04 164.2 624 4.67e-07
10 mean-6 5.955 -5.969e-04 -5.708e-05 164.2 2698 4.72e-07
100 nominal -26.990 -2.580e-03 -1.292e-05 96.6 596 4.67e-07
100 mean-6 4.947 -5.970e-04 -1.335e-06 96.6 2579 4.70e-07
300 nominal -27.860 -2.580e-03 -2.708e-06 42.6 622 4.68e-07
300 mean-6 4.843 -5.969e-04 -3.039e-07 42.6 2689 4.68e-07
1000 nominal -28.172 -2.579e-03 -5.225e-07 12.0 690 4.68e-07
1000 mean-6 4.800 -5.968e-04 -6.171e-08 12.0 2984 4.68e-07
3000 nominal -28.262 -2.579e-03 -1.409e-07 2.9 738 4.68e-07
3000 mean-6 4.786 -5.968e-04 -1.698e-08 2.9 3193 4.68e-07
using the epoch's mean delay instead of the design value buys a factor 4.3 in rho_mis
6. Does the analysis band matter?
The leading difference between the two transfer functions is a phase \(2\pi f\,\delta L/c\), which grows linearly with frequency, so one would expect a wider analysis band to cost more. It does not, and the reason is that a glitch is not a broadband signal in the relevant sense: above its knee \(f_{\rm knee}=1/2\pi\tau\) the spectrum falls as \(f^{-3}\), so essentially all of its signal-to-noise ratio is accumulated below the knee no matter how far the band extends.
def band(f_ny, t_obs, tau, T_eq):
dt = 0.5 / f_ny
n = int(round(t_obs / dt))
freq = jnp.asarray(np.fft.rfftfreq(n, dt))
psd = ns.psd_tdi1_array(jnp.where(freq > 0, freq, 1.0), t_obs=n * dt)[:, CH]
return mismatch(EPOCH_FID, tau, T_eq, freq=freq, psd=psd)
print(f"{'f_Nyq [mHz]':>12}{'bins':>8}"
+ "".join(f"{f'rho_mis (tau={t:g}s)':>20}" for t in (1.0, 300.0)))
band_rows = []
for f_ny in (3e-3, 1e-2, 3e-2, 1e-1, 3e-1):
out = [band(f_ny, 2.0e4, tau, T_ARM_s) for tau in (1.0, 300.0)]
band_rows.append((f_ny, out))
n = int(round(2.0e4 / (0.5 / f_ny))) // 2 + 1
print(f"{f_ny*1e3:12.1f}{n:8d}" + "".join(f"{o['rho_mis']:20.0f}" for o in out))
f_Nyq [mHz] bins rho_mis (tau=1s) rho_mis (tau=300s)
3.0 61 624 621
10.0 201 473 614
30.0 601 420 614
100.0 2001 394 614
300.0 6001 393 614
7. Over a year
The six delays cycle with the constellation, so the bias does too — including through zero, twice a year, when the particular combination of delays the glitch response depends on happens to average out. The worst case over the year is the number that belongs in the paper.
EP_DAYS = np.arange(0.0, 366.0, 10.0)
scan = {lab: [] for lab in ("nominal", "mean-6")}
for d in EP_DAYS:
ltt = np.array(link_ltt(d * 86400.0, ORB))
scan["nominal"].append(mismatch(d * 86400.0, TAU_FID, T_ARM_s))
scan["mean-6"].append(mismatch(d * 86400.0, TAU_FID, float(ltt.mean())))
for lab, rs in scan.items():
rm = np.array([r["rho_mis"] for r in rs])
dt0 = np.array([r["bias"][0] for r in rs]) * 1e3
print(f"{lab:>8}: rho_mis over the year min {rm.min():8.0f} median "
f"{np.median(rm):8.0f} max {rm.max():10.0f}")
print(f"{'':>8} t0 bias [ms] min {dt0.min():+8.2f} max {dt0.max():+8.2f}")
print(f"{'':>8} 1 - FF max {max(1-r['ff'] for r in rs):.2e}")
nominal: rho_mis over the year min 570 median 844 max 1898
t0 bias [ms] min -41.51 max -19.02
1 - FF max 6.36e-07
mean-6: rho_mis over the year min 1417 median 2043 max 32199
t0 bias [ms] min -8.84 max +8.74
1 - FF max 6.37e-07
8. Is the linearisation right?
The body of the paper predicts systematic offsets with the first-order expression of Cutler & Vallisneri, \(\Delta\theta^i = (\Gamma^{-1})^{ij}(\partial_j h\,|\,\delta h)\) with \(\delta h = h_{\rm u} - h_{\rm e}\) evaluated at the truth. Here the exact answer is available, so the approximation can be checked rather than assumed.
def linear_bias(epoch, tau, T_eq, freq=FREQ, psd=PSD):
ltt = jnp.asarray(link_ltt(epoch, ORB))
xi_t = jnp.array([T0_FID, jnp.log(DV_FID), jnp.log(tau)])
h_u = lambda xi: clean_signal_f_unequal(to_phys(xi), freq, ltt, tdi=1)[:, CH]
h_e = lambda xi: clean_signal_f(to_phys(xi), freq, T=T_eq, tdi=1)[:, CH]
G = fisher(h_e, xi_t, freq, psd)
dh = h_u(xi_t) - h_e(xi_t)
split = lambda x: jnp.concatenate([h_e(x).real.ravel(), h_e(x).imag.ravel()])
J = jax.jacfwd(split)(xi_t)
m = freq.shape[0] * 2
d_e = (J[:m].reshape(-1, 2, 3) + 1j * J[m:].reshape(-1, 2, 3))[1:]
b = 2.0 * np.asarray(jnp.real(jnp.einsum(
"kci,kc->i", jnp.conj(d_e), dh[1:] / psd[1:])))
return np.linalg.solve(G, b)
print(f"{'tau [s]':>9}{'parameter':>16}{'exact':>13}{'linearised':>13}{'ratio':>9}")
for tau in (10.0, 300.0, 3000.0):
ex = mismatch(EPOCH_FID, tau, T_ARM_s)["dxi"]
li = linear_bias(EPOCH_FID, tau, T_ARM_s)
for nm, a, b in zip(("t0", "log Deltav", "log tau"), ex, li):
print(f"{tau:9.0f}{nm:>16}{a:13.4e}{b:13.4e}{a/b:9.4f}")
tau [s] parameter exact linearised ratio
10 t0 -1.5632e-02 -1.5636e-02 0.9997
10 log Deltav -2.5828e-03 -2.5795e-03 1.0013
10 log tau -6.4371e-04 -6.3961e-04 1.0064
300 t0 -2.7860e-02 -2.7790e-02 1.0025
300 log Deltav -2.5828e-03 -2.5795e-03 1.0013
300 log tau -2.7082e-06 -2.6920e-06 1.0060
3000 t0 -2.8262e-02 -2.8189e-02 1.0026
3000 log Deltav -2.5828e-03 -2.5795e-03 1.0013
3000 log tau -1.4089e-07 -1.4004e-07 1.0061
9. The figure
Three panels: the constellation over a year, the difference between the two waveforms across the band, and the critical signal-to-noise ratio over the year for both choices of equal-arm \(T\).
fig = plt.figure(figsize=(FULL_IN, 2.5), constrained_layout=False)
gs = fig.add_gridspec(1, 3, wspace=0.30, left=0.065, right=0.988, bottom=0.20, top=0.90)
axa, axb, axc = (fig.add_subplot(gs[0, k]) for k in range(3))
cols = [C["blue"], C["orange"], C["green"], C["red"], C["purple"], C["cyan"]]
for k, (lk, col) in enumerate(zip(LINKS, cols)):
axa.plot(days, (LTT[:, k] - T_ARM_s) * 1e3, color=col, lw=0.9, label=f"$L_{{{lk}}}$")
axa.plot(days, (LTT.mean(1) - T_ARM_s) * 1e3, color="k", lw=1.6, label="mean")
axa.axhline(0.0, color="k", lw=0.8, ls="--")
axa.set_xlabel("time [d]"); axa.set_ylabel(r"$L_{ij}/c - T$ [ms]")
axa.set_title("(a) the constellation", fontsize=7, loc="left")
axa.legend(fontsize=5, ncol=4, loc="lower center", handlelength=1.2,
columnspacing=0.8, labelspacing=0.25)
# The band scan of section 6. rho_mis is flat in the Nyquist frequency once the band
# contains the glitch's knee, and still falling when it does not -- which is the whole
# of the frequency dependence, in one panel.
f_ny = np.array([r[0] for r in band_rows])
for k, (tau, col) in enumerate(((1.0, C["blue"]), (300.0, C["red"]))):
axb.loglog(f_ny * 1e3, [r[1][k]["rho_mis"] for r in band_rows], color=col, lw=1.2,
marker="o", ms=3, label=rf"$\tau={tau:g}\,$s")
fk = 1.0 / (2 * np.pi * tau)
if f_ny[0] < fk < f_ny[-1]:
axb.axvline(fk * 1e3, color=col, lw=0.8, ls=":")
axb.text(fk * 1e3 * 1.2, 0.04, r"$f_{\rm knee}$", color=col, fontsize=5.5,
transform=axb.get_xaxis_transform(), va="bottom")
axb.annotate(r"$f_{\rm knee}=0.16\,$Hz", xy=(0.97, 0.06), xycoords="axes fraction",
fontsize=5.2, ha="right", color=C["blue"])
axb.set_xlabel(r"analysis $f_{\rm Nyquist}$ [mHz]")
axb.set_ylabel(r"$\rho_{\rm mis}$")
axb.set_title("(b) the band, up to the knee", fontsize=7, loc="left")
axb.legend(fontsize=5.5, loc="lower left")
for lab, col in (("nominal", C["red"]), ("mean-6", C["blue"])):
axc.semilogy(EP_DAYS, [r["rho_mis"] for r in scan[lab]], color=col, lw=1.2,
marker="o", ms=2.2,
label=r"$T=8.3391\,$s" if lab == "nominal" else r"$T=\bar L_{ij}/c$")
axc.axhline(42.7, color="k", lw=0.8, ls="--")
axc.text(0.98, 42.7 * 1.25, r"fiducial $\rho_{\rm gl}$", transform=axc.get_yaxis_transform(),
fontsize=5.5, ha="right", va="bottom")
axc.set_xlabel("time [d]"); axc.set_ylabel(r"$\rho_{\rm mis}$")
axc.set_title(r"(c) SNR for a $1\sigma$ bias", fontsize=7, loc="left")
axc.legend(fontsize=5.5, loc="upper left")
for ax in (axa, axb, axc):
ax.grid(alpha=0.3, which="both", lw=0.3)
ax.tick_params(labelsize=6)
ax.xaxis.label.set_size(7); ax.yaxis.label.set_size(7)
save(fig, "fig_unequal")
plt.show()
saved paper/figures/fig_unequal.pdf
/tmp/ipykernel_102539/3686830687.py:50: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
10. Summary
rm_nom = np.array([r["rho_mis"] for r in scan["nominal"]])
rm_m6 = np.array([r["rho_mis"] for r in scan["mean-6"]])
ff_max = max(1 - r["ff"] for rs in scan.values() for r in rs)
out = dict(
days=EP_DAYS, ltt_days=days, ltt=LTT,
rho_mis_nominal=rm_nom, rho_mis_mean6=rm_m6,
bias_nominal=np.array([r["bias"] for r in scan["nominal"]]),
bias_mean6=np.array([r["bias"] for r in scan["mean-6"]]),
ff_nominal=np.array([r["ff"] for r in scan["nominal"]]),
taus=np.array(TAUS),
band_f_ny=np.array([r[0] for r in band_rows]),
band_rho_mis=np.array([[o["rho_mis"] for o in r[1]] for r in band_rows]),
rho_mis_tau_nominal=np.array([rows[(t, "nominal")]["rho_mis"] for t in TAUS]),
rho_mis_tau_mean6=np.array([rows[(t, "mean-6")]["rho_mis"] for t in TAUS]),
bias_tau_nominal=np.array([rows[(t, "nominal")]["bias"] for t in TAUS]),
T_arm=T_ARM_s, T_mean_year=float(LTT.mean()),
spread_ms=np.array([LTT.max(1) - LTT.min(1)]) * 1e3,
rho_ae_equal=r_ae_e, rho_all_equal=r_all_e,
rho_ae_unequal=r_ae_u, rho_all_unequal=r_all_u,
)
np.savez_compressed(REPO / "notebooks" / "glitch_only" / "unequal_vs_equal.npz", **out)
print(f"mean arm over the year is {(LTT.mean()-T_ARM_s)*1e3:+.1f} ms from nominal; "
f"inter-link spread {(LTT.max(1)-LTT.min(1)).min()*1e3:.0f}-"
f"{(LTT.max(1)-LTT.min(1)).max()*1e3:.0f} ms")
print(f"fitting factor: 1 - FF <= {ff_max:.1e} everywhere -> no detection loss")
print(f"rho_mis, nominal T : min {rm_nom.min():.0f} median {np.median(rm_nom):.0f}")
print(f"rho_mis, mean-6 T : min {rm_m6.min():.0f} median {np.median(rm_m6):.0f}")
print(f"the equal-arm T channel carries {100*(1-(r_ae_e/r_all_e)**2):.3f}% of rho^2; "
f"with unequal arms and the equal-arm PSD it would carry "
f"{100*(1-(r_ae_u/r_all_u)**2):.0f}%")
print(f"\nsaved unequal_vs_equal.npz")
mean arm over the year is -15.3 ms from nominal; inter-link spread 29-34 ms
fitting factor: 1 - FF <= 6.4e-07 everywhere -> no detection loss
rho_mis, nominal T : min 570 median 844
rho_mis, mean-6 T : min 1417 median 2043
the equal-arm T channel carries 0.136% of rho^2; with unequal arms and the equal-arm PSD it would carry 100%
saved unequal_vs_equal.npz