7. TIMING ANALYSIS; XRONOS

7.1 Overview of XRONOS

ライトカーブを描いたり、時系列解析を行なう場合は、HEADAS に含まれる XANADU package の XRONOS ツール群を使用する。「あすか」の頃は、X線解析の ための独立したツール群となっていて、XRONOS 環境に入ってから使う仕様になっ ていたが、今は、個別のftoolsとして、コマンドラインで使えるようになってる。 詳細は、HEASARC の大元のページである、

http://heasarc.gsfc.nasa.gov/docs/xanadu/xronos/xronos.html

を参考のこと。記述に不一致があれば、HEASARC の方を信じ、 suzaku_local@crab.riken.jp に連絡ねがいます。大まかには下記のようなこと ができる。

• ライトカーブの描画 (lcurve)
• パワースペクトルの作成 (powspec)
• 周期的変動の探査 (efsearch)
• 畳み込みライトカーブの作成 (efold)

7.2 Barycentric Correction for Suzaku

衛星は地球のまわりをまわっており、地球は太陽のまわりをまわっている。 光は地球半径分すすむのに20 msecかかり、太陽まで8分かかる。 他の衛星と時刻を合わせたい場合は、 衛星の軌道運動による効果を打ち消すために、時刻を地球中心の値に変換し (Geocentric correction)、 さらに、観測時期による違い(地球の公転運動の効果)を打ち消すために、 太陽中心での時刻に補正する必要がある。これを Heliocentric correctionと呼ぶ。 周期の速いX線パルサーや電波パルサー、ミリ秒パルサーの Pulse の到着時刻について、 他の時期の観測と比較するためには、さらに、惑星の運動をすべて加味した太陽系重心への変換が必要となる。 これをbarycentric correctionと呼ぶ。

XRONOS package にも barycentric correction tool は含まれているが、 「すざく」の軌道要素が読めるツールとして aebarycen を HEADAS v6.1.x から リリースするに至った^7.1。 ただし、まだ、絶対時刻にバグがある可能性が高いので全面的には信用しないと共に、 なにかバグをみつけたら suzaku_help@crab.riken.jp まで報告してほしい。 用意するファイルは、時刻補正したいイベントファイル(aeXXXXX.evt)と 軌道ファイル(aeXXXXXX.orb)である。

aebarycen version 2006-07-24
Written by Y.Terada (RIKEN), T.Enoto (UT), Y.ISHISAKI (TMU)
Input/output file or @filelist to be corrected[] ae101005040hxd_0_wel_uf2.evt<--入力File
Orbit ephemeris file for Suzaku[]ae101005040.orb                            <--軌道要素
leapsec file name [] : /usr/local/lheasoft/6.06/headas/ftools/refdata/leapsec.fits 
                      <--- 閏秒のファイル( HEADAS をインストールした場所にある)
                           (CALDB と書けば CALDB から適切に取って来てくれる)
Right Ascension (HeaderKey/NNhNNmNNs/deg): 500 <--- 天体のRA
Declination (HeaderKey/+NNdNNmNNs/deg): 50    <--- 天体のDEC 
        (デフォルト値は用意されていない。自分の興味のある天体の座標をいれる。)
        (このように定義外の値を与えると, event FITS のヘッダの RA_NON, DEC_NOM が使われる)
....

こうして出来た FITS file を見てみると、

HISTORY -----------------------------------------------------
HISTORY  aebarycen version 2006-07-24 at 2006-08-02T07:34:53
HISTORY -----------------------------------------------------
HISTORY   filelist='ae101005040hxd_0_wel_uf2.evt'
HISTORY   orbit='ae101005040.orb'
HISTORY   leapfile='/usr/local/lheasoft/6.06/headas/ftools/refdata/leapsec.fits'
HISTORY   (RA,DEC) = (500,500) = ( 300.0000 , 40.0000 )
HISTORY   time_col='TIME', start_col='START', stop_col='STOP'
PLEPHEM = 'JPL-DE200'          / Solar System ephemeris used for baryctr corr.

といった履歴がかかれ、

TIMESYS = 'TDB     '           / Coodinate Reference System
TIMEREF = 'SOLARSYSTEM'        / Times are pathlength-corrected to barycenter
TSTART  =     206432797.694330 / time start   (+243.395389 s barycen corrected)
TSTOP   =     206487866.716041 / time stop    (+244.418146 s barycen corrected)
DATE    = '2006-08-02T07:34:53' / file creation date (YYYY-MM-DDThh:mm:ss UT)

のように TSTARTやTSTOPの値などが書き換わり、TIMEREF が、'LOCAL' から 'SOLARSYSTEM'に 変わっているはずである。もちろん、event 本体の時刻や GTI も変更されているはずである。 確認のため、補正前の時刻との差で最もよく效く値として、 Hericentric correction の値を計算したい場合、天体の座標(RA,DEC)と地球の位置とから、 天体へのベクトルを太陽公転面の垂線に射影した距離から図7.1を参考にするとよい。

図 7.1: 太陽中心補正で変更される時間差

7.3 Make Light Curve

下記の要領でライトカーブが描ける。

% lcurve 
 lcurve 1.0 (xronos5.21)
 
Number of time series for this task[1] 1
                    <--------- (1 つのファイルだけをプロットするなら1, 2つ以上ならその数)
Ser. 1 filename +options (or @file of filenames +options)[] crabhxdnor_hxd_wel_pin.evt 
                    <--------- (Event FITS もしくは light curve FITSの名をいれる)
.....
 Default Newbin Time is:   39.632403    (s) (to have 1 Intv. of     512 Newbins)
 Type INDEF to accept the default value

Newbin Time or negative rebinning[100] 100.0
                     <------------ ライトカーブの一ビンの時間幅を入力。今の場合100秒。
                                   当然、最小の時間分解能より小さくすることは出来ない。
 Newbin Time ......     100.00000      (s)
 Maximum Newbin No.               203
 Default Newbins per Interval are:         203
 (giving       1 Interval of          203 Newbins)
                     <----------- イベントのカバーする時間範囲と、先程入力した時間ビンから、
                                  トータルで何ビンのライトカーブになるかを教えてくれている。

Number of Newbins/Interval[239] 203
                     <----------- 描画する時間領域(Interval) あたり何ビンを表示するか入力。
                                  今は、一画面でみたいので、先に教えてもらった 203 を入力。
                                  たとえば、20 を指定すると、21枚分 描画される。

Maximum of       1 Intvs. with          203 Newbins of       100.000     (s)
                     <----------- 単なる確認。1 Interval となっていますね。
Name of output file[default]   <------------- このあたりはdefault 、リターンでOK
Do you want to plot your results?[yes] <----- リターン
Enter PGPLOT device[/XW]  <-----------------  リターン
....
PLT> q    <------------------------ QDP に入る。自由にどうぞ。

複数のファイルを 1 シリーズとしてまとめる時は、ファイルをリストアップした flist というファイルを作り、 Ser.1. ... [] の答えを @flist とすればOK。また、2 シリーズ以上を与えた時は、最後にプロットする段階で、 複数プロットモードか、比を取るモードか、などを聞かれるので、素直に答えてやればよい。

7.4 Make Power Spectrum

周期的な変動がある場合は、 ライトカーブの各周期毎のパワーを描いた「パワースペクトル」を描いたときに、 その特定の周波数に高いパワーがあらわれる。 「すざく」の場合は、SAA や COR、地球の蝕などによって、 データが離散的であるので、人為的にでた「嘘のパワー」を見極め、 天体固有の周期変動を探査する必要がある。

% powspec 
Ser. 1 filename +options (or @file of filenames +options)[] crab_hxd_wel_pin.evt
                    <--------- (Event FITS もしくは light curve FITSの名をいれる)
...
Name of the window file ('-' for default window)[] -
                    <--------- デフォルトでよい。リターン
...
Default Newbin Time is:   2.4743037    (s) (to have 1 Intv. of    8192 Newbins)
                    <--------- 時間分解能の適当な値を教えてくれるが、参考程度にするだけ。
Newbin Time or negative rebinning[0.1] 0.01
                    <--------- パワースペクトルを描くための、元となるライトカーブの時間分解
                               能を入力。ミリ秒を探査したいのに 大きな値 10.0 sec として
                               もダメだし、500 sec程度の周期変動を探査したいのに、あまり
                               細かな時間ビンを指定しても、一ビンあたりの統計が悪く、やたら
                               ビン数も増えるばかりでよろしくない。
 Newbin Time ......    0.10000000E-01  (s)
 Maximum Newbin No.           2026802
                    <--------- 指定した時間分解能だと、全体で何ビンのライトカーブかを教えて
                               くれる。(参考値) 



つぎに、「計算の元となるライトカーブ」を何枚分用意するか、を考える。
 Default Newbins per Interval are:        8192
 (giving     248 Intervals of         8192 Newbins each)
                    <--------- 参考値。8192 bin づつのライトカーブなら 248枚できる、と
                               言っている。
Number of Newbins/Interval[8192] 204800
                    <--------- 一枚(1 Interval)あたりのビン数を入力。ここでは、約10枚
                               (10 Interval)を予定して、(8192x248/10) 程度。

 **** Warning: No. of Newbins/Intv. reset to   262144
 Maximum of       8 Intvs. with       262144 Newbins of      0.100000E-01 (s)
Default intervals per frame are:         8
                    <--------- 適切な調整がなされ、だいたい予定どおり、8 Interval の計
                               算となった。




最後に、Interval 毎のパワーの計算を何枚(Frame)づつ重ねて描画するか、を決める。
Number of Intervals/Frame[25] 8
                    <--------- 8 Interval 全部を 1枚(1 Frame) のパワースペクトルにま
                               とめたいので、ここでは 8 を指定。たとえば 2 とすると 
                               8/2 = 4 Frame 結果がでる。

Rebin results? (>1 const rebin, <-1 geom. rebin, 0 none)[0] <------ デフォルト
Name of output file[default]  <------------------------------------ デフォルト
Do you want to plot your results?[yes]  <-------------------------- デフォルト
Enter PGPLOT device[/XW] <----------------------------------------- デフォルト
...
PLT> q    <------------------------ QDP に入る。自由にどうぞ。

こんな感じでやると、図7.2のようなパワースペクト ルがかける。

図 7.2: Crab パルサーのパワースペクトル

7.5 Perform Periodic Search

前の節で、およその周期が分かったので、その精度をあげる。 ベストの周期で、天体のライトカーブを畳み込んでやると、 ひじょうに変動の大きなライトカーブになるはずである。 逆に、無関係な周期で畳み込むと、平らなライトカーブとなるであろう。 サーチプログラムでは、 次節7.6で述べるような畳み込みのライトカーブを、 いろいろな周期で描き、それが一定なのか否かを統計的に$\chi^2$検定し、 その$\chi^2$が最大なるような周期を算出する。 下記の例では、Crab パルサーに対し、 33.5810 msec を中心に 0.0002 msec の分解能で +/- 128 個分のトライアルを行なってみる。

% efsearch 
Ser. 1 filename +options (or @file of filenames +options)[file1] evt/crab_hxd_wel_pin.evt 
                    <--------- (Event FITS もしくは light curve FITSの名をいれる)
... 
Name of the window file ('-' for default window)[-] 
                    <--------- デフォルトでよい。リターン
...

計算の元となる 畳み込みライトカーブの時刻原点と、探査したい周期を与える
 Default Epoch is:  13628.00000
 Epoch format is days.
Epoch[34 234.23] 13628  <------------- 原点: MJD から 40000 引いた値を用いる(*注)
 Period format is seconds.
Period[88.87] 0.0335810 <------------- 周期: default設定 なら 秒単位で与える
 Expected Cycles ..         603555.91 <---- 何周期分のデータに相当するか表示(参考値)
...

畳み込みライトカーブの 0 〜 1 phase を何ビンに分けるかを設定
 Default phase bins per period are:         8      <--- 参考値
Phasebins/Period {value or neg. power of 2}[-3] 32 <--- 自分の設定値
...

次に、ひとつの計算で用いる時間範囲を指定する
 (giving       1 Interval of     19313790 Newbins) <----- 全体を1 Interval にす
 Type INDEF to accept the default value                   るならこの値(参考値)
Number of Newbins/Interval[10] 19313790  <--------------- 設定値
 Maximum of       1 Intvs. with     19313790 Newbins of      0.104941E-02 (s)
...

最後に探査する時間分解能とビン数を指定する。
Resolution for period search {value or neg. power of 2}[.03] 0.0000002
                                  <------------ 0.0002 ミリ秒ステップで計算
Number of periods to search[100] 128 <--------- 前後 128 x 0.0002 msec。
...
Name of output file[default] <----------------- デフォルトでOK
Do you want to plot your results?[yes] <------- デフォルトでOK
Enter PGPLOT device[/XW] <--------------------- デフォルトでOK
PLT> q    <------------------------ QDP に入る。自由にどうぞ。

結果は下記のような感じででる。これを繰りかえし、時間分解能をあげて探査するとよい。

図 7.3: Crab パルサーの周期探査の結果画面

周期的変動をサーチする際、検出器の時間分解能以上の時間が探査できることに疑問を持つかもしれない。 たとえば、1.00 秒しか時間分解能のない検出器があったとしても、 100秒の間に 100 pulse が含まれるか、99 pulse しかふくまれないか、が区別できれば、 入射イベントが 1.00 秒周期なのか 1.01 秒周期なのかの区別がつく。 すなわち、周期的な変動の探査には、検出器の時間分解能だけでなく、 観測時間が何サイクル分カバーしているかも効いてくる。

(*注) XRONOSのデフォルトの時間軸の単位は TJD (Truncated Julian Days TJD=JD-2440000.5) である。一方、MJD (Modified Julian Days) = JD-2400000.5 [day]) なので、 TJD = MJD - 40000 [day] となり、40000 days のずれがでる。

7.6 Make Folded Light Curve

周期が求まったので、次に、その周期のなかで天体がどのような波形で変動しているのかを知るために、 この時間で畳み込んだライトカーブを描いてみる。 周期 $P$で変動すると考える天体からある時刻 $T$ にイベントが到来し、 その時刻$T$ が、ある時刻原点$T_0$ から数えて$N$(整数) + $x$(0..1) サイクル目に 相当する場合( $T = T_0 + P \times (N+x)$)、 x は phase と呼ばれ、0 から 1 の値をとる。 ふつうのライトカーブは横軸 $T$、縦軸 カウントのヒストグラムであるが、 横軸 $x$ 縦軸カウントとしたライトカーブは、 畳み込みライトカーブと呼ぶ。(周期 $P$で折り畳んで表示しているようなもの) 解析では下記のように efold というツールを用いて描画できる。

% efold
Number of time series for this task[1] 1
                    <--------- (1 つのファイルだけをプロットするなら1, 2つ以上ならその数)
Ser. 1 filename +options (or @file of filenames +options)[] crabhxdnor_hxd_wel_pin.evt 
                    <--------- (Event FITS もしくは light curve FITSの名をいれる)
... 
Name of the window file ('-' for default window)[-] 
                    <--------- デフォルトでよい。リターン
...

 Default Epoch is:  13628.00000
 Epoch format is days.
Epoch[] 13628 <----------------------- 畳み込みの原点: Phase = 0 の値: MJD - 40000
 Period format is seconds.
Period[] 33.5809E-3    <-------------- 畳み込み周期: default設定 なら 秒単位で与える
Period derivative       [0.00] <------ P_dot (defaultなら聞かれない)

Phasebins/Period {value or neg. power of 2}[17] 257 
                       <--- 畳み込みライトカーブの 0 〜 1 phase を何ビンに分けるかを設定
...
 
 Default Newbins per Interval are:   155114332
 (giving       1 Interval of    155114332 Newbins)
                       <---- 計算を すべてまとめて 1枚(1 Interval) でおさめる場合の
                             ビン数が表示される(参考値)
 Number of Newbins/Interval[10260482] 155114332
                       <----- 計算は、1 Interval で行ないたいので、参考値そのまま。
Default intervals per frame are:         1 <---- 描画の枚数(Frame)の参考値
                                                 今は 1 Intervalしかないので1枚。
Number of Intervals/Frame[1] 1 <------------ 1 枚にまとめて表示と入力

...
Name of output file[]  <-------------------- デフォルトでOK
Do you want to plot your results?[yes] <---- デフォルトでOK
Enter PGPLOT device[/xw] <------------------ デフォルトでOK
PLT> q    <------------------------ QDP に入る。自由にどうぞ。

結果は、次の図7.4に示すようなかたちで表示される。 一番はじめの箇所で、複数のシリーズを指定した場合は、画面が縦に分割されて表示される。

図 7.4: Crab パルサーのライトカーブを 33.5809 ミリ秒で畳み込んだライトカーブ

7.7 Additional Comments

次の段階としては、図7.4の phase で分割してスペクトルを 描きたい、などの要求がでてくるだろう。詳細は xselect の章 ($\S$5)にまかせる が、

xselect 上で filterする際、
xsel:SUZAKU-HXD-WELL_PIN > filter PHASE
> Enter what filter? >[TIME] PHASE
> Enter the epoch for phase selection ( in MJD ) >[] 53628.0
> Enter the period for phase selection ( in DAYS ) >[] 3.88668e-07
> Enter phases from 0-1 (e.g. 0.1-0.25,0.3-0.45) >[] 0-0.5

といった感じでイベントを抜きだせる。なお、時刻原点は XRONOS と MJD では 値が40000ずれて いること(§7.5注)、 Period が秒ではなく Day であることに注意する。 また、Period があまりに短いものは xselect 上では切り出せないことがあるようである。 (bug でしょう) 一度 filter したものを event におとし、efold などで正しく切り出されているかを 確認して次のステップにすすむのがよかろう。