Swift/BAT 15–50 keV モニター入門：Crab 単位と daily/orbit ライトカーブを Python で可視化する

Last updated at 2026-02-20Posted at 2024-10-13

はじめに

本記事では、Neil Gehrels Swift Observatory に搭載された Burst Alert Telescope (BAT) の 15–50 keV バンドの公開モニターデータからライトカーブを作成する方法を、matplotlib と plotly で紹介する。

対象は主に大学院生レベルを想定し、

Swift/BAT モニターが物理的に何を測っているのか
「count rate」や「Crab 単位」が物理量として何を意味するのか
FITS 形式（daily / orbit）の違いによるハマりどころ

も含めて、実務で困らない程度の背景を整理する。

公式情報や代表的な文献は以下を参照されたい：

公式の BAT Transient Monitor

Swift ミッション概要（NASA）

BAT 装置論文
Barthelmy et al. 2005, Space Sci. Rev., 120, 143

BAT Transient Monitor 概要
Krimm et al. 2013, ApJS, 209, 14

日本語による Swift/BAT 解説（例）
青山学院大学坂本研究室「ガンマ線バースト探査衛星」

Swift 衛星の軌道はこちらから確認できます。

Perigee: 399.7 km
Apogee: 403.3 km
Period: 92.4 minutes
Launch date: November 20, 2004

なので、だいたい 1 orbit が 90分の低軌道(LEO)衛星です。

本記事は「可視化」を主目的とする。
絶対フラックス換算や応答行列を用いた厳密解析は、上記文献および HEASoft/XSPEC のドキュメントを参照されたい。

プログラムの基礎部分は下記記事と同様である：

Google Colab でも実行できるノートブック例：

Swift/BAT と 15–50 keV モニターの位置づけ

Swift 衛星とは

Swift は 2004 年打ち上げの多波長衛星で、主目的はガンマ線バースト（GRB）の迅速追跡観測である。搭載装置は主に：

BAT（15–150 keV 領域の広視野硬 X 線検出器、コーデッドマスク方式）
XRT（0.3–10 keV の集光型 X 線望遠鏡）
UVOT（紫外～可視の望遠鏡）

本記事で扱うのは、BAT が提供している 15–50 keV モニターバンドの長期ライトカーブである。

BAT モニターの性質と限界

BAT は

広視野（空のかなり広い部分を常時監視）
コーデッドマスク方式（マスクパターンから像を逆問題として再構成）
バックグランド支配型検出（低フラックス源ではノイズは検出器＋宇宙線バックグラウンドで決まる）

という性質を持つ。そのため、BAT トランジェントモニターは

「長期変動（状態遷移・アウトバースト）」の検出・追跡に特化した硬 X 線モニタ

と理解するのがよい。個々のスペクトルを精密に議論したい場合は、同時期の XRT や他衛星（MAXI, XRISM, NuSTAR など）と組み合わせる必要がある。

また、

Telescope PSF 17 arcmin

が目安であり、この範囲内に２つの点源がいると区別はできないことに注意。あるいは、近傍に明るい天体がいる場合も漏れこみには注意が必要(空間的な漏れ込みはBATに限った話ではなく、いつでも確認が必要。)

count rate と物理量の関係

公開されている 15–50 keV モニターデータの RATE は概ね

背景差し引き
マスク重み付け後

の count rate（通常は count/cm^2/s）として与えられている。

一般的には

C = \int_{15}^{50} R(E)\, F(E)\, dE

$F(E)$：光子フラックス [ph cm⁻² s⁻¹ keV⁻¹]
$R(E)$：有効面積・応答を含んだ検出効率 [cm²]
$C$：観測される count rate [count s⁻¹] または [count cm⁻² s⁻¹]

であり、BAT モニターで扱っている RATE はこの $C$ に対応する。

重要な点は、

「BAT の RATE は、応答関数 $R(E)$ で重み付けされた「観測量」であり、
そのまま理論スペクトルと 1 対 1 に比較できるエネルギーフラックス（erg cm⁻² s⁻¹）ではない」

ということである。

エネルギーフラックスに換算したい場合は、

スペクトル形（典型的にはパワーロー＋指数カットオフ等）
BAT の応答行列（RMF/ARF）

を仮定し、HEASoft/XSPEC 等で評価する必要がある。

あるいは、近似的な方法としては、1光子あたりのエネルギーを仮定し、距離を仮定して光度に戻す。

Crab 単位とは何か？

硬 X 線・ガンマ線の観測では、明るく比較的安定である Crab 星雲を基準天体として

1 Crab
1 mCrab（＝10⁻³ Crab）

といった単位が広く用いられている。装置ごとに有効面積やエネルギーバンドが異なるため、Crab 単位で表現することで「装置間のざっくりした比較」がしやすくなる。

Swift/BAT 15–50 keV における代表値

BAT Transient Monitor の Crab ライトカーブや Krimm et al. (2013) の図から読み取ると、15–50 keV バンドではおおよそ：

$C_{\rm Crab} \sim 0.22\ \rm{count~cm^{-2}~s^{-1}}$
対応するエネルギーフラックスは
$F_{\rm Crab} \sim 1.3\times10^{-8}\ \rm{erg~cm^{-2}~s^{-1}}$
（Crab スペクトル仮定に依存する近似値）

程度とみなされている（オーダーとしての代表値）。

fig2 より

実際に BAT Transient Monitor の Crab ページ：

の 15–50 keV グラフを眺めると、明るい状態で 0.2 count/cm²/s 程度であることが確認できる。

したがって、この記事では便宜的に

1\ \mathrm{Crab} \simeq 0.22\ \mathrm{count~cm^{-2}~s^{-1}}

とし、

1\ \mathrm{mCrab} = 10^{-3}\ \mathrm{Crab} \simeq 2.2\times10^{-4}\ \mathrm{count~cm^{-2}~s^{-1}}

という換算を用いる。

注意）
この値は「代表値」であり、厳密なフラックスやスペクトル解析には BAT の最新キャリブレーションと応答行列を用いた評価が必要である。

データのダウンロードと FITS 形式の違い（daily / orbit）

ダウンロード方法

公式ページ：

あるいは以下の Python スクリプトでも取得できる：

対象天体（例：Cyg X-1）の 15–50 keV モニターファイルをダウンロードすると、通常は

CygX-1.lc.fits（daily light curve）
CygX-1_orbit.lc.fits（orbit light curve）など

2 種類の FITS が存在することがある。

例えば、下記のように Daily light curves をクリックすると、

CygX-1.lc.fits というFITS形式のファイルがダウンロードされるはずです。

ここが Swift/BAT のちょっとしたハマりどころで、FITS ヘッダの定義が daily/orbit で微妙に異なる。

daily light curve（*.lc.fits）

Daily LC の例（RapidBurster.lc.fits）では、ヘッダは概ね：

TIMEUNIT = 'd'
TIMESYS = 'MJD'
BINTABLE カラム：
- TIME：MJD（Modified Julian Date）
- RATE：count/cm²/s
- ERROR：その誤差
- 別途 MJD カラムは無い

すなわち、

daily LC では TIME カラムがそのまま MJD（日）

として解釈できる。

orbit light curve（*.orbit.lc.fits）

一方、Orbit LC の例（RapidBurster.orbit.lc.fits）では：

TIMEUNIT = 's'
TIMESYS = 'TT'
ヘッダに MJDREFI, MJDREFF が存在
BINTABLE カラム：
- TIME：MET（Mission Elapsed Time；MJDREF からの経過秒）
- MJD：別カラムとして MJD（日）が存在
- RATE, ERROR：定義は同様

つまり、

orbit LC では TIME は MET（秒）であり、そのまま MJD として使ってはいけない。
MJD を使うなら、MJD カラムを使うか、MJDREFI/MJDREFF と TIME から変換する必要がある。

本記事の可視化スクリプトでは、

BINTABLE のカラム名を調べ、MJD カラムの有無で
- MJD が無ければ daily LC → TIME を MJD と見なす
- MJD があれば orbit LC → TIME を MET と見なし、MJDREFI/MJDREFF を使って MJD に変換

というロジックで daily/orbit を自動判定している（詳細は Appendix 参照）。

時間変換は

\mathrm{MJD} = \mathrm{MJDREFI} + \mathrm{MJDREFF} + \frac{\mathrm{TIME} + \mathrm{TIMEZERO}}{86400}

という近似式を用いている。
MJDREF の定義や TT/UTC の厳密な扱い（うるう秒など）は長期ライトカーブの可視化レベルでは無視しうる近似である。

プロット方法

使用スクリプト：

このスクリプトは、以下を行う：

daily / orbit LC を自動判定して読み込み
RATE, ERROR に対して
- rate > 0、NaN/Inf 除去
- 相対誤差カット（ERROR/RATE <= threshold、デフォルト 50%）
matplotlib または plotly で 15–50 keV ライトカーブを描画
--yscale により linear/log を切り替え
matplotlib の場合、右軸に Crab 単位（1 Crab = 0.22 count/cm²/s）を表示
--bin_days により 10, 30, 60 日などの重い時間ビニングも実行可能

matplotlib を使う場合（論文向き）

縦軸リニア（線形）

python swift_plot_lc.py CygX-1.lc.fits \
  --plotter matplotlib \
  --highlight 2023-11-03T23:51:04 2023-12-05T14:01:04 \
  --yscale linear \
  --relerr_thresh 80

左軸：count/cm²/s（15–50 keV）
右軸：Crab 単位（1 Crab = 0.22 count/cm²/s と仮定）
fill_between による期間ハイライト
出力：CygX-1.lc_matplotlib.png（など）

縦軸ログ

python swift_plot_lc.py CygX-1.lc.fits \
  --plotter matplotlib \
  --highlight 2023-11-03T23:51:04 2023-12-05T14:01:04 \
  --yscale log \
  --relerr_thresh 80

縦軸を log スケールにした例
matplotlib.secondary_yaxis を使って、log スケールでも 右軸の Crab 単位が左軸と正しく対応するように調整している

暗い天体（Rapid Burster など）では、--relerr_thresh を 80–100 程度まで緩めた方がデータ点が残りやすい。

plotly を使う場合（対話型）

python swift_plot_lc.py CygX-1.lc.fits \
  --plotter plotly \
  --highlight 2023-11-03T23:51:04 2023-12-05T14:01:04 \
  --yscale linear \
  --relerr_thresh 80

ブラウザ上でズーム・マウスオーバー可能なインタラクティブ表示
HTML として保存 (*_plotly.html)

重い時間ビニング（10, 30, 60 日など）

暗い天体では、daily LC ではほとんどがノイズに埋もれてしまう。その場合は

python swift_plot_lc.py RapidBurster.lc.fits \
  --plotter matplotlib \
  --yscale linear \
  --relerr_thresh 80 \
  --bin_days 10 30 60

のように --bin_days オプションを使うことで、

10 日
30 日
60 日

などのビニングを適用したライトカーブを追加で生成することができる（各ビン内では $1/\sigma^2$ 重み付き平均を取っている）。

有名なブラックホール連星 Cyg X-1 の 15–50 keV 帯で、low/hard 状態のときに BAT モニター上で 1 Crab 程度（RATE≈0.2 count/cm²/s） で見えていることを確認してみると、Crab 単位の感覚が掴みやすい。

研究利用時の注意点

BAT モニターは「広帯域・広視野・長期変動」に強いが、個々のスペクトル形状を議論するには向かない
RATE は応答関数込みの量であり、そのまま理論フラックス（erg cm⁻² s⁻¹）とは比較できない
mCrab 表記を用いる場合は、必ず
- 「BAT 15–50 keV での mCrab」
- もしくは「Swift/BAT 15–50 keV」
  と、装置とエネルギーバンドを明示すること
TIME が MJD か MET かは FITS 形式（daily/orbit）によって異なるので、ヘッダの TIMEUNIT, TIMESYS, MJDREFI, MJDREFF を必ず確認する

本記事のスクリプトと Appendix は、こうした罠をある程度吸収するよう設計しているが、研究利用時には各々の天体の原著論文やファイルのヘッダなどを必ず自分で確認することを強く推奨する。

まとめ

Swift/BAT 15–50 keV モニターは、硬 X 線変動天体の長期追跡に非常に有用
公開データの RATE は count rate（応答込みの観測量）であり、物理量への換算にはスペクトル仮定と応答行列が必要
代表値として、BAT 15–50 keV において 1 Crab ≈ 0.22 count/cm²/s として扱える
FITS 形式は daily / orbit で異なり、TIME が MJD か MET かはヘッダを見ないと分からない
可視化自体は Python（matplotlib, plotly）で比較的簡単に実行できるが、背後の物理と FITS 仕様を理解しておくとハマりどころを避けやすい

Appendix: Swift/BAT ライトカーブ可視化コードの概要

ここでは、本記事で用いた Python スクリプト

の中で、特に

FITS データ形式（daily / orbit）の扱い
load_and_filter_data() によるフィルタリング
時間ビニング
Crab 副軸と --yscale の関係
CLI オプションの意味

を簡単に解説する。

A.1 daily / orbit 形式の自動判定

load_and_filter_data() の冒頭では、まず BINTABLE のカラム名を調べて

colnames = [name.upper() for name in data.names]

if 'MJD' not in colnames:
    # daily LC: TIME = MJD [d]
else:
    # orbit LC: TIME = MET [s], MJDREFI/MJDREFF を使って MJD を計算

という形で daily/orbit を判定している。

daily LC
MJD カラムが存在しない → TIME が MJD（TIMEUNIT='d', TIMESYS='MJD'）
→ mjd = data['TIME'] としてよい。
orbit LC
MJD カラムが存在する → TIME は MET（秒）
ヘッダにある MJDREFI, MJDREFF, TIMEZERO を用いて
```
\mathrm{MJD} = \mathrm{MJDREFI} + \mathrm{MJDREFF} + \frac{\mathrm{TIME} + \mathrm{TIMEZERO}}{86400}
```
として近似的に MJD に変換している。

変換後の MJD は astropy.time.Time(mjd, format='mjd') を通じて datetime へ変換し、プロット用の横軸（dates）として用いている。

A.2 `load_and_filter_data()` のフィルタリング手順

関数の概略は以下の通り：

FITS を開き、daily/orbit を判定
mjd, rate, error を NumPy 配列に変換

有限値かつ正の flux であるものを残す

finite_mask = np.isfinite(rate) & np.isfinite(error)
pos_mask    = rate > 0
base_mask   = finite_mask & pos_mask

相対誤差フィルタ（オプション）を適用
```
relative_error = error / rate
mask = (relative_error <= relative_error_threshold/100.0)
```
※ --no_relerr_filter を付けた場合はこのステップをスキップ
フィルタ後の dates, mjd, rate, error をまとめて返す

このとき、標準出力には

データ形式（daily/orbit）
MJD の範囲
RATE の min/max/mean/median
rate<=0 の個数
相対誤差の min/max
フィルタ前後のデータ数

などが必ず表示されるようにしてあり、どのフィルタでデータが消えているかが簡単にトレースできるようになっている。

下記が Cyg X-1 で実行した時の標準出力の例。

python swift_plot_lc.py CygX-1.lc.fits --plotter matplotlib --yscale linear --relerr_thresh 80   --highlight 2023-11-03T23:51:04 2023-12-05T14:01:04
=== load_and_filter_data ===
filename            : CygX-1.lc.fits
relerr_threshold[%] : 80.0
apply_relerr_filter : True
columns             : ['TIME', 'RATE', 'ERROR', 'YEAR', 'DAY', 'STAT_ERR', 'SYS_ERR', 'DATA_FLAG', 'TIMEDEL_EXPO', 'TIMEDEL_CODED', 'TIMEDEL_DITHERED']
mode                : one-day (TIME = MJD [d])
N total rows        : 7112
MJD range           : 53416.000 – 61068.000
RATE stats          : min=-2.094e-03, max=7.151e-01, mean=1.286e-01, median=1.446e-01
ERROR stats         : min=1.801e-05, max=5.394e-02, mean=5.409e-03
N rate<=0           : 3
N finite(rate,err)  : 7112
N rate>0            : 7109
N after finite&>0   : 7109
Rel.err stats       : min=0.035, max=8.961
N rel.err <= 0.800 : 7108 / 7109
N after relerr filt.: 7108
=== end of load_and_filter_data ===
Saving matplotlib figure to CygX-1.lc_matplotlib.png

--debug オプションを付けると、相対誤差の percentile（P10, P50, P90）など、もう少し詳細な統計も表示する。

A.3 相対誤差フィルタの物理的意味

相対誤差

\frac{\sigma}{F} = \frac{\mathrm{ERROR}}{\mathrm{RATE}}

が大きい点は、count が少なく S/N が低いデータ点である。可視化目的では

「相対誤差が例えば 50% を超える点は、見た目のノイズとして扱って捨てる」

という簡易な品質カットがよく使われる。

本スクリプトの --relerr_thresh はこのしきい値（%）であり、例えば

--relerr_thresh 50 → ERROR/RATE <= 0.5 を採用
--relerr_thresh 80 → ERROR/RATE <= 0.8 を採用

といった形で制御している。
暗い天体では、しきい値を 80–100% 程度まで緩めないとデータがほぼ全滅することもある。

論文でこのフィルタを使う場合は、

「どのしきい値を採用したか」
「それによってどの程度の fraction のデータが残っているか」

を本文中で明示しておくとよい。

A.4 時間ビニング（`bin_flux_by_days`）

bin_flux_by_days() 関数は、指定した日数 bin_size_days で時間ビニングしたライトカーブを生成する。

ビン境界を np.arange(t_min, t_max + bin_size_days, bin_size_days) で定義
各ビン内の rate は誤差 $\sigma$ を使って $1/\sigma^2$ 重み付き平均
```
\bar{F} = \frac{\sum_i F_i / \sigma_i^2}{\sum_i 1/\sigma_i^2}
```

ビン内の誤差は

\bar{\sigma} = \frac{1}{\sqrt{\sum_i 1/\sigma_i^2}}

という形で計算している。ビン中心の時刻は、そのビンに含まれる MJD の単純平均である。

この重み付けにより、誤差が小さい点（S/N が高い点）がより強く効くようになり、長期ビニングでも統計的に自然な平均値が得られる。

A.5 y 軸スケールと Crab 副軸

create_matplotlib_subplots() では

左軸：rate [count/cm²/s]
右軸：Crab 単位（1 Crab = 0.22 count/cm²/s）

を同時に表示するために matplotlib.secondary_yaxis を使用している。

CRAB_RATE_15_50_KEV = 0.22

def counts_to_crab(y_counts):
    return y_counts / CRAB_RATE_15_50_KEV

def crab_to_counts(y_crab):
    return y_crab * CRAB_RATE_15_50_KEV

ax.set_yscale(yscale)  # 'linear' or 'log'

ax_crab = ax.secondary_yaxis(
    'right',
    functions=(counts_to_crab, crab_to_counts),
)
ax_crab.set_ylabel('Flux [Crab] (15–50 keV)')

これにより、

左軸を --yscale linear / --yscale log で切り替えても
右軸は 常に同じ y 位置が「同じ物理量」を表すように 自動で変換される

ようになっている。twinx で自前同期するよりも安全で、ズームや自動スケーリングにも追従する。

A.6 コマンドラインオプションのまとめ

主なオプションを表にしておく：

オプション	意味
`filename`	入力 FITS ファイル（daily/orbit どちらでも可）
`--plotter`	`matplotlib` / `plotly` の選択（デフォルト：matplotlib）
`--start_date`	プロットする開始日時（ISOT 形式）
`--end_date`	プロットする終了日時（ISOT 形式）
`--highlight START END`	指定期間をハイライト（複数指定可）
`--relerr_thresh`	相対誤差しきい値 [%]（デフォルト 50）
`--no_relerr_filter`	付けると相対誤差フィルタを無効化
`--yscale`	`linear` / `log`（デフォルト：linear）
`--bin_days`	重い時間ビニング（日単位）。例：`--bin_days 10 30 60`
`--debug`	相対誤差分布など、追加の統計情報を表示

実際の解析では、対象天体の明るさや目的に応じて

--relerr_thresh を調整し、
faint source では --bin_days を併用、
表示の好みに応じて --yscale を linear / log で切り替える

といったチューニングを行うことになる。

まとめ

Swift/BAT の 15–50 keV モニターは、「Crab 単位でざっくり明るさを見る」「daily/orbit の違いを理解したうえで長期変動を追う」といった、硬 X 線観測の実務感覚を身につけるうえで非常に良い教材である。本記事のスクリプトと解説を足がかりに、たとえば BAT と MAXI を重ねてプロットして状態遷移を比較したり、特定の期間を切り出して XRT や XRISM のスペクトル解析と結びつけたり、といった発展的な解析にもぜひ挑戦してほしい。Swift/BAT のライトカーブを「安全に・正しく」読めるようになることは、次世代ミッションを創造することに直結するはずと期待。

Swift XRT の lightcurve は下記から生成できます。

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up