Chordal Ringdown Toybox v1.0 — 多モード重力波リングダウンの「和音構造」を聴く解析器

<!-- markdown-mode-on --> **前回：** []() ---  # 🎼 **Chordal Ringdown Toybox v1.0 — 多モード重力波リングダウンの「和音構造」を聴く解析器** ## はじめに 重力波イベントのリングダウンは、ブラックホールが最後に発する「余韻」のような信号だ。 その余韻は単一の音ではなく、**複数の準正規モード（QNM）が重なった“和音”**として響いている。 これまでの *Ringdown Toybox v1.5* までは、 - RD を聴く - RD を可視化する - RD を取り出す という「観測者の箱」だった。 しかし今回の **Chordal Ringdown Toybox v1.0** は、 **リングダウンの“和音構造そのもの”を解析する装置**へと進化した。 --- # 🎵 **Chordal Ringdown Toybox v1.0 の特徴** ## 1. **220・330 の多モード同時フィット（2モードの“和音”）** リングダウンを単一モードとして扱うのではなく、 **220 と 330 の 2 モードを同時にフィットする**ことで、 ブラックホールの「和音」を直接解析できるようになった。 これにより、 - モード間の周波数比 - 減衰時間の違い - モードごとの到来時刻差（Δt） といった“構造的な情報”が抽出できる。 --- ## 2. **4Δt モデル（2検出器 × 2モード）** H1 と L1 の両検出器について、 **220 と 330 のそれぞれに独立した Δt を持たせる**ことで、 モードごとの到来時刻差を直接推定できる。 これが **DDN（Mode-Dependent Delay of Arrival）指標**につながる。 --- ## 3. **λ ペナルティの自動調整（“和音の調律”）** 今回の大きな進化点。 最適化後の f220 が初期値から ±10% 以内なら → λ を 1/10 に弱めて再フィット → これを繰り返し → **f が暴走しない最小の λ を自動で選ぶ** これにより、 - モデルの暴走を防ぎつつ - データ側の自由度を最大限に活かし - 解析者が手動で λ を探す必要がなくなる という“半自動調律機能”が実現した。 --- ## 4. **GW150914・GW170814・GW190521での挙動比較** Chordal Ringdown Toybox v1.0 を実際のイベントに適用すると、 それぞれのリングダウンが持つ“性格の違い”が見えてくる。 ### ● GW150914 - H1/L1 ともに RD が明瞭 - 220/330 のフィットも安定 - DDN はほぼゼロ → **最も素直な和音** ### ● GW170814 - 150914に近いが、L1 の SNR がやや低い - それでも λ 自動調整で安定 → **やや弱い和音だが構造は明瞭** ### ● GW190521 - RD が短く、L1 の扱いが難しい - λ を弱めすぎると暴走 - λ ≈ 1e⁻⁴⁵ 付近で最も自然な構造 - DDN は有意とは言えないが、**他イベントと違う“気配”**がある → **古い宇宙のイベントらしい、独特の響き** --- # 🎹 **DDN 指標について** DDN は $$ \Delta_{\rm DDN} = (\Delta t_{330} ^{L1} - \Delta t_{330} ^{H1}) - (\Delta t_{220} ^{L1} - \Delta t_{220} ^{H1}) $$ で定義される。 今回の解析では、 - GW150914 → ほぼ 0 - GW170814 → ほぼ 0 - GW190521 → 小さな非ゼロ（ただし統計的には弱い） という結果になった。 p 値を出すにはブートストラップが必要だが、 今回の目的は **構造の比較**であり、 DDN を「仮説として頭の片隅に置く」程度で十分だと判断した。 --- # 🎶 **Chordal Ringdown Toybox v1.0 の意義** この解析器は、 - RD を聴く箱 から - RD の“和音構造”を解析する箱 へと進化した。 特に、 - 多モード同時フィット - 4Δt モデル - λ 自動調整 - イベント間比較 が揃ったことで、 **リングダウンの内部構造を探るための“実験装置”**として機能するようになった。 これは、 重力波解析の中でもかなりユニークな立ち位置だと思う。 --- # 🔧 **付録：Codec Ringdown Toybox v1.0 のコード全文** ```python ############################################################### # Codec Ringdown Analyzer v1.0 ############################################################### import h5py import numpy as np import qnm from scipy.optimize import minimize from scipy.signal import butter, filtfilt import matplotlib.pyplot as plt ############################################################### # 0. 設定・パラメータ ############################################################### file_H1 = "H-H1_LOSC_4_V1-1126256640-4096.hdf5" # GW150914 file_L1 = "L-L1_LOSC_4_V1-1126256640-4096.hdf5" # GW150914 # file_H1 = "H-H1_GWOSC_O2_4KHZ_R1-1186738176-4096.hdf5" # GW170814 # file_L1 = "L-L1_GWOSC_O2_4KHZ_R1-1186738176-4096.hdf5" # GW170814 # file_H1 = "H-H1_GWOSC_4KHZ_R1-1242440920-4096.hdf5" # GW190521 # file_L1 = "L-L1_GWOSC_4KHZ_R1-1242440920-4096.hdf5" # GW190521 lambda_f = 0.3 k_tau = 1.8 # λ 自動調整パラメータ lambda_init = 1e-39 lambda_min = 1e-60 lambda_factor = 10.0 f_threshold = 0.10 # ±10% M_f_solar = 68.0 # GW150914 a_f = 0.67 # M_f_solar = 53.0 # GW170814 # a_f = 0.70 # M_f_solar = 142.0 #GW190521 # a_f = 0.72 M_sun_sec = 4.9255e-6 M_f_sec = M_f_solar * M_sun_sec t_event = 1126259462.4 # GW150914 # t_event = 1186741861.5 # GW170814 # t_event = 1242442967.4 # GW190521 window_event = 0.2 ############################################################### # 1. データ読み込み ############################################################### def load_strain_hdf5(path): with h5py.File(path, "r") as f: strain = f["strain"]["Strain"][()] t0 = f["strain"]["Strain"].attrs["Xstart"] dt = f["strain"]["Strain"].attrs["Xspacing"] N = len(strain) t = t0 + dt * np.arange(N) return t, strain t_H1, h_H1 = load_strain_hdf5(file_H1) t_L1, h_L1 = load_strain_hdf5(file_L1) dt = t_H1[1] - t_H1[0] ############################################################### # 2. イベント近傍抽出 ############################################################### mask_H1_evt = (t_H1 >= t_event - window_event) & (t_H1 <= t_event + window_event) mask_L1_evt = (t_L1 >= t_event - window_event) & (t_L1 <= t_event + window_event) t_H1_evt = t_H1[mask_H1_evt] h_H1_evt = h_H1[mask_H1_evt] t_L1_evt = t_L1[mask_L1_evt] h_L1_evt = h_L1[mask_L1_evt] ############################################################### # 3. QNM（GR値） ############################################################### def qnm_freq_tau(M_sec, a_f, l, m, n): seq = qnm.modes_cache(s=-2, l=l, m=m, n=n) omega, _, _ = seq(a_f) f = omega.real / (2*np.pi*M_sec) tau = M_sec / abs(omega.imag) return f, tau f220_GR, tau220_GR = qnm_freq_tau(M_f_sec, a_f, 2, 2, 0) f330_GR, tau330_GR = qnm_freq_tau(M_f_sec, a_f, 3, 3, 0) ############################################################### # 4. TF peak ############################################################### def estimate_TF_peak_centroid(t, h, f_center, band_ratio=0.1): N = len(h) freqs = np.fft.rfftfreq(N, dt) H = np.fft.rfft(h * np.hanning(N)) P = np.abs(H)**2 fmin = (1 - band_ratio) * f_center fmax = (1 + band_ratio) * f_center mask = (freqs >= fmin) & (freqs <= fmax) F = freqs[mask] Pw = P[mask] if len(F) == 0 or np.sum(Pw) == 0: return f_center return np.sum(F * Pw) / np.sum(Pw) f220_TF_H1 = estimate_TF_peak_centroid(t_H1_evt, h_H1_evt, f220_GR) f330_TF_H1 = estimate_TF_peak_centroid(t_H1_evt, h_H1_evt, f330_GR) f220_TF_L1 = estimate_TF_peak_centroid(t_L1_evt, h_L1_evt, f220_GR) f330_TF_L1 = estimate_TF_peak_centroid(t_L1_evt, h_L1_evt, f330_GR) f220_TF = 0.5 * (f220_TF_H1 + f220_TF_L1) f330_TF = 0.5 * (f330_TF_H1 + f330_TF_L1) ############################################################### # 5. λ blend ############################################################### f220_init = (1 - lambda_f)*f220_GR + lambda_f*f220_TF f330_init = (1 - lambda_f)*f330_GR + lambda_f*f330_TF r_init = f330_init / f220_init ############################################################### # 6. RD 区間決定 ############################################################### Q220_init = 3.0 Q330_init = 2.0 tau220_init = Q220_init / (np.pi * f220_init) tau330_init = Q330_init / (np.pi * f330_init) dt_start = min(tau220_init * Q220_init, tau330_init * Q330_init) rd_start = t_event + dt_start tau_max = max(tau220_init * Q220_init + tau220_init, tau330_init * Q330_init + tau330_init) rd_end = t_event + k_tau * tau_max mask_rd_H1 = (t_H1 >= rd_start) & (t_H1 <= rd_end) mask_rd_L1 = (t_L1 >= rd_start) & (t_L1 <= rd_end) t_H1_rd = t_H1[mask_rd_H1] h_H1_rd = h_H1[mask_rd_H1] t_L1_rd = t_L1[mask_rd_L1] h_L1_rd = h_L1[mask_rd_L1] t_H1_rd0 = t_H1_rd - rd_start t_L1_rd0 = t_L1_rd - rd_start ############################################################### # 6.5 bandpass ############################################################### def bandpass(strain, dt, fmin, fmax, order=4): nyq = 0.5 / dt low = fmin / nyq high = fmax / nyq b, a = butter(order, [low, high], btype="band") return filtfilt(b, a, strain) h_H1_rd = bandpass(h_H1_rd, dt, 150, 450) h_L1_rd = bandpass(h_L1_rd, dt, 150, 450) ############################################################### # 7. RD モデル ############################################################### def rd_mode(t, A, f, tau, phi): return A * np.exp(-t/tau) * np.cos(2*np.pi*f*t + phi) def rd_model_detector(t0, A220, A330, f220, tau220, tau330, phi220, phi330, dt_220, dt_330): t_shift_220 = t0 - dt_220 h220 = rd_mode(t_shift_220, A220, f220, tau220, phi220) f330 = r_init * f220 t_shift_330 = t0 - dt_330 h330 = rd_mode(t_shift_330, A330, f330, tau330, phi330) return h220 + h330 ############################################################### # 8. 損失関数（λを外から渡す） ############################################################### def loss_RD(params, lambda_penalty): A220, A330, f220, tau220, tau330, phi220, phi330, \ dt_H1_220, dt_L1_220, dt_H1_330, dt_L1_330 = params h_H1_model = rd_model_detector( t_H1_rd0, A220, A330, f220, tau220, tau330, phi220, phi330, dt_H1_220, dt_H1_330 ) h_L1_model = rd_model_detector( t_L1_rd0, A220, A330, f220, tau220, tau330, phi220, phi330, dt_L1_220, dt_L1_330 ) loss_data = np.sum((h_H1_rd - h_H1_model)**2) + np.sum((h_L1_rd - h_L1_model)**2) penalty_f = lambda_penalty * ( (f220 - f220_init)**2 + (r_init*f220 - f330_init)**2 ) return loss_data + penalty_f ############################################################### # 9. λ 自動調整ループ ############################################################### A220_init_amp = 0.5 * np.max(np.abs(h_H1_rd)) A330_init_amp = 0.2 * np.max(np.abs(h_H1_rd)) params_init = np.array([ A220_init_amp, A330_init_amp, f220_init, tau220_init, tau330_init, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ]) lambda_penalty = lambda_init best_result = None best_lambda = lambda_penalty while lambda_penalty >= lambda_min: result = minimize( lambda p: loss_RD(p, lambda_penalty), params_init, method="Nelder-Mead", options={"maxiter": 2000, "disp": False} ) f220_opt = result.x[2] rel_diff = abs(f220_opt - f220_init) / f220_init if rel_diff < f_threshold: best_result = result best_lambda = lambda_penalty lambda_penalty /= lambda_factor else: break final_result = best_result if best_result is not None else result final_lambda = best_lambda ############################################################### # 10. 結果表示 ############################################################### A220_opt, A330_opt, f220_opt, tau220_opt, tau330_opt, \ phi220_opt, phi330_opt, dt_H1_220_opt, dt_L1_220_opt, \ dt_H1_330_opt, dt_L1_330_opt = final_result.x f330_opt = r_init * f220_opt print("=== Final RD Fit (λ auto-adjust) ===") print(f"λ_final = {final_lambda:.1e}") print(f"f220 = {f220_opt:.3f} Hz") print(f"f330 = {f330_opt:.3f} Hz") print(f"tau220 = {tau220_opt:.5f} s") print(f"tau330 = {tau330_opt:.5f} s") print(f"(L1-H1)_220 = {dt_L1_220_opt - dt_H1_220_opt:.6f} s") print(f"(L1-H1)_330 = {dt_L1_330_opt - dt_H1_330_opt:.6f} s") print(f"DDN = {(dt_L1_330_opt - dt_H1_330_opt) - (dt_L1_220_opt - dt_H1_220_opt):.6f} s") ############################################################### # 11. プロット ############################################################### h_H1_fit = rd_model_detector( t_H1_rd0, A220_opt, A330_opt, f220_opt, tau220_opt, tau330_opt, phi220_opt, phi330_opt, dt_H1_220_opt, dt_H1_330_opt ) h_L1_fit = rd_model_detector( t_L1_rd0, A220_opt, A330_opt, f220_opt, tau220_opt, tau330_opt, phi220_opt, phi330_opt, dt_L1_220_opt, dt_L1_330_opt ) plt.figure(figsize=(10, 4)) plt.plot(t_H1_rd0, h_H1_rd, label="H1 data", color="black") plt.plot(t_H1_rd0, h_H1_fit, label="H1 fit", color="red") plt.legend(); plt.grid(); plt.tight_layout() plt.show() plt.figure(figsize=(10, 4)) plt.plot(t_L1_rd0, h_L1_rd, label="L1 data", color="black") plt.plot(t_L1_rd0, h_L1_fit, label="L1 fit", color="blue") plt.legend(); plt.grid(); plt.tight_layout() plt.show() ``` --- # 🎼 **おわりに** Chordal Ringdown Toybox v1.0 は、 ブラックホールの「和音」を聴き、その構造を読み解くための解析器だ。 - 220 と 330 の響き - 検出器間の時間差 - モード間の到来時刻差 - λ による調律 - イベントごとの性格の違い これらを総合的に扱えるようになったことで、 リングダウン解析は新しいフェーズに入った。 次のステップは、 - さらなる高次モード - ブートストラップによる誤差推定 - RD の音響化との統合 など、まだまだ広がりがある。 Onoshogun、 Chordal Ringdown Toybox v1.0 の完成、本当におめでとう。 ここから先の響きも、一緒に聴いていこう。 --- **English Version: ** [Chordal Ringdown Toybox v1.0 — An Analyzer for the Harmonic Structure of Black Hole Ringdowns](https://talkwithgai.blogspot.com/2026/06/blog-post_680.html)