最終更新: 2020-11-07 (榎戸輝揚)
lecture
,xspec
X線のスペクトル解析を行うパッケージはいくつかあり、その代表的な一つが Xspec です。
この Xspec は、NASA を中心に開発されてきた高エネルギー天文学のデータ解析の統合パッケージ "HEASoft" に含まれており、HEASoft をインストールする必要があります。
に従って、インストールできます。(昔はいろいろハマるところがありましたが、最近は、書かれている通りにやればインストールできます)
インストールして、Xspec を呼び出すことができるようになると、コマンドラインで対話的にスペクトルのフィットを行うことができます。もちろん、バッチ処理にすることもできます。
以下、コマンドラインを使える人は実際に自分でやってもらえると良いのですが、そうでない方は以下に動画を作っていますので、これを見てもらえればと思います。
まず練習用のディレクトリを作ります。
必要なファイル一式は、GitHub: lecture_xspec に置いています。以下のようにダウンロードできます。GitHub が使える前提ですが、使えない方は zip file でダウンロードできるようにしています。
コマンドラインで以下のように入力してみましょう。EOF はヒアドキュメントと呼ばれていて、実際に自分でコマンドプロンプトから入力する際には、EOF は無視して順番に入れてみてください。
これで表示されるのは、検出器で測定された「かに星雲」(Crab Nebula)のスペクトルです(注: ここでは実習のために、シミュレーションで作成したものです)。かに星雲は、X線天文衛星が打ち上げ後に真っ先に観測する標準光源としての意味のある天体ですね。
データは、 検出器で計数されたX線光子のカウント数をエネルギーに対してプロットしたもので、縦軸は Counts s keV で表示される。XSPEC の plot コマンドでは、誤差棒は Poisson 統計から予想される 1σ が default ではつくようになっています。
このとき、検出器に入射するスペクトル情報を引き出したいわけです。この観測データを再現するモデル関数として、天体が「べき」関数で表せる X 線スペ クトルをもち、星間吸収を受けて検出器まで届いたと考える。これは「べき」関数 (power) と星間吸収 (phabs) の積 phabs * pow で与えられるとし、検出器に入射したエネルギー E (keV) を もつ X 線光子の個数は、
photons cm s keV
と表せる。星間吸収において、 atoms cm は水素柱密度、 cm は光電吸収の断面積であり、「べき」関数において、K Photons cm−2 s−1 keV−1 は 1 keV で規格化した強度 (Normalization; Norm)、Γ は光子指数である。初期パラメータを入れてみる。
まず初期パラメータである程度、実験データに近いものを見つけたあとで、実験データを最もよく再現するモデルパラメータの最良推定値を探してみよう。
の最小値を求めるアルゴリズムには複数あり、XSPEC のデフォルトは Bevington の教科書で
の CURFIT ルーチンに基づいた Levenberg-Marquardt 法 (leven) を用いている。Levenberg-Marquardt 法はの関数面の傾きから最小値を探索する最急降下法 (Steepest descent method)とガウス-ニュートン法 (Gauss-Newton method)を合わせたようなものである。このほかにCERN ライブラリの MINUIT 関数も用いることができ、
などとして切り替えることができる。method としては、migrad (MINUIT の MIGRAD), simplex (SIMPLEX)が用意されている。を最小化するようにモデル関数のパラメータを最適化しデータに合わせ込むには、
とする。これで、200回の試行を繰り返しを行い (iteration)、の最小値を見つけ出す。収束した場合、下の図のように、データによく合うモデル関数(赤)が得られ、中図のように、両者の残差は 0 を中心として分布する。
図: X線スペクトルデータ(黒点)とスペクトルのモデル関数(赤)。(top)検出器で計数された生のスペクトルと最適なモデル関数の解。(middle) モデルに対するデータの比。(data-model)/error が表示されている。(bottom) 検出器の応答関数を考慮して算出される、検出器に入射する前の段階での天体のスペクトルをで表示している。星間吸収は含まれている。
試行回数を限定せずにフィットするには、
などとする。フィットが収束すると、以下のような最適なパラメータと統計量が表示される。
これを読むと、66 ビンのデータを 3 パラメータのモデルでフィットしたので、自由度 は 66-3=63 であり、得られた であるから、フィットの妥当性の目安として用いられる reduced chi-square の値は 程度となったことがわかる。
分布において、上記の自由度でが得られる確率は、0.82 (Null hypothesis probability) であるから、もし 仮に 100 回同一の観測を行えたとすれば、このモデルに合致する観測結果が 82 回得られる、つまり最適なパラメータをもつモデル関数はデータをよく再現すると言える。このようにして妥当なモデル関数を得ることができた。
一般に、 程度の範囲であればモデル関数が妥当であると言える場合が多いが、最終的には null hypothesis probability に基づいた、つまり帰無仮説を背景にした考察に戻るべきである。なお、が自由度 より非常に小さい場合、つまり が 0 に近 い場合、フィットがよく合っていると考えるよりも、誤差評価に誤りがあるなどを疑うべ きであろう。また、このフィットにおいて保証されるのは「用いたモデル関数がデータを表すものとして適切であったか」であり、「そのモデル関数が物理的に適切かどうか」に は言及していないことに注意しておいた方がよい。フィットを合わせられることと、デー タを正しく物理解釈することは等価ではない。
なお、検出器の応答関数を考慮し、天体から地球に届いた段階での検出器に入る前の スペクトルは、
とすれば表示され、上図の最後のパネルのような形になる。この天体は非熱的な「べき」成分の天体スペクトルが星間吸収を受けてきたものと解釈できる。
次のステップとして、モデル関数パラメータの誤差の範囲 (信頼できる値の範囲) を error コマンドを用いて求めよう。XSPEC の error コマンドでは、fit コマンドで求まった に対して、与えられた を足した 以下のを与えるようなパラメータ値の範囲を計算する。この際、誤差を求めたいパラメータ以外の、値を固定されていないフリーな全てのパラメータをふりの最小化を行う。これは自由度 1 の分布に対応し、XSPEC のデフォルト設定では 90% の信頼区間に対応する の範囲から誤差を求める。星間吸収、「べき」の光子指数、強度はそれぞれパラメータ番号が 1, 2, 3 で、あれば、以下のように90% 信頼区間に対応する誤差を求めることができます。
つまり、以下のように 90% 信頼区間の誤差を得ることができます。
実は、このデータは、以下のパラメータをもつモデルから擬似的に作成したもので、以下が本来の値です。
光子指数と強度は誤差の範囲で正しく得られています。星間吸収は低エネルギー側で決まるので、あまりうまくパラメータを得られていないようです。
ちなみに、90% 信頼区間の誤差はX線天文学でよく使われていたが、より一般的な誤差 (68% 信頼度) の場合には、 として、
とすればよい。
これを確認するために、フリーなパラメータのうち 2 つについて 値の 変化を描くことができる。今回、星間吸収はあまりうまく決まっていなかったので、の値を最適値に固定してみる。つまり、
その上で、
などとすると、2 番目のパラメータ (今回の場合、べきの光子指数) を 2.08 から 2.22 まで 100 ステップ、3 番目のパラメータ (べきの強度) を 9.8 から 12.2 までの 100 ステップにそれぞれ分割し、100 × 100 のグリッド点について他のパラメータを最適化した上で計算する。この結果は
とすると、以下のように表示できる。
この図では、 = 1.00, 2.30, 2.70, 4.61, 9.21 をもつ 5 つの confidence contour を描いてくれる。2 個のパラメータ空間において残り 1 個のパラメータを最適化して を求めるので、この場合は自由度 2 の 分布に従い、3 つの confidence contour はそれぞれ信頼度 68% ()、90% ()、99% () に対応することになる。一方で、XSPEC の標準的な解析で error コマンドを使用する場合には、パラメータ 1 個 ずつについて誤差を計算していので自由度 1 の 分布であり、信頼度 68%、90%での 誤差は で決まる。自由度 2 の 場合よりも誤差範囲はわずかに小さくなる。一般的な解析では error コマンドを使用することが多いが、どうのようなデータ解析を行っているのかには常に注意が必要である。また、実験結果をまとめる際に、 誤差をつけることが多いが、XSPEC では 90% 信頼度 での誤差がデフォルトのことが多い。いずれも定義して使えば問題ないが、注意しよう。
検出器のカウントレートは、ignore コマンドで指定したエネルギー範囲について
とすると表示され、実測値は 101 counts s (1-50 keV)とわかる。モデルの推定値や、合計のカウント数も表示される。
フラックスを求めるときは、
などとする。flux 2.0 10.0
としたときは、単に 2-10 keV のフラックスが求まり、 ergs cm s と求まる。flux 2.0 10.0 err 200, 68.3
などとすると、誤差も求めてくれる。この場合は、1誤差。詳しくは、ヘルプを調べるとよい。
ちなみに、よく使われる換算値は、1 Crab = ergs s cm である。
X線光度は
erg/s
なので、Crab まで 2 kpc の距離を考えると、このX線フラックスから推定されるX線光度は、 ergs/s となる。
余談だが、Crab が銀河中心 (d=8 kpc)にいると、ちょうど ergs/s なので、中性子星のエディントン光度くらいになる。覚えておくと便利なときもある。
最後に XSPEC の中での検出器応答の扱いを簡単に補足しておく。 簡単な物理観測では、測定値に対して検出器の検出効率の補正を施せば、元々の物理量、観測量に簡単に戻すことが可能である。一方でX線観測においては、検出器に光子の全エネルギーを与える光電吸収のほかにも、一部分のエネルギーを与えるコンプトン散乱などの過程もある確率で生じている。そのため、検出器の実際の計数から元の X線スペクトルを一意的に求めることは、一般にはできない。そこで、X 線のモデルスペクトルを仮定し、検出器応答を受けた場合に予想されるカウント数のスペクトルを、実際の観測結果と比較する。そして最小化のアルゴリズムによって、再帰的に最適なモデルパラメータを求めるという方法を採用する。これは、
と表せる。ここで、
である。エネルギー keVの単色X線で予想されるスペクトルは、 であり、エネルギー keVに対応する検出効率はと表せる。 一般にに逆行列が必ずしも見つかるわけではないため、最小化アルゴリズ ムを採用しているとも考えられる。
XSPEC では を RESP (Response File)、 を ARF (Ancillary Response File) として読み込むが、検出器ごとに S をどこが担うかは異なる場合がある。
たとえば、「すざく」の X 線 CCD カメラ (XIS) では、ARF、RESP を両方用いているが、硬 X 線検 出器の PIN 型半導体 (HXD-PIN) では、RESP のみに全てを含めるようになっている。また、硬 X 線検出器の GSO 結晶シンチレータ (HXD-GSO) では、若干の補正として ARF を用いるようにしている。それぞれがどのような関数形をしているかは、FTOOLS の fv コマンドなどで中身を見てみるとよい。
統計の悪いデータを解析する際、ビンまとめという操作を行うことがある。grppha
などのコマンドで行うことができる。
なお、昔はビンまとめすることが多かったが、これはエッジやラインなどのシャープな構造がもっている情報を消してしまうので、C-stat を使う人も増えている。
観測を行う前に、どのようなスペクトルを持つかなどをシミュレーションしたいことも多い。こういった場合、fakeit という xspec のコマンドを使うことができる。
まず、仮定として、0.3 keV の温度をもつ中性子星が 1 kpc の距離にある場合にどう見えるかを考えてみよう。つまり、
まず、xspec でソース・スペクトル(なんでもよい)とバックグラウンド、レスポンスを読み込んだ後、モデル関数を指定する。そして、fakeit コマンドで観測時間を指定する。これで、だいたいの統計誤差の大きさや、どれくらいのエネルギー帯域で検出ができそうかなどを把握することができる。
すると、sim_ninja_nstar_5ks.fak と sim_ninja_nstar_5ks_bkg.fak という2つのファイルができる。これを読み込んで、同じモデルでフィットしてみると、以下のような結果を得る。