astropy を使ってドップラートラッキングをする

背景

はじめに

• 天の川銀河をただよう恒星や星間ガスは、銀河中心の周りを回転しています。

  • かっこいいシミュレーション動画

• 私たちがこういった構造を観測しようとする場合、例えば、電波分光装置を使って星間ガスを観測します。

(Credit: NASA/JPL-Caltech/ESO/R. Hurt, いらすとや)

天体の速度と、観測者の速度

• 観測で得られた輝線は、周波数がズレており、ズレの分だけガスは速度を持っています (ドップラー効果)。

• ただし、ドップラー効果は、天体の速度と、観測者の速度の 2 つが影響します。

  • 天体の速度は 100〜10 万年くらいは変わらないのに対し、
  • 観測者は、地球の自転 (1日周期)、公転 (1年周期) で速度が変わります。

• 天体の速度だけを記録したいので、地球の速度とは関係のない、普遍的な静止系を定義してやると便利です。

• 静止系として、電波天文学で用いられているのは、局所静止基準 (LSR)です。

  • 局所静止基準は、太陽系の中心を原点とした静止系です。
  • ただし、局所静止基準の速度は、太陽系の速度を採用しませんでした。
    • なぜなら、局所静止基準には、太陽系の速度ではなく、天の川銀河の中で、太陽の位置での、銀河としての大局的な速度を採用したいからです。
    • 太陽系の速度は、この場所での銀河としての大局的な速度からは外れており、適しません。
  • 局所静止基準の速度には、太陽系の位置での銀河としての大局的な速度 (太陽系周辺のいくつかの恒星の速度から求める) が使われます。

天体速度を取得する (ドップラートラッキング)

• したがって、「天体の速度を記録する」とは、局所静止基準に対する速度 $V_{LSR}$ を記録する、ということになります。
  • 観測者が地球上で観測した天体信号から、LSR に対する観測者の速度 $V_{observer}$ を引けば $V_{LSR}$ が得られます。
    • $V_{observer}$ は、
      • LSRに対する太陽速度 : $V_{sun}$
      • 太陽に対する地球の公転速度 : $V_{earth}$
      • 地球に対する観測者の自転速度 : $V_{location}$
      • からなります。
  • $V_{observer}$ を引く処理のことを、ドップラートラッキングと呼びます。

太陽速度が変わってしまう...！

• LSR に対する太陽速度 $V_{sun}$ は、観測精度が向上したり、LSR の定義の違いによって、変わってしまいます。

  • 例えば、最近の研究 Schönrich et al. (2010) によれば、$(U, V, W) = (11.1, 12.24, 7.25)$ ${\rm km,s}^{-1}$ とされました。
    • $U$ : 太陽速度のうち、銀河中心方向への成分 (銀経 = $0^\circ$ の方向)
    • $V$ : 太陽速度のうち、銀河面に沿って、銀河中心へ垂直な成分 (銀経 = $90^\circ$ の方向)
    • $W$ : 太陽速度のうち、銀河面に垂直な成分

• これでは、採用する結果が異なれば、同じ天体の観測であっても $V_{LSR}$ が変わってしまって不便です。

• そこで、電波天文学における $V_{sun}$ は次のように定義することになっています。

  • 太陽は局所静止系に対して1900年分点で赤経$18{\rm h}$、赤緯$+30^{\circ}$の方向に$+20,{\rm km,s}^{-1}$で運動している (天文学辞典より引用)
  • この時、$(U, V, W) = (10.27, 15.32, 7.74)$ ${\rm km,s}^{-1}$ です。

astropy で計算してみよう！

それでは、astropy を使って、実際に $V_{observer}$ を計算し、ドップラートラッキングをして見ましょう。

太陽速度 Vsun = (U, V, W) を計算する

import astropy.coordinates
from astropy.units import km, s, deg
from numpy import cos, sin

# 太陽の運動方向を定義
dir_sun = astropy.coordinates.SkyCoord(
    ra = 18 * 15 * deg,                   # R.A. = 18 h
    dec = 30 * deg,                       # Dec. = 30 deg
    frame = 'fk4',                        
    equinox = astropy.time.Time('B1900'), # 1900 年分点を指定する
).galactic

# 太陽の速さを定義
v_sun = 20 * km / s

# 太陽速度を、銀河面に直交した座標系に成分分解し、U, V, W を計算
U = v_sun * cos(dir_sun.b) * cos(dir_sun.l)
V = v_sun * cos(dir_sun.b) * sin(dir_sun.l)
W = v_sun * sin(dir_sun.b)

# CartesianDifferential 型にしておく
v_bary = astropy.coordinates.CartesianDifferential(U, V, W)

# >>> v_bary
# <CartesianDifferential (d_x, d_y, d_z) in km / s
#     (10.27059164, 15.31741091, 7.73898381)>

• これで、$V_{sun}$ (変数名は v_bary) が計算できました。
• とても簡単ですね。

Vobserver を計算する

import time
import astropy.coordinates
import astropy.time
from astropy.units import deg, m
from numpy import cos, sin

# 観測者の地球上の座標を設定します。例は長野県野辺山です。
loc1p85 = astropy.coordinates.EarthLocation(
        lon = 138.472153 * deg,
        lat = 35.940874 * deg,
        height = 1386 * m,
)

# 観測時刻を定義します
tobs = astropy.time.Time(time.time(), format='unix')

# 観測方向を定義します。例は Ori-KL です。
target = astropy.coordinates.SkyCoord(
    ra = 83.809 * deg,
    dec = -5.372639 * deg,
    frame = 'icrs',
    obstime = tobs,     # 時刻を設定する
    location = loc1p85, # 観測者座標も設定する
)

# V_observer (LSR 系からみた観測者の速度) の計算
v_observer = astropy.coordinates.SkyCoord(
    # 観測者位置を、観測時刻での GCRS 座標系に変換する
    loc1p85.get_gcrs(target.obstime)
).transform_to(
    # それを、LSR 座標系に変換する
    astropy.coordinates.LSR(v_bary=v_bary) # 太陽速度を指定する
                                           # デフォルトでは、Schönrich et al. (2010) が使われてしまう
).velocity

# V_obs のうち、天体方向の成分を取り出す
v_correction = v_observer.d_x * cos(target.icrs.dec) * cos(target.icrs.ra) + \
               v_observer.d_y * cos(target.icrs.dec) * sin(target.icrs.ra) + \
               v_observer.d_z * sin(target.icrs.dec)

# tobs, ra, dec には配列も渡せます

• あらかじめ計算しておいた $V_{sun}$ (変数名 v_bary) を用いて、$V_{observer}$ を計算し、さらに、$V_{observer}$ のうち、観測天体の方向の成分を v_correction として計算しました。
• これで、観測された速度から v_correction を引けば $V_{LSR}$ を取得できます。
• コードの大半は座標の定義です。計算はとても簡単ですね。astropy 最高。

メモ

• WARNING: failed to download http://maia.usno.navy.mil/ser7/finals2000A.all が出る場合
  • IERS のテーブルの保存元が、アメリカ海軍天文台から、NASAに変わった
    • 新しいダウンロード先リンク : ftp://cddis.gsfc.nasa.gov/pub/products/iers/finals2000A.all
  • 新しい astropy では対応されているので、アップデートしましょう