はじめに

前回の記事では、SS 433 の歳差ジェットを三次元空間の速度ベクトルとして表しました。

SS 433 のジェットは、光速の約 $0.26$ 倍という高速で飛び出しており、さらにジェット軸が歳差運動しています。
そのため、ジェットの向きは時間とともに変わり、三次元空間ではらせん状のような構造を描きます。

ただ、私たちが実際に観測するのは、三次元空間そのものではありません。

観測画像で見えるのは、天球面に投影された二次元の構造です。
つまり、三次元的なジェット方向を計算しただけでは、まだ「観測画像の上でどう見えるか」までは見えてきません。

そこで今回は、前回の歳差ジェットモデルをもとにして、ジェットの軌道を空の上に投影してみます。

ここで特に重要になるのが、光の到達時間効果です。

SS 433 のジェットは相対論的な速度で運動しているため、単純に

位置 = 速度 \times 時間

として空に投影するだけでは、少し足りません。

ジェットがこちらに向かって動く成分を持つ場合、後から放たれた光は、より観測者に近い位置から出発します。
そのため、先に放たれた光との差が縮まり、見かけ上、ジェットが大きく伸びて見えます。

この記事では、この光の到達時間効果を入れて、SS 433 の見かけのジェット軌道を Python で可視化してみます。

実装コードは Google Colab からもお試しいただけます：
https://colab.research.google.com/drive/1nOjwGxrOUxVwFFFVKJse4N5k6tyOf63-?usp=sharing

観測されるのは「今のジェット」だけではない

まず大事なのは、観測画像に写っているジェットの各点が、すべて同じ時刻に放出されたものではないということです。

中心に近い部分は、比較的最近放出されたジェットです。
一方、中心から離れた部分は、もっと昔に放出されたジェットです。

観測時刻を

t_{\rm obs}

ジェット要素の見かけ上の年齢を

t_{\rm age}

とします。

ここでの $t_{\rm age}$ は、観測時刻から見て「どれくらい前に放出されたジェット要素として見えているか」を表す量です。

このとき、そのジェット要素が放出された時刻は

t_{\rm emit} = t_{\rm obs} - t_{\rm age}

です。

SS 433 のジェットは歳差運動しているため、ジェットの向きは現在の観測時刻 $t_{\rm obs}$ だけで決まるわけではありません。
それぞれのジェット要素が放出された時刻 $t_{\rm emit}$ における歳差位相で決まります。

そのため、空の上に見えるジェットは単純な直線にはならず、過去の歳差運動の履歴を反映した曲がった軌跡になります。

今回のコードでは、この考え方に基づいて、各ジェット要素ごとに放出時刻を計算し、その時刻に対応する歳差位相を求めています。

三次元速度を空の上に投影する

前回の記事に倣い、ジェット速度を光速で割ったものを

\boldsymbol{\beta}
=
(\beta_x, \beta_y, \beta_z)

と書くことにします。

ここでは、

$x$ 方向：視線方向
$y$ 方向：Dec 方向
$z$ 方向：RA 方向

として考えます。

この座標系では、$\beta_x$ が視線方向の速度成分、$\beta_y$ と $\beta_z$ が天球面上の速度成分です。

SS 433 のジェットは歳差運動しているため、速度ベクトルの向きは時間とともに変わります。
歳差位相を $\psi(t)$ とすると、観測者方向を基準にした速度成分は次のように書けます。

\left\{
\begin{aligned}
\beta_x(t) &= \pm\beta
\left[
\cos i\cos\theta
+
\sin i\sin\theta\cos\psi(t)
\right],
\\[4pt]
\beta_y(t) &= \pm\beta
\left[
\cos\chi
\left(
\sin i\cos\theta
-
\cos i\sin\theta\cos\psi(t)
\right)
+
s_{\rm rot}\sin\chi\,\sin\theta\sin\psi(t)
\right],
\\[4pt]
\beta_z(t) &= \pm\beta
\left[
\sin\chi
\left(
\sin i\cos\theta
-
\cos i\sin\theta\cos\psi(t)
\right)
-
s_{\rm rot}\cos\chi\,\sin\theta\sin\psi(t)
\right].
\end{aligned}
\right.

ここで、

$\beta = v/c$：ジェット速度を光速で割った値
$i$：歳差軸の傾斜角
$\theta$：歳差円錐の半開口角
$\chi$：歳差軸の位置角
$\psi(t)$：歳差位相
$s_{\rm rot}$：歳差回転の向きを表す符号

です。

また、$\pm$ は双方向ジェットを表しています。
ここでは便宜的に、$+$ を approaching-side jet、$-$ を receding-side jet として描いています。

観測画像上の位置に直接関係するのは、天球面上の成分である

\beta_y(t),\quad \beta_z(t)

です。

もし光の到達時間効果を考えなければ、ジェット要素の天球面上の位置は、おおまかに

y \sim \beta_y(t)c t,
\quad
z \sim \beta_z(t)c t

となります。

つまり、天球面上の速度成分に時間をかければ、空の上での位置が決まる、という考え方です。

ただし、SS 433 のようにジェットが高速で動いている場合には、視線方向の運動も見かけの位置に効いてきます。

その視線方向の運動を表しているのが

\beta_x(t)

です。

次に見るように、この $\beta_x(t)$ が、光の到達時間効果を通して見かけの位置を変えます。

光の到達時間効果

ここからが今回の記事の中心です。

天球面上の位置だけを考えるなら、$\beta_y(t)$ と $\beta_z(t)$ に時間をかければよさそうに見えます。
しかし実際には、視線方向速度 $\beta_x(t)$ によって、光が観測者に届くタイミングが変わります。

その効果を入れると、見かけの位置は次のように書けます。

y_{\rm app}
=
\frac{\beta_y(t_{\rm emit}) c t_{\rm age}}{1-\beta_x(t_{\rm emit})},
\quad
z_{\rm app}
=
\frac{\beta_z(t_{\rm emit}) c t_{\rm age}}{1-\beta_x(t_{\rm emit})}.

ここで重要なのが、分母の

1-\beta_x(t_{\rm emit})

です。

もし

\beta_x > 0

なら、ジェットは観測者に近づく成分を持ちます。

このとき、

1-\beta_x < 1

となるため、見かけの位置は大きくなります。
つまり、こちらに向かう成分を持つジェットは、空の上で引き伸ばされて見えます。

一方で、

\beta_x < 0

なら、ジェットは観測者から遠ざかる成分を持ちます。

このとき、

1-\beta_x > 1

となるため、見かけの位置は小さくなります。

このように、同じ速さで反対向きに飛んでいる双方向ジェットであっても、観測画像の上では少し非対称に見えることがあります。

このあたりが、相対論的ジェットを眺めるときの面白いところだと思います。

なお、ここでは見かけの位置の式をそのまま使いましたが、見かけ速度の式がどこから出てくるのかについては、こちらの記事でも整理しています。よろしければ合わせてご覧ください。

距離を仮定して角度スケールに変換する

ここまでで、ジェット要素の見かけの位置を

y_{\rm app},
\quad
z_{\rm app}

として求めました。

これは物理的な長さの単位、たとえば cm で表された位置です。
しかし、観測画像で見る量は、RA offset や Dec offset のような角度です。

そのため、SS 433 までの距離を $D$ とすると、天球面上の角度は小角度近似で

\Delta {\rm Dec}
=
\frac{y_{\rm app}}{D},
\quad
\Delta {\rm RA}
=
\frac{z_{\rm app}}{D}

として求められます。

今回のコードでは、SS 433 までの距離として

D = 5.5\,{\rm kpc}

を仮定しています。

したがって、最終的な見かけの位置は

\Delta {\rm Dec}
=
\frac{y_{\rm app}}{D}
\times 206265\ {\rm arcsec},
\quad
\Delta {\rm RA}
=
\frac{z_{\rm app}}{D}
\times 206265\ {\rm arcsec}

として、arcsec 単位に変換しています。

ここで、距離 $D$ を変えると、同じ物理的なジェット構造でも、天球面上での見かけのスケールが変わります。
距離が大きければ、同じ物理サイズの構造は小さな角度に見えます。
逆に距離が小さければ、同じ構造はより大きな角度に見えます。

SS 433 の距離についてはいくつか議論がありますが、Blundell & Bowler (2004) では、歳差ジェットが作るジグザグ、あるいはコルクスクリュー状の構造の対称性と角周期性から、$5.5 \pm 0.2\,{\rm kpc}$ という距離が示されています。

今回の可視化でも、$D=5.5\,{\rm kpc}$ を仮定することで、計算した物理的な軌道を RA/Dec offset の arcsec スケールに変換しています。

つまり、この図の形は歳差運動と光の到達時間効果で決まり、図の角度スケールは仮定した距離にも依存します。

RA/Dec offset と天球座標としての RA/Dec の違い

ここで、図の軸名について少し補足しておきます。

今回の図では、横軸を

RA offset [arcsec]

縦軸を

Dec offset [arcsec]

としています。

ただし、ここでの RA offset / Dec offset は、天球座標としての絶対的な RA, Dec を計算しているという意味ではありません。

あくまで、SS 433 の中心を原点にした局所的な接平面上で、

$z$ 方向の変位を RA 方向の見かけのオフセット
$y$ 方向の変位を Dec 方向の見かけのオフセット

として、arcsec 単位で表示しているものです。

実際の天球座標としての RA, Dec を扱う場合、RA 方向の角距離には赤緯に依存する補正が入ります。
たとえば、同じ RA 差であっても、天球上の実際の角距離は赤緯によって少し変わります。

そのため、本格的に観測画像の WCS と対応させる場合には、RA 方向の $\cos\delta$ 補正や、投影法、WCS 変換を別途扱う必要があります。この辺りの詳細は、こちらの記事を参照してください。

一方、今回の可視化では、SS 433 の近傍の数 arcsec スケールを考えているため、局所的な接平面上の小角度近似として扱っています。

つまり、ここで表示している RA offset / Dec offset は、

厳密な天球座標変換ではなく、モデルから得られた見かけの変位を arcsec 単位で表示したもの

です。

今回の目的は、実データの WCS に直接重ねることではなく、歳差運動・視線方向速度・光の到達時間効果によって、ジェットが空の上でどのように見えるかを直感的に見ることです。

そのため、図は「観測画像に近い向きで表示したモデル軌道」として見るのがよいと思います。

approaching jet と receding jet

SS 433 では、双方向にジェットが出ているため、この記事では

approaching-side jet
receding-side jet

の2つを描いています。

ただし、ここで少し注意があります。

名前としては approaching jet / receding jet と呼んでいますが、歳差運動しているため、軌道上のすべての場所で常に同じ視線方向速度を持つわけではありません。

つまり、approaching-side jet であっても、軌道の一部では視線方向速度が小さくなったり、場合によっては符号が変わったりします。

今回の図では、視線方向速度 $\beta_x$ の符号によって線種を分けています。

実線：こちらに向かう成分を持つ部分
点線：遠ざかる成分を持つ部分

このように分けることで、空の上の軌道だけでなく、その場所の視線方向運動も同時に見ることができます。

単に「青が approaching、赤が receding」と見るだけでなく、軌道上のどの部分がこちら向きで、どの部分が遠ざかる向きなのかを見ると、SS 433 のジェット構造が少し立体的に感じられます。

静止画として見る

まずは、ある1つの歳差位相におけるジェット軌道を描いてみます。

コードはこちら

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator

# ==========================================================
# SS 433 projected jet trajectory
# comparison:
#   upper panel : without light-travel-time effect
#   lower panel : with light-travel-time effect
#
# age_day      : apparent age of each jet element
# obs_time_day : observation time
# t_emit_day   : ejection time
#
# t_emit_day = obs_time_day - age_day
# ==========================================================

# -----------------------------
# physical constants
# -----------------------------
c = 2.99792458e10          # [cm/s]
day = 86400.0              # [s]
kpc = 3.085677581e21       # [cm]
rad_to_arcsec = 206265.0

# -----------------------------
# SS 433 parameters
# -----------------------------
D_kpc = 5.5
D = D_kpc * kpc

beta = 0.26

# geometry
inc_deg = 78.83
theta_deg = 19.85
chi_deg = 98.2

# precession
P_prec_day = 162.375
s_rot = -1

# convert to rad
i = np.deg2rad(inc_deg)
theta = np.deg2rad(theta_deg)
chi = np.deg2rad(chi_deg)


# ==========================================================
# helper functions
# ==========================================================

def precession_phase(age_day, obs_time_day=0.0):
    """
    Precession phase of each jet element.

    A jet element observed with age_day was ejected at

        t_emit_day = obs_time_day - age_day

    Therefore, the jet direction should be evaluated
    at t_emit_day.
    """

    t_emit_day = obs_time_day - age_day
    psi = 2 * np.pi * t_emit_day / P_prec_day

    return psi


def beta_components(age_day, sign=+1, obs_time_day=0.0):
    """
    Velocity components in the observer frame.

    sign = +1 : approaching-side jet
    sign = -1 : receding-side jet

    x : line-of-sight direction
    y : Dec direction
    z : RA direction
    """

    psi = precession_phase(
        age_day=age_day,
        obs_time_day=obs_time_day
    )

    bx = sign * beta * (
        np.cos(i) * np.cos(theta)
        + np.sin(i) * np.sin(theta) * np.cos(psi)
    )

    by = sign * beta * (
        np.cos(chi)
        * (
            np.sin(i) * np.cos(theta)
            - np.cos(i) * np.sin(theta) * np.cos(psi)
        )
        + s_rot * np.sin(chi) * np.sin(theta) * np.sin(psi)
    )

    bz = sign * beta * (
        np.sin(chi)
        * (
            np.sin(i) * np.cos(theta)
            - np.cos(i) * np.sin(theta) * np.cos(psi)
        )
        - s_rot * np.cos(chi) * np.sin(theta) * np.sin(psi)
    )

    return bx, by, bz


def projected_coordinates(age_day, sign=+1, obs_time_day=0.0, light_time=True):
    """
    Projected coordinates on the sky.

    If light_time=False:
        simple projected position is used:

            beta_perp * c * age

    If light_time=True:
        light-travel-time effect is included:

            beta_perp * c * age / (1 - beta_x)
    """

    age_sec = age_day * day

    bx, by, bz = beta_components(
        age_day=age_day,
        sign=sign,
        obs_time_day=obs_time_day
    )

    if light_time:
        denom = 1.0 - bx
    else:
        denom = 1.0

    # projected physical positions
    y_app = by * c * age_sec / denom
    z_app = bz * c * age_sec / denom

    # angular offsets
    dec_arcsec = (y_app / D) * rad_to_arcsec
    ra_arcsec  = (z_app / D) * rad_to_arcsec

    return ra_arcsec, dec_arcsec, bx


def split_by_los(age_day, sign=+1, obs_time_day=0.0, light_time=True):
    """
    Split trajectory by line-of-sight velocity.

    Solid line:
        bx > 0, toward observer

    Dashed line:
        bx <= 0, away from observer

    NaN is inserted so matplotlib does not connect
    separated segments.
    """

    ra, dec, bx = projected_coordinates(
        age_day=age_day,
        sign=sign,
        obs_time_day=obs_time_day,
        light_time=light_time
    )

    mask_toward = bx > 0
    mask_away   = bx <= 0

    ra_toward  = np.where(mask_toward, ra, np.nan)
    dec_toward = np.where(mask_toward, dec, np.nan)

    ra_away  = np.where(mask_away, ra, np.nan)
    dec_away = np.where(mask_away, dec, np.nan)

    return ra_toward, dec_toward, ra_away, dec_away


def plot_jet_trajectory(ax, age_day, obs_time_day=0.0, light_time=True):
    """
    Plot approaching and receding jet trajectories on a given axis.
    """

    # approaching-side jet
    ra_a_toward, dec_a_toward, ra_a_away, dec_a_away = split_by_los(
        age_day=age_day,
        sign=+1,
        obs_time_day=obs_time_day,
        light_time=light_time
    )

    # receding-side jet
    ra_r_toward, dec_r_toward, ra_r_away, dec_r_away = split_by_los(
        age_day=age_day,
        sign=-1,
        obs_time_day=obs_time_day,
        light_time=light_time
    )

    # ------------------------------------------------------
    # approaching-side jet
    # ------------------------------------------------------
    ax.plot(
        ra_a_toward,
        dec_a_toward,
        color="tab:blue",
        lw=2.5,
        linestyle="-",
        label="approaching jet, toward"
    )

    ax.plot(
        ra_a_away,
        dec_a_away,
        color="tab:blue",
        lw=2.5,
        linestyle="--",
        alpha=0.8,
        label="approaching jet, away"
    )

    # ------------------------------------------------------
    # receding-side jet
    # ------------------------------------------------------
    ax.plot(
        ra_r_toward,
        dec_r_toward,
        color="tab:red",
        lw=2.5,
        linestyle="-",
        label="receding jet, toward"
    )

    ax.plot(
        ra_r_away,
        dec_r_away,
        color="tab:red",
        lw=2.5,
        linestyle="--",
        alpha=0.8,
        label="receding jet, away"
    )

    # central source
    ax.scatter(
        0,
        0,
        color="black",
        s=100,
        marker="*",
        zorder=10
    )

    # axis settings
    ax.set_aspect("equal")

    # astronomical convention:
    # RA increases to the left in sky images
    ax.invert_xaxis()

    ax.set_xlim(3.5, -3.5)
    ax.set_ylim(-1.5, 1.5)

    ax.xaxis.set_major_locator(MultipleLocator(0.5))
    ax.yaxis.set_major_locator(MultipleLocator(0.5))

    ax.grid(alpha=0.3)


# ==========================================================
# single-frame settings
# ==========================================================

# visible jet age range
age_day = np.linspace(0, 600, 3000)

# observation time
# 0 <= prec_phase < 1
# prec_phase = 0 means that the central approaching-side jet
# is closest to the observer in this convention.
prec_phase = 0.0
obs_time_day = P_prec_day * prec_phase


# ==========================================================
# plot comparison
# ==========================================================

fig, axes = plt.subplots(
    2,
    1,
    figsize=(8, 10),
    sharex=True,
    sharey=True
)

# upper panel: without light-travel-time effect
plot_jet_trajectory(
    ax=axes[0],
    age_day=age_day,
    obs_time_day=obs_time_day,
    light_time=False
)

axes[0].set_title(
    "Without light-travel-time effect"
)

# lower panel: with light-travel-time effect
plot_jet_trajectory(
    ax=axes[1],
    age_day=age_day,
    obs_time_day=obs_time_day,
    light_time=True
)

axes[1].set_title(
    "With light-travel-time effect"
)

# labels
axes[0].set_ylabel("Dec offset [arcsec]")
axes[1].set_ylabel("Dec offset [arcsec]")
axes[1].set_xlabel("RA offset [arcsec]")

# legend only once
axes[0].legend(fontsize=8, loc="upper left")

fig.suptitle(
    "SS 433 projected jet trajectory\n"
    f"precession phase = {prec_phase:.3f}, "
    f"t = {obs_time_day:.1f} day, "
    f"D = {D_kpc:.1f} kpc",
    y=0.98
)

plt.tight_layout()
plt.show()

図では、中心の黒い星印が SS 433 本体を表しています。
そこから左右に伸びている曲線が、歳差運動の履歴を反映したジェット軌道です。

比較のために、上段には光の到達時間効果を入れない場合、下段には光の到達時間効果を入れた場合を示しています。

上段は実際の見かけの軌道というより、単純に天球面上の速度成分に時間をかけた比較用の図です。
この場合、双方向に同じ速さで飛ぶジェットをそのまま投影しているため、青と赤のジェットは対称に見えます。

一方、下段では、視線方向速度によって光が届くタイミングが変わる効果を含めています。
その結果、こちらに向かう成分は見かけ上伸び、遠ざかる成分は相対的に縮んで見えます。

特に receding 側では、光の到達の遅れによって軌道が圧縮され、場所によっては構造が少し重なって見えます。
このように、光の到達時間効果を入れるだけでも、単純な対称性からずれた見かけの構造が現れます。

中心に近い部分は比較的最近放出されたジェット、外側はより過去に放出されたジェットに対応します。
そのため、ジェットの曲がり方には、過去数百日分の歳差運動の履歴が含まれています。

この記事では、歳差位相 $\psi=0$ を、approaching-side jet が最も観測者側に傾く位相として扱っています。
つまり、中心付近のジェットがこの位相にあるとき、approaching-side jet は視線方向に最も近づく向きを向いています。

なお、横軸は天文画像の慣例に合わせて、RA が左向きに増えるように反転しています。

アニメーションとして見る

次に、観測時刻を少しずつ変えて、1歳差周期分のアニメーションを作ります。

コードはこちら

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from matplotlib.ticker import MultipleLocator
from IPython.display import Image, display

# ==========================================================
# SS 433 projected jet trajectory
# one precession-cycle animation
# with light-travel-time effect
#
# age_day      : elapsed time after ejection
# obs_time_day : observation time
# t_emit_day   : ejection time
#
# t_emit_day = obs_time_day - age_day
# ==========================================================

# -----------------------------
# physical constants
# -----------------------------
c = 2.99792458e10          # [cm/s]
day = 86400.0              # [s]
kpc = 3.085677581e21       # [cm]
rad_to_arcsec = 206265.0

# -----------------------------
# SS 433 parameters
# -----------------------------
D_kpc = 5.5
D = D_kpc * kpc

beta = 0.26

# geometry
inc_deg = 78.83
theta_deg = 19.85
chi_deg = 98.2

# precession
P_prec_day = 162.375
s_rot = -1

# convert to rad
i = np.deg2rad(inc_deg)
theta = np.deg2rad(theta_deg)
chi = np.deg2rad(chi_deg)

# ==========================================================
# helper functions
# ==========================================================

def precession_phase(age_day, obs_time_day=0.0):
    """
    Precession phase of each jet element.

    A jet element observed with age_day was ejected at

        t_emit_day = obs_time_day - age_day

    Therefore, the jet direction should be evaluated
    at t_emit_day.
    """

    t_emit_day = obs_time_day - age_day
    psi = 2 * np.pi * t_emit_day / P_prec_day

    return psi


def beta_components(age_day, sign=+1, obs_time_day=0.0):
    """
    Velocity components in the observer frame.

    sign = +1 : approaching-side jet
    sign = -1 : receding-side jet

    x : line-of-sight direction
    y : Dec direction
    z : RA direction
    """

    psi = precession_phase(
        age_day=age_day,
        obs_time_day=obs_time_day
    )

    bx = sign * beta * (
        np.cos(i) * np.cos(theta)
        + np.sin(i) * np.sin(theta) * np.cos(psi)
    )

    by = sign * beta * (
        np.cos(chi)
        * (
            np.sin(i) * np.cos(theta)
            - np.cos(i) * np.sin(theta) * np.cos(psi)
        )
        + s_rot * np.sin(chi) * np.sin(theta) * np.sin(psi)
    )

    bz = sign * beta * (
        np.sin(chi)
        * (
            np.sin(i) * np.cos(theta)
            - np.cos(i) * np.sin(theta) * np.cos(psi)
        )
        - s_rot * np.cos(chi) * np.sin(theta) * np.sin(psi)
    )

    return bx, by, bz


def apparent_coordinates(age_day, sign=+1, obs_time_day=0.0):
    """
    Apparent coordinates on the sky.

    Light-travel-time effect is included as

        beta_perp * c * age / (1 - beta_x)
    """

    age_sec = age_day * day

    bx, by, bz = beta_components(
        age_day=age_day,
        sign=sign,
        obs_time_day=obs_time_day
    )

    # light-travel-time effect
    denom = 1.0 - bx

    # projected physical positions
    y_app = by * c * age_sec / denom
    z_app = bz * c * age_sec / denom

    # angular offsets
    dec_arcsec = (y_app / D) * rad_to_arcsec
    ra_arcsec  = (z_app / D) * rad_to_arcsec

    return ra_arcsec, dec_arcsec, bx


def split_by_los(age_day, sign=+1, obs_time_day=0.0):
    """
    Split trajectory by line-of-sight velocity.

    Solid line:
        bx > 0, toward observer

    Dashed line:
        bx <= 0, away from observer

    NaN is inserted so matplotlib does not connect
    separated segments.
    """

    ra, dec, bx = apparent_coordinates(
        age_day=age_day,
        sign=sign,
        obs_time_day=obs_time_day
    )

    mask_toward = bx > 0
    mask_away   = bx <= 0

    ra_toward  = np.where(mask_toward, ra, np.nan)
    dec_toward = np.where(mask_toward, dec, np.nan)

    ra_away  = np.where(mask_away, ra, np.nan)
    dec_away = np.where(mask_away, dec, np.nan)

    return ra_toward, dec_toward, ra_away, dec_away


# ==========================================================
# animation settings
# ==========================================================

# visible jet age range
age_day = np.linspace(0, 600, 3000)

# one full precession cycle
n_frames = 162
obs_time_grid = np.linspace(
    0,
    P_prec_day,
    n_frames,
    endpoint=False
)

gif_name = "ss433_one_precession_cycle.gif"

# ==========================================================
# figure
# ==========================================================

fig, ax = plt.subplots(figsize=(8, 4.5), dpi=70)

# initial observation time
obs_time_day = obs_time_grid[0]

# initial trajectories
ra_a_toward, dec_a_toward, ra_a_away, dec_a_away = split_by_los(
    age_day=age_day,
    sign=+1,
    obs_time_day=obs_time_day
)

ra_r_toward, dec_r_toward, ra_r_away, dec_r_away = split_by_los(
    age_day=age_day,
    sign=-1,
    obs_time_day=obs_time_day
)

# ----------------------------------------------------------
# line objects
# ----------------------------------------------------------

line_a_toward, = ax.plot(
    ra_a_toward,
    dec_a_toward,
    color="tab:blue",
    lw=2.5,
    linestyle="-",
    label="approaching jet, toward"
)

line_a_away, = ax.plot(
    ra_a_away,
    dec_a_away,
    color="tab:blue",
    lw=2.5,
    linestyle=":",
    alpha=0.8,
    label="approaching jet, away"
)

line_r_toward, = ax.plot(
    ra_r_toward,
    dec_r_toward,
    color="tab:red",
    lw=2.5,
    linestyle="-",
    label="receding jet, toward"
)

line_r_away, = ax.plot(
    ra_r_away,
    dec_r_away,
    color="tab:red",
    lw=2.5,
    linestyle=":",
    alpha=0.8,
    label="receding jet, away"
)

# central source
ax.scatter(
    0,
    0,
    color="black",
    s=120,
    marker="*",
    zorder=10
)

# ==========================================================
# axis settings
# ==========================================================

ax.set_xlabel("RA offset [arcsec]")
ax.set_ylabel("Dec offset [arcsec]")

ax.set_aspect("equal")

# astronomical convention:
# RA increases to the left in sky images
ax.invert_xaxis()

# display range
ax.set_xlim(3.5, -3.5)
ax.set_ylim(-1.5, 1.5)

# 0.5 arcsec tick spacing
ax.xaxis.set_major_locator(MultipleLocator(0.5))
ax.yaxis.set_major_locator(MultipleLocator(0.5))

ax.grid(alpha=0.3)
ax.legend(fontsize=8, loc="upper left")

plt.tight_layout()

# ==========================================================
# update function
# ==========================================================

def update(frame):

    obs_time_day = obs_time_grid[frame]
    prec_phase = (obs_time_day / P_prec_day) % 1.0

    # approaching-side jet
    ra_a_toward, dec_a_toward, ra_a_away, dec_a_away = split_by_los(
        age_day=age_day,
        sign=+1,
        obs_time_day=obs_time_day
    )

    # receding-side jet
    ra_r_toward, dec_r_toward, ra_r_away, dec_r_away = split_by_los(
        age_day=age_day,
        sign=-1,
        obs_time_day=obs_time_day
    )

    line_a_toward.set_data(ra_a_toward, dec_a_toward)
    line_a_away.set_data(ra_a_away, dec_a_away)

    line_r_toward.set_data(ra_r_toward, dec_r_toward)
    line_r_away.set_data(ra_r_away, dec_r_away)

    ax.set_title(
        "SS 433 projected jet trajectory\n"
        f"one precession cycle: phase = {prec_phase:.3f}, "
        f"t = {obs_time_day:.1f} day"
    )

    return (
        line_a_toward,
        line_a_away,
        line_r_toward,
        line_r_away
    )

# ==========================================================
# create and save animation
# ==========================================================

anim = FuncAnimation(
    fig,
    update,
    frames=n_frames,
    interval=80,
    blit=True
)

anim.save(
    gif_name,
    writer=PillowWriter(fps=15)
)

plt.close(fig)

print(f"Saved: {gif_name}")

# ==========================================================
# display in Google Colab
# ==========================================================

display(Image(filename=gif_name))

観測時刻を変えると、中心付近に新しく放出されるジェットの向きが変わります。
それと同時に、外側の構造も、異なる放出時刻のジェットとして更新されます。

その結果、空の上のジェット軌道が、歳差周期に従ってゆっくり変化していく様子を見ることができます。

今回の可視化で見えてくること

今回の可視化では、SS 433 のジェットが空の上でどのように見えるかを、歳差運動と光の到達時間効果を入れて眺めてみました。

三次元空間でのジェットの向きは歳差運動によって決まりますが、観測される見かけの位置には、視線方向速度による光の到達時間効果も関わります。

また、物理的な位置を RA/Dec offset の arcsec スケールに変換するためには、SS 433 までの距離も仮定する必要があります。

今回の図は、厳密な WCS 変換を行った観測座標ではなく、SS 433 近傍の局所的な小角度近似に基づくモデル軌道です。
それでも、歳差運動や光の到達時間効果によって、ジェットの見かけの形がどのように変化するのかを見るには十分興味深いと思います。

特に、実際の観測で見えているジェット構造を考えるときには、このような幾何学的・相対論的な効果も含めて考える必要があります。

おわりに

今回は、SS 433 の歳差ジェットを空の上に投影し、光の到達時間効果を含めた見かけの軌道を可視化しました。

前回の記事では、ジェットの向きを三次元ベクトルとして表しました。
今回はそこから一歩進めて、観測者が実際に見る二次元の構造として描いてみました。

特に、こちらに向かう成分を持つ部分は見かけ上伸び、遠ざかる成分を持つ部分は相対的に縮んで見える、という点が重要です。

ただし、今回扱ったのはかなりシンプルな歳差モデルです。
実際の SS 433 では、約162日の歳差運動に加えて、数日スケールの nodding motion も知られています。最近の XRISM 観測でも、X線スペクトルのドップラーシフトという観点から、このような短い時間スケールの速度変動が議論されています（e.g., Shidatsu et al. 2025; Sakai et al. 2026）。

つまり、ジェット方向は単純な歳差円錐上をなめらかに回るだけではなく、さらに小さく揺れるような運動も含んでいます。
このような運動を含めた軌道モデルについては、たとえば Collins et al. (2002) などが参考になります。

私自身もこのあたりはまだ整理しきれていない部分がありますが、実際の観測で見えている複雑なジェット構造を理解するためには、こうした時間変動や幾何学的効果を少しずつ含めて考えていく必要があるのだと思います。

今回はまず、基本となる歳差運動と光の到達時間効果だけに絞って、SS 433 の見かけのジェット軌道を可視化してみました。

参考