はじめに

XSPEC では flux コマンドに err オプションを付けることで、フラックスの 信頼区間（上下限） を計算できます。対話的解析ではおなじみですが、

PyXspec からこの結果をどう受け取ればいいのか？

については、意外とまとまった情報がありません。
この記事では、PyXspec で flux err を実行し、その結果を Python 側で受け取る方法を簡単に整理します。

前提となる最小コード例

以下のように、スペクトルを読み込み、モデルを定義してフィットが終わっている状態を前提とします。

import xspec as xs

spec = xs.Spectrum("example.pha")
model = xs.Model("powerlaw")

xs.Fit.perform()

（※モデルの種類は本記事の主題ではないため、ここでは一例として powerlaw を用いています。）

ここで押さえておくのは次の点だけです。

xs は PyXspec モジュール
spec は Spectrum オブジェクト
flux コマンドの結果は この spec に保存される

基本：XSPEC の `flux err` 出力

XSPEC で次のように打つと、

XSPEC> flux 2.0 7.0 err 1000 90

標準出力には次のような結果が表示されます。

出力例

Parameter distribution is derived from fit covariance matrix.
 Model Flux 0.0047934 photons (2.8108e-11 ergs/cm^2/s) range (2.0000 - 7.0000 keV)
     Error range  0.004750 - 0.004844    (2.783e-11 - 2.843e-11)  (90.00% confidence)

ここで XSPEC は、

photon flux（値＋下限＋上限）
energy flux（値＋下限＋上限）

を 同時に計算しています。
このコマンドは一例ですが、「2–7 keV の範囲で、1000 回の Monte Carlo を行い、90% 信頼区間の誤差を求める」という意味です。

誤差計算のための Monte Carlo のデフォルト値について

本記事では誤差計算の詳細な方法については割愛しますが、flux err の Monte Carlo 回数を省略した場合、XSPEC ではデフォルトで 100 回の試行が行われているようです（「XSPEC の flux err のデフォルト挙動が気になったのでソースコードを確認したメモ」参照）。

ただし、試行回数が少ない場合には上下限のばらつきが大きくなることがあるため、解析の目的に応じて明示的に回数を指定する方が安全です。

PyXspec での実行方法

PyXspec でも同じことができます。

import xspec as xs

xs.AllModels.calcFlux("2.0 7.0 err 1000 90")

このコマンドは XSPEC 内部で flux err を実行し、その結果を実行対象となっている Spectrum オブジェクトに保存します。
PyXspec では、この結果を Spectrum.flux 属性から取得できます。

ここで指定しているエネルギー範囲（2.0–7.0 keV）、Monte Carlo 回数、信頼区間はいずれも一例です。解析の目的に応じて適宜変更してください。

`spec.flux` 属性に入っているもの

flux コマンド実行後、spec.flux には 6 つの要素が並んでおり、以下のように受け取ることができます。

flux_erg, flux_erg_lo, flux_erg_hi, flux_ph, flux_ph_lo, flux_ph_hi = spec.flux

spec.flux の並びは 次の順序で固定されています。

index	内容	単位
0	energy flux	erg/cm²/s
1	energy flux (lower)	erg/cm²/s
2	energy flux (upper)	erg/cm²/s
3	photon flux	photons/cm²/s
4	photon flux (lower)	photons/cm²/s
5	photon flux (upper)	photons/cm²/s

つまり、

energy flux → photon flux の順で、値・下限・上限

です。

XSPEC の標準出力では photons → ergs の順に表示されるため、直感と逆になりやすい点には注意が必要です。

実用的な使い方例

実務向けに、結果を見やすい形で出力する一例を示します。

flux_erg, flux_erg_lo, flux_erg_hi, flux_ph, flux_ph_lo, flux_ph_hi = spec.flux

print(f"Energy flux: {flux_erg:.3e} "
      f"({flux_erg_lo:.3e} – {flux_erg_hi:.3e}) erg/cm^2/s")

print(f"Photon flux: {flux_ph:.3e} "
      f"({flux_ph_lo:.3e} – {flux_ph_hi:.3e}) photons/cm^2/s")

# ± 表記（非対称誤差）
flux_erg_err_minus = flux_erg - flux_erg_lo
flux_erg_err_plus  = flux_erg_hi - flux_erg

flux_ph_err_minus = flux_ph - flux_ph_lo
flux_ph_err_plus  = flux_ph_hi - flux_ph

print(
    f"Energy flux: {flux_erg:.3e} "
    f"-{flux_erg_err_minus:.3e} "
    f"+{flux_erg_err_plus:.3e} erg/cm^2/s"
)

print(
    f"Photon flux: {flux_ph:.3e} "
    f"-{flux_ph_err_minus:.3e} "
    f"+{flux_ph_err_plus:.3e} photons/cm^2/s"
)

出力例

Energy flux: 2.811e-11 (2.783e-11 – 2.843e-11) erg/cm^2/s
Photon flux: 4.793e-03 (4.750e-03 – 4.844e-03) photons/cm^2/s
Energy flux: 2.811e-11 -2.799e-13 +3.215e-13 erg/cm^2/s
Photon flux: 4.793e-03 -4.363e-05 +5.047e-05 photons/cm^2/s

補足：モデルの成分ごとの flux を求めるには

本記事で扱った flux は、モデル全体に対するフラックスです。
モデルの特定成分（例：powerlaw 成分のみ）のフラックスを求めたい場合は、XSPEC の cflux コンポーネントを使用します。PyXspec でも同様に、モデルに組み込むことで扱えます。

cflux はモデルの一部として組み込んで使用するため、flux コマンドとは設計思想が異なります。本記事の範囲を超えるため詳細は扱いませんが、誤差計算の考え方自体は通常のモデルパラメータ誤差と同じです。
PyXspec におけるモデルパラメータの誤差計算については、以下の記事を参照していただけると幸いです。

まとめ

PyXspec では、XSPEC と同じ書式で flux err を実行でき、その結果は Spectrum.flux から取得できます。
この際、energy flux と photon flux の両方が返り、順序が固定されているため、単位と並びには注意が必要です。

参考