はじめに

この記事では、SS 433のジェットにおけるドップラーシフトの解析方法について説明します。SS 433は連星系を持つ天体であり、非常に高速で噴出するジェットを持つマイクロクエーサーとして知られています。この連星系の中心にはブラックホールまたは中性子星が存在すると考えられていますが、現時点ではどちらであるかは明確に特定されていません。ジェットは光速の約26%の速度で噴出されており、この高速運動が観測されるドップラーシフトに大きな影響を与えます。

ジェットの動きに伴うドップラーシフトを解析することで、ジェットの速度、方向、および周期性を理解できます。特に、相対論的なジェットを持つマイクロクエーサーやガンマ線バーストでは、非相対論的な天体（例: 超新星残骸）とは異なり、観測者の方向の速度と角度が相対論的な効果において重要な要素となるため、注意が必要です。

SS 433の周期性

SS 433には以下のような主要な周期性があります：

軌道周期 (Orbital Period): 連星系において、2つの星が共通の重心を周回する周期です。SS 433では、この周期は約13.1日とされています。この軌道周期は、特に可視光とX線観測において、食に伴う増減光として観測されます¹。
自転周期 (Spin Period): 主星（ブラックホールあるいは中性子星）の自転に対応する周期です。ブラックホールの場合、一般的にこの周期を直接測定することは非常に困難です。
歳差周期 (Precession Period): ジェットの自転軸がゆっくりと円を描くように変化する周期で、SS 433では約162日とされています²。
章動周期 (Nutation Period): 自転軸の歳差運動に加えて、軸がさらに小さな円を描くように揺れ動く周期で、約6.29日とされています。これは歳差運動に加えて追加的な変動を引き起こします³。

ドップラーシフトのモデル

この記事では、SS 433のジェットのドップラーシフトを解析するために、以下の2つのモデルを使用します：

Six-Parameterモデル: 主に歳差運動に関連するモデルで、ジェットの基本的な動きを記述します。このモデルでは、ジェットの速度や傾斜角度など、6つの主要なパラメータが考慮されます。
Nodding-Motionモデル: このモデルは、Six-Parameterモデルにさらに章動や追加的な揺動を考慮したものです。これにより、より複雑なドップラーシフトの変動を説明することができます。

これらのモデルは、Katz et al. (1982)およびGies et al. (2002)などの研究で使用されています。この記事では、Katz et al. (1982)のモデルに基づき、$\dot{\Omega}_s=0, \dot{\omega}, \xi^\prime, \theta^\prime=0$ といった近似を用いてドップラーシフトの計算を行います。これにより、SS 433のジェットがどのように観測されるかを解析します。

モデルの数式表現

SS 433のジェットのドップラーシフトを解析するために使用される2つの主要なモデル、Six-ParameterモデルとNodding-Motionモデルについて説明します。

Six-Parameterモデルとは

Six-Parameterモデルは、SS 433のジェットの基本的な力学モデルです。このモデルでは、以下の6つのパラメータに基づいてジェットのドップラーシフトを計算します。

z \equiv \frac{\Delta \lambda}{\lambda} = \gamma \left[ \pm \beta \sin i \sin \theta \cos \left( \xi - \frac{1}{2}\pi \right) \pm \beta \cos i \cos \theta + 1 \right] - 1

この式は、ジェットの相対論的ドップラー効果による赤方偏移または青方偏移を表しています。各項の意味は次のとおりです：

正負の符号: 正の符号（+）は後退ジェット（Receding Jet）を、負の符号（−）は接近ジェット（Approaching Jet）を表します。
$z$: 赤方偏移または青方偏移の量を示し、これは波長の変化率 $\frac{\Delta \lambda}{\lambda}$ として定義されます。
$\gamma$: ローレンツ因子（$\gamma=1/\sqrt{1-\beta^2}$）です。
$\beta$: 光速に対するジェットの速度の割合で表されます。
$i$: ジェット軌道の傾斜角で、観測者から見たジェットの軌道面の傾きを示します。
$\theta$: ジェットが進む方向と軌道面との間の角度です。
$\xi$: ジェットの歳差運動の位相です。加速度項を無視すると、$$\xi = \Omega_s(0)t + \xi_0$$ と表されます。
- $t$: 観測時刻です。
- $\xi_0$: 基準時刻の位相です。
- $\Omega_s(0)$: 歳差運動の角速度で、時間とともに歳差運動の位相 $\xi$ が変化します。
- $P_s(0)$: 歳差運動の周期で、歳差運動の角速度から計算されます。$$P_s(0) \equiv \frac{2\pi}{\Omega_s(0)}$$ と定義されます。

Six-Parameterモデルは、ジェットの歳差運動の影響を考慮し、ドップラーシフトの基本的な変動を計算しますが、章動などジェット軸自体が周期的に揺れるNoddingの効果は含まれていません。

Nodding-Motionモデルとは

Nodding-Motionモデルは、Six-Parameterモデルに加え、Noddingの影響を考慮したモデルです。Katz et al. 1982によると、観測的に6.3日の章動周期の他に、5.8日周期の微弱な変動が確認されています。数式で次のように表されます。

\delta z = \pm A_{5.8} \cos \left[ \omega_{5.8}(t - t_0^{5.8}) - 3\xi_0 + \frac{\pi}{2} \right] 
         \pm A_{6.3} \cos \left[ \omega_{6.3}(t - t_0^{6.3}) - \xi_0 - \frac{\pi}{2} \right]

正負の符号: 正の符号（+）は後退ジェット（Receding Jet）を、負の符号（−）は接近ジェット（Approaching Jet）を表します。
下付き上付きの添え字の$5.8$、$6.3$: それぞれ5.8日、6.3日周期の変動成分のパラメータに対応します。
$\delta z$: Six-Parameterモデルによる基本的なドップラーシフト $z$ に対する変動量です。
$A$: 振幅を示します。
$\omega$: 角速度を示します。
$t_0$: 基準時刻を表します。
$\xi_0$, $t$: Six-Parameterモデルと共通のパラメータを表します。

このモデルにより、歳差運動に加えた章動の影響を考慮した、より精密なドップラーシフトの解析が可能となります。

実装

実装コードは

からもご覧いただけます。

SS 433ドップラーシフト解析

import numpy as np
import matplotlib.pyplot as plt
from astropy.time import Time
from datetime import datetime, timedelta


def calc_ss433_doppler_shift(t_jd, nodding_flag=True,
                             beta=0.2591, theta_deg=19.81, 
                             i_deg=78.88, xi0_deg=-35.2, ps=168.72, 
                             a_58=3.53e-3, t0_58=2443551.20, p_58=5.841, 
                             a_63=-6.82e-3, t0_63=2443551.58, p_63=6.3057):
    """
    SS 433のドップラーシフトを計算する関数。
    デフォルト値: Katz et al. 1982
    https://ui.adsabs.harvard.edu/abs/1982ApJ...260..780K/abstract
    
    Parameters:
    - t_jd: JD単位の計算する時刻
    - nodding_flag: ノッディングの影響を考慮するかどうか
    - beta: 相対速度 (光速の単位)
    - theta_deg: 放射軸の傾き角 (度)
    - i_deg: 系の傾斜角 (度)
    - xi0_deg: 初期位相 (度)
    - ps: 歳差周期 (日)
    - a_58, t0_58, p_58: ノッディングの5.8日成分のパラメータ
    - a_63, t0_63, p_63: ノッディングの6.3日成分のパラメータ
    
    Returns:
    - z_r: 後退ジェットのドップラーシフト
    - z_b: 接近ジェットのドップラーシフト
    """
    
    # 角度をラジアンに変換
    theta = np.radians(theta_deg)
    i = np.radians(i_deg)
    xi0 = np.radians(xi0_deg)

    # ローレンツ因子の計算
    gamma = 1 / np.sqrt(1 - beta**2) 
    
    # 位相の計算
    xi = 2 * np.pi * (t_jd / ps) + xi0
    
    # レッドシフトとブルーシフトの計算
    z_r = gamma * (beta * np.sin(i) * np.sin(theta) * np.cos(xi - 0.5 * np.pi) + beta * np.cos(i) * np.cos(theta) + 1) - 1
    z_b = gamma * (-beta * np.sin(i) * np.sin(theta) * np.cos(xi - 0.5 * np.pi) - beta * np.cos(i) * np.cos(theta) + 1) - 1

    if nodding_flag:
        # ノッディングの項を計算
        def calculate_nodding(t, a, t0, p, xi0_shift):
            omega = 2 * np.pi / p
            return a * np.cos(omega * (t - t0) + xi0_shift)

        delta_z_58 = calculate_nodding(t_jd, a_58, t0_58, p_58, -3 * xi0 + np.pi / 2)
        delta_z_63 = calculate_nodding(t_jd, a_63, t0_63, p_63, -xi0 - np.pi / 2)

        delta_z = delta_z_58 + delta_z_63

        z_r += delta_z
        z_b -= delta_z

    return z_r, z_b


def convert_z_ene_funcs(ene_rest):
    def z_to_ene(z):
        ene_obs = ene_rest / (1 + z)
        return ene_obs

    def ene_to_z(ene_obs):
        z = (ene_rest / (ene_obs+1e-15)) - 1
        return z
    
    return z_to_ene, ene_to_z

def main():
    # 期間の設定
    start_date = datetime(2024, 1, 1)
    end_date = datetime(2025, 1, 1)
    delta = timedelta(days=0.1)
    num_days = (end_date - start_date).days / 0.1

    # 計算に用いる時刻配列の準備
    dates = np.array([start_date + timedelta(days=i * 0.1) for i in range(int(num_days) + 1)])
    jd_ts = Time(dates).jd

    # zの計算
    z_rs, z_bs = calc_ss433_doppler_shift(jd_ts, nodding_flag=False)
    z_r_nods, z_b_nods = calc_ss433_doppler_shift(jd_ts, nodding_flag=True)

    # プロットの作成
    fig, ax = plt.subplots(figsize=(12, 6))

    ax.plot(dates, z_rs, label=f'z_r six-parameter model', c='r', ls='--')
    ax.plot(dates, z_bs, label=f'z_b six-parameter model', c='b', ls='--')
    ax.plot(dates, z_r_nods, label=f'z_r nodding model', c='r')
    ax.plot(dates, z_b_nods, label=f'z_b nodding model', c='b')
    ax.set_xlabel('Date')
    ax.set_ylabel('Doppler Shift')
    ax.grid(True)

    # 特定の期間に対してグレーのマスクを追加
    mask_start_date = datetime(2024, 8, 10, 1, 1, 1)  # 観測開始時刻
    mask_end_date = datetime(2024, 8, 15, 1, 1, 1)  # 観測終了時刻
    ax.axvspan(mask_start_date, mask_end_date, label='observation period', color='gray', alpha=0.7)

    # エネルギーの二次軸(Fe XXV Kαの例) 
    fe_xxv_kalpha_ene = 6700
    z_to_ene, ene_to_z = convert_z_ene_funcs(fe_xxv_kalpha_ene)
    secax_y = ax.secondary_yaxis('right', functions=(z_to_ene, ene_to_z))
    secax_y.set_ylabel('Fe XXV kα Energy (eV)')

    # # JD時刻の二次軸
    jd_offset = 58484 + 2400000.5
    secax_x = ax.secondary_xaxis('top')
    secax_x.set_xlabel(f'JD - {jd_offset}')
    jd_labels = np.linspace(min(jd_ts), max(jd_ts), len(ax.get_xticks()))
    secax_x.set_xticks(ax.get_xticks())
    secax_x.set_xticklabels([f'{jd - jd_offset:.0f}' for jd in jd_labels])

    plt.legend(loc='upper left', bbox_to_anchor=(1.07, 1))
    fig.tight_layout()
    plt.savefig('ss433_phase.png')
    plt.show()

if __name__ == '__main__':
    main()

出力結果

JD時刻に基づいて上図はSix-Parameterモデル（破線）とNodding-Motionモデル（実線）の結果を表示しています。
SS 433の2つのジェット（z_rが後退ジェット、z_bが接近ジェット）のドップラーシフトとそれに対応するFe XXV Kα線のエネルギーを2軸で示しています。
Nodding-Motionモデルでは軽微な周期的な変動も見えているのがわかります。

補足

パラメータについて

この記事で使用しているパラメータKatz et al. (1982)およびGies et al. (2002)に基づいていますが、Six-ParameterモデルにおいてGies et al. (2002)のパラメータを使用することが多いです。

両者のパラメータの大きな違いは歳差周期であり、Katz et al. (1982)では約168.72日としていますが、Gies et. al (2002)では162.15日と推定しています。また、この軌道周期をサポートする解析結果として、Eikenberry et al. (2001)による162.375日やDavydov et al. (2008)の162.250日などがあります。

ただし、モデルの選択は解析結果の精度に影響を与えるため、どのモデルを使用したのかを十分に考慮する必要があります。例えば、Eikenberry et al. (2001)は、約20年にわたる長期観測データに基づき、スペクトルを平滑化することで章動の影響を抑えたと説明し、Six-Parameterモデルを採用しています。一方、Davydov et al. (2008)は、章動の影響を考慮した別のモデルを使用しているため、両者の間でパラメータに若干の違いが生じています。この点については、Davydov et al. (2008)で言及されています。

さまざまな輝線の同時プロット

SS 433では、Fe He-like Kα（Fe XXV）やH-like Lyα（Fe XXVI）、Ni He-like Kα（Ni XXVII）、Ni H-like Lyα（Ni XXVIII）など、さまざまな輝線が観測されています⁴。ここでは、これらの複数の輝線を表示するコードを示します。

なお、輝線のエネルギーは検出器の分解能により広がりを持つため、複数のスピン状態や量子状態から発生するさまざまな輝線が観測されます。この点については、こちらを参照してください。

また、輝線のエネルギーを調べる際には、NISTのAtomic Spectra Databaseが便利です。"Lines"を選択し、ページ遷移後に"Search for Photon Energy"を選びます。続いて、"Spectrum"欄に調べたい輝線（例: Fe XXV）を入力することで、対応するエネルギーを確認できます。本記事では、このサイトの対応する遷移の中で最も強度の高い強度のエネルギーを使用しています。

複数の輝線のドップラーシフトの描画

def main():
    # 期間の設定
    start_date = datetime(2024, 1, 1)
    end_date = datetime(2025, 1, 1)
    delta = timedelta(days=0.1)
    num_days = (end_date - start_date).days / 0.1

    # 計算に用いる時刻配列の準備
    dates = np.array([start_date + timedelta(days=i * 0.1) for i in range(int(num_days) + 1)])
    jd_ts = Time(dates).jd

    # zの計算
    z_rs, z_bs = calc_ss433_doppler_shift(jd_ts, nodding_flag=False)
    z_r_nods, z_b_nods = calc_ss433_doppler_shift(jd_ts, nodding_flag=True)

    # 輝線のエネルギー
    fe_xxv_kalpha_ene = 6700
    z_to_ene, _ = convert_z_ene_funcs(fe_xxv_kalpha_ene)
    fe_xxv_kalpha_ene_rs = z_to_ene(z_rs)
    fe_xxv_kalpha_ene_bs = z_to_ene(z_bs)
    fe_xxv_kalpha_ene_rs_nods = z_to_ene(z_r_nods)
    fe_xxv_kalpha_ene_bs_nods = z_to_ene(z_b_nods)

    fe_xxvi_lyalpha_ene = 6952
    z_to_ene, _ = convert_z_ene_funcs(fe_xxvi_lyalpha_ene)
    fe_xxvi_lyalpha_ene_rs = z_to_ene(z_rs)
    fe_xxvi_lyalpha_ene_bs = z_to_ene(z_bs)
    fe_xxvi_lyalpha_ene_rs_nods = z_to_ene(z_r_nods)
    fe_xxvi_lyalpha_ene_bs_nods = z_to_ene(z_b_nods)

    ni_xxvii_kalpha_ene = 7806
    z_to_ene, _ = convert_z_ene_funcs(ni_xxvii_kalpha_ene)
    ni_xxvii_kalpha_ene_rs = z_to_ene(z_rs)
    ni_xxvii_kalpha_ene_bs = z_to_ene(z_bs)
    ni_xxvii_kalpha_ene_rs_nods = z_to_ene(z_r_nods)
    ni_xxvii_kalpha_ene_bs_nods = z_to_ene(z_b_nods)

    ni_xxviii_lyalpha_ene = 8073
    z_to_ene, _ = convert_z_ene_funcs(ni_xxviii_lyalpha_ene)
    ni_xxviii_lyalpha_ene_rs = z_to_ene(z_rs)
    ni_xxviii_lyalpha_ene_bs = z_to_ene(z_bs)
    ni_xxviii_lyalpha_ene_rs_nods = z_to_ene(z_r_nods)
    ni_xxviii_lyalpha_ene_bs_nods = z_to_ene(z_b_nods)

    # プロットの作成
    fig, ax = plt.subplots(figsize=(12, 6))

    # Feのプロット
    ax.plot(dates, fe_xxv_kalpha_ene_rs, label=f'Fe XXV Kα red jet six-parameter model', c='r', ls='--')
    ax.plot(dates, fe_xxv_kalpha_ene_bs, label=f'Fe XXV Kα blue jet six-parameter model', c='b', ls='--')
    ax.plot(dates, fe_xxv_kalpha_ene_rs_nods, label=f'Fe XXV Kα red jet nodding model', c='r')
    ax.plot(dates, fe_xxv_kalpha_ene_bs_nods, label=f'Fe XXV Kα blue jet nodding model', c='b')

    ax.plot(dates, fe_xxvi_lyalpha_ene_rs, label=f'Fe XXVI Lyα red jet six-parameter model', c='r', ls='--', alpha=0.4)
    ax.plot(dates, fe_xxvi_lyalpha_ene_bs, label=f'Fe XXVI Lyα blue jet six-parameter model', c='b', ls='--', alpha=0.4)
    ax.plot(dates, fe_xxvi_lyalpha_ene_rs_nods, label=f'Fe XXVI Lyα red jet nodding model', c='r', alpha=0.4)
    ax.plot(dates, fe_xxvi_lyalpha_ene_bs_nods, label=f'Fe XXVI Lyα blue jet nodding model', c='b', alpha=0.4)

    # Niのプロット
    ax.plot(dates, ni_xxvii_kalpha_ene_rs, label=f'Ni XXVII Kα red jet six-parameter model', c='m', ls='--', alpha=0.3)
    ax.plot(dates, ni_xxvii_kalpha_ene_bs, label=f'Ni XXVII Kα blue jet six-parameter model', c='c', ls='--', alpha=0.3)
    ax.plot(dates, ni_xxvii_kalpha_ene_rs_nods, label=f'Ni XXVII Kα red jet nodding model', c='m', alpha=0.3)
    ax.plot(dates, ni_xxvii_kalpha_ene_bs_nods, label=f'Ni XXVII Kα blue jet nodding model', c='c', alpha=0.3)

    ax.plot(dates, ni_xxviii_lyalpha_ene_rs, label=f'Ni XXVIII Lyα red jet six-parameter model', c='m', ls='--', alpha=0.12)
    ax.plot(dates, ni_xxviii_lyalpha_ene_bs, label=f'Ni XXVIII Lyα blue jet six-parameter model', c='c', ls='--', alpha=0.12)
    ax.plot(dates, ni_xxviii_lyalpha_ene_rs_nods, label=f'Ni XXVIII Lyα red jet nodding model', c='m', alpha=0.12)
    ax.plot(dates, ni_xxviii_lyalpha_ene_bs_nods, label=f'Ni XXVIII Lyα blue jet nodding model', c='c', alpha=0.12)

    ax.set_xlabel('Date')
    ax.set_ylabel('Energy (eV)')
    ax.grid(True)

    # 特定の期間に対してグレーのマスクを追加
    mask_start_date = datetime(2024, 8, 10, 1, 1, 1)  # 観測開始時刻
    mask_end_date = datetime(2024, 8, 15, 1, 1, 1)  # 観測終了時刻
    ax.axvspan(mask_start_date, mask_end_date, label='observation period', color='gray', alpha=0.7)
    ax.legend()

    # # jdの二次軸
    jd_offset = 58484 + 2400000.5
    secax_x = ax.secondary_xaxis('top')
    secax_x.set_xlabel(f'JD - {jd_offset}')
    jd_labels = np.linspace(min(jd_ts), max(jd_ts), len(ax.get_xticks()))
    secax_x.set_xticks(ax.get_xticks())
    secax_x.set_xticklabels([f'{jd - jd_offset:.0f}' for jd in jd_labels])

    plt.legend(loc='upper left', bbox_to_anchor=(1.02, 1))
    fig.tight_layout()
    plt.savefig('ss433_phase_Fe_Ni.png')
    plt.show()

if __name__ == '__main__':
    main()

出力結果

上図では、複数の輝線が交わる時期があることが確認できます。このように、SS 433のスペクトル解析は非常に複雑であり、解析が難しいことがわかります。図中では、各輝線を任意の透明度でプロットしていますが、輝線の強度に応じて厳密に明度や透明度を調整することで、解析がより視覚的にわかりやすくなると思います。

まとめ

この記事では、SS 433のジェットにおけるドップラーシフトの解析方法について解説しました。SS 433は非常に高速で噴出するジェットを持つ連星系で、そのジェットのドップラーシフトを解析することで、ジェットの速度、方向、および運動に関連する周期性を理解することができます。

主に使用したモデルは、Six-ParameterモデルとNodding-Motionモデルであり、これらを用いてSS 433のドップラーシフトを数値的に解析しました。さらに、Fe XXV Kα線や他の輝線のエネルギーシフトをプロットし、ジェットの挙動に伴うエネルギー変動を視覚化しました。

また、NISTのAtomic Spectra Databaseを利用して輝線のエネルギーを調べる方法についても説明し、最も強度の高い輝線のエネルギーを用いて解析を行いました。本記事がSS 433のドップラーシフトの解析に役立つことを願っています。

軌道周期とライトカーブの関係: e.g., Cherepashchuk et al. (2009), Goranskij et al. (2011), Burenin et al. (2011) ↩
歳差周期の推定値: e.g., Margon et al. (1979): 164日, Katz et al. (1982): 168.72日, Eikenberry et al. (2001): 162.375日, Gies et al. (2002): 162.15日, Davydov et al. (2008): 162.250日 ↩
章動周期の観測事例: e.g., Vermeulen et al. (1993), Davydov et al. (2008) ↩
さまざまなX線観測衛星で観測された輝線: e.g., ASCA: Kotani et al. (1994) XMM-Newton: Medvedev et al. (2018), Suzaku: Kubota et al. (2010), Chandra: Medvedev et al. (2019) ↩