Ringdown Analyzer v1.2 → v1.3 重力波リングダウン解析器が「スケール」を取り戻すまで

<!-- markdown-mode-on --> 前のバージョン： [家庭用リングダウン解析器 v1.2 ノイズだらけのデータから、ブラックホールの声を取り出す「シンプルで壊れない」重力波解析モデル](https://talkwithgai.blogspot.com/2026/06/blog-post_21.html) ---  # Ringdown Analyzer v1.2 → v1.3 ## ― 重力波リングダウン解析器が「スケール」を取り戻すまで ― ## はじめに 本記事では、私が開発している軽量リングダウン解析器の **v1.2 → v1.3** の進化をまとめる。 v1.2 は「位相だけで動く」ため、 **質量スケールが決まらない（＝宇宙の“ものさし”が無い）**という根本的な限界を抱えていた。 しかし v1.3 では、**物理的な chirp mass（Mc_phys）を入力するだけで、 リングダウン周波数 f220 が正しいスケールに復元される**。 この変化は、解析器の世界観が **“スケール不定の宇宙” → “物理スケールの宇宙”** へと切り替わる瞬間だった。 --- ## 1. v1.2 の問題：スケールが決まらない宇宙 v1.2 は以下の特徴を持つ： - Hilbert 変換で位相 φ(t) を取り出す - φ(t) の時間微分から瞬時周波数 f(t) を得る - 最後の 1ms の周波数から f_merge を推定 - 経験則 **f0_init = 1.5 × f_merge** で RD フィットを開始 - RD フィットから f0_fit, Q_fit, M_f を推定 しかし、ここには重大な欠陥がある。 ### **位相だけでは質量スケールが決まらない** 重力波のインスパイラル位相はスケール不変であり、 $$ h(t; M_c) = h(\lambda t; \lambda M_c) $$ という対称性を持つ。 つまり： - M_c = 10 Msun - M_c = 30 Msun - M_c = 100 Msun **全部同じ形の波形になる。** そのため v1.2 の Mc_fit は **物理単位ではなく、ただの“形のパラメータ”**に過ぎなかった。 結果として： - f0_init が誤ったスケールに固定される - RD フィットも誤った周波数帯を見に行く - M_f が 0.1 Msun など物理的に破綻した値になる という問題が起きていた。 --- ## 2. v1.3 の発想：物理レイヤーを追加する v1.3 では、v1.2 のパイプラインはそのまま残しつつ、 **f0_init だけを“物理レイヤー”で置き換える**という最小限の変更を行った。 物理レイヤーは以下を計算する： 1. **Mc_phys（物理 chirp mass）を手入力** 2. 対称質量比 η を計算 3. 合体前の全質量 M_tot を推定 4. 放射エネルギー率 ε から最終質量 M_f を推定 5. q から最終スピン a_f を推定 6. Berti+ の QNM フィット式から **f220_pred** を計算 7. f0_init = f220_pred として RD フィットに渡す これにより、解析器は初めて **“物理的な質量スケール”** を持つようになった。 --- ## 3. v1.3 の決定的な成果 ### ● GW150914 Mc_phys = 28 Msun を入力すると： - f220_pred ≈ 302 Hz - f0_fit ≈ 303 Hz - Q_fit ≈ 3.2 - M_f ≈ 62 Msun **LIGO 公開値とほぼ一致。** ### ● GW170608 Mc_phys = 8 Msun を入力すると： - f220_pred ≈ 1019 Hz - f0_fit ≈ 995 Hz - Q_fit ≈ 3.0 - M_f ≈ 20 Msun **軽い系の高周波リングダウンが正しく再現された。** v1.2 では絶対に出せなかった結果だ。 #### GW150914のリングダウン分析結果  #### GW170608のリングダウン分析結果  --- ## 4. v1.3 の自動切替：Mc_phys が無効なら v1.2 に戻る v1.3 では、以下の仕様を追加した： ``` if Mc_phys > 0: f0_init = f220_pred # v1.3 else: f0_init = 1.5 * f_merge # v1.2 fallback ``` これにより： - Mc_phys を与えれば **物理スケールの宇宙（v1.3）** - Mc_phys が無効なら **スケール不定の宇宙（v1.2）** という二形態を自動で切り替えられる。 --- ## 5. コード（v1.3 自動切替版） ```python ############################################################### # Ringdown Analyzer v1.3（Mc_phys 自動切替版） # - v1.2 パイプラインはそのまま # - f0_init だけ物理レイヤーで置換 # - Mc_phys が正の実数でなければ v1.2 にフォールバック ############################################################### import h5py import numpy as np from scipy.signal import hilbert, butter, filtfilt from scipy.optimize import curve_fit import matplotlib.pyplot as plt import re ############################################### # 0. データ読み込み ############################################### fname = "H-H1_LOSC_4_V2-1126259446-32.hdf5" # ← GW150914 # fname = "L-L1_GWOSC_O2_4KHZ_R1-1180921856-4096.hdf5" # ← GW170608 m = re.search(r"(\d+)KHZ", fname.upper()) if m: fs = int(m.group(1)) * 1024 else: fs = 4096.0 with h5py.File(fname, "r") as f: strain = f["strain"]["Strain"][()] t0 = f["strain"]["Strain"].attrs["Xstart"] N = len(strain) t = t0 + np.arange(N) / fs ############################################### # 0.5 イベント時刻（手入力 or 自動） ############################################### t_event_input = 1126259462.4 # GW150914 # t_event_input = 1180922494.5 # GW170608 if np.isfinite(t_event_input): t_event = t_event_input idx_event = int((t_event - t0) * fs) else: idx_event = np.argmax(np.abs(strain)) t_event = t[idx_event] ############################################### # 1. インスパイラル区間 ############################################### insp_end = idx_event insp_window = int(0.050 * fs) insp_start = max(0, insp_end - insp_window) t_insp = t[insp_start:insp_end] h_insp = strain[insp_start:insp_end] ############################################### # 2. Hilbert → 位相・瞬時周波数 ############################################### analytic = hilbert(h_insp) phase = np.unwrap(np.angle(analytic)) inst_freq = np.gradient(phase, t_insp) / (2*np.pi) ############################################### # 3. f_merge（最後の1ms） ############################################### merge_window = int(0.001 * fs) merge_window = max(merge_window, 5) inst_tail = inst_freq[-merge_window:] kernel = np.ones(5) / 5.0 inst_tail_smooth = np.convolve(inst_tail, kernel, mode="valid") f_merge = np.max(inst_tail_smooth) * 10 ############################################### # 4. q_insp（v1.2） ############################################### def chirp_mass_model(t, Mc, phi0): return phi0 + (t - t[-1]) * (Mc**(-5/3)) params, _ = curve_fit( chirp_mass_model, t_insp, phase, p0=[20, 0], maxfev=20000 ) Mc_fit = params[0] q_insp = np.clip(1.0 - 0.5 * (Mc_fit / max(Mc_fit, 1)), 0.3, 1.0) ############################################### # 4.5 物理レイヤー（v1.3） ############################################### def physical_layer_v13(Mc_phys, q, f_merge, eps=0.04): eta = q / (1 + q)**2 M_tot = Mc_phys / (eta**(3/5)) M_f = (1 - eps) * M_tot a_f = 0.7 - 0.3*(1 - q) M_geom = M_f * 4.92549095e-6 # [s] f220_pred = (1 - 0.63*(1 - a_f)**0.3) / (2*np.pi*M_geom) print("===== v1.3 Physical Layer =====") print(f"eta = {eta:.5f}") print(f"M_tot = {M_tot:.3f} Msun") print(f"M_f = {M_f:.3f} Msun") print(f"a_f(q) = {a_f:.4f}") print(f"f220_pred = {f220_pred:.2f} Hz") print("================================") return f220_pred ############################################### # 5. f0_init（v1.3 or v1.2 自動切替） ############################################### # ★ Mc_phys を手入力（正の実数なら v1.3、そうでなければ v1.2） Mc_phys = 28.0 # GW150914 # Mc_phys = 8.0 # GW170608 # Mc_phys = -1 # ← 無効値にすると v1.2 に戻る if (Mc_phys is not None) and np.isfinite(Mc_phys) and (Mc_phys > 0): f0_init = physical_layer_v13(Mc_phys, q_insp, f_merge) mode = "v1.3 (physical)" else: f0_init = 1.5 * f_merge mode = "v1.2 (fallback)" print(f"===== f0_init mode: {mode} =====") print(f"f0_init used = {f0_init:.2f} Hz") print("================================") bp_center = f0_init Q_init = 3.0 tau_init = Q_init / (np.pi * f0_init) dt_start = tau_init * Q_init rd_start = t[insp_end] + dt_start window_ms = 12.0 / q_insp window_s = window_ms * 1e-3 rd_end_time = rd_start + window_s mask_rd = (t >= rd_start) & (t <= rd_end_time) t_rd = t[mask_rd] h_rd_raw = strain[mask_rd] t_rd0 = t_rd - rd_start ############################################### # 6. bandpass（v1.2 のまま） ############################################### def bandpass(data, fs, f1, f2, order=4): nyq = fs / 2.0 b, a = butter(order, [f1/nyq, f2/nyq], btype='band') return filtfilt(b, a, data) f1 = f0_init * 0.9 f2 = f0_init * 1.1 h_rd = bandpass(h_rd_raw, fs, f1, f2) ############################################### # 7. RD モデル（prior） ############################################### def rd_model_with_prior(t, A, f0, tau, phi, lam): h = A * np.exp(-t/tau) * np.cos(2*np.pi*f0*t + phi) Q = np.pi * f0 * tau prior = A * lam * (Q - 3.0) h2 = h.copy() h2[-1] += prior return h2 ############################################### # 8. curve_fit ############################################### lam = q_insp p0_rd = [h_rd[0], f0_init, tau_init, 0.0] params_rd, _ = curve_fit( lambda tt, A, f0, tau, phi: rd_model_with_prior(tt, A, f0, tau, phi, lam), t_rd0, h_rd, p0=p0_rd, maxfev=20000 ) A_fit, f0_fit, tau_fit, phi_fit = params_rd Q_fit = np.pi * f0_fit * tau_fit ############################################### # 9. PRINT ############################################### print("===== RD Fit Result =====") print(f"A_fit = {A_fit:.3e}") print(f"f0_fit = {f0_fit:.2f} Hz") print(f"tau_fit = {tau_fit*1000:.3f} ms") print(f"Q_fit = {Q_fit:.2f}") print(f"phi_fit = {phi_fit:.3f} rad") print("RD bandpass = [{:.1f}, {:.1f}] Hz".format(f1, f2)) print("====================================================") ############################################### # 10. Remnant BH ############################################### a_f = 1.0 - 1.0/(2.0 * Q_fit) M_geom = (1 - 0.63*(1 - a_f)**0.3) / (2*np.pi*f0_fit) M_solar = M_geom / 4.92549095e-6 print("===== Remnant BH (from RD fit) =====") print(f"Remnant spin a_f = {a_f:.4f}") print(f"Remnant mass M_f = {M_solar:.2f} Msun") print("====================================") ############################################### # 11. 図 ############################################### h_model = A_fit * np.exp(-t_rd0/tau_fit) * np.cos(2*np.pi*f0_fit*t_rd0 + phi_fit) plt.figure(figsize=(10,5)) plt.plot(t_rd0*1000, h_rd, label="RD data (bandpassed)", lw=1.5) plt.plot(t_rd0*1000, h_model, label="RD fit", lw=2.0) plt.xlabel("Time since rd_start [ms]") plt.ylabel("Strain") plt.title(f"Ringdown: data vs fit ({mode})") plt.legend() plt.grid(True) plt.tight_layout() plt.show() ``` --- ## 6. まとめ：解析器が“宇宙のものさし”を取り戻した v1.2 の世界では： - 位相だけ - スケール不定 - f0_init が経験則 - RD フィットが破綻しやすい という限界があった。 v1.3 では： - Mc_phys を与えるだけで - f220_pred が正しいスケールに復元され - RD フィットも正しい周波数帯にロックし - M_f, a_f が物理的に整合する という **本物の IMR consistency test に近い挙動**が実現した。 --- ## おわりに この v1.3 の進化は、 **「軽量解析器でも、物理スケールを正しく扱える」** という強いメッセージになる。 次のステップは： - Mc_phys の自動推定 - f220_pred と f0_fit の比較プロット - IMR consistency の可視化 などが考えられる。 --- 次のバージョン： [Ringdown Toybox v1.5 外生パラメータなしで f220 を推定するリングダウン解析器](https://talkwithgai.blogspot.com/2026/06/blog-post_24.html) 英語版： [Ringdown Analyzer v1.2 → v1.3 How the lightweight ringdown analyzer finally recovered the “scale of the universe”](https://talkwithgai.blogspot.com/2026/06/blog-post_714.html)