はじめに - CMB とは -

現代の標準宇宙論によれば、ビッグバンによって宇宙が誕生した後、加速膨張により急激に温度が下がった結果、（ビッグバンから約40万年後に）電子と陽子が再結合することにより光が宇宙を直進できるようになった頃の初期宇宙の放射が存在すると考えられています。これを「宇宙の晴れ上がり」と呼びます。この放射場は、宇宙が膨張するにつれて赤方偏移を受けて波長が伸び、現在（138億年）では、マイクロ波として宇宙至る所に存在していると考えられています。これを宇宙マイクロ波背景放射（Cosmic Microwave Background / 通称 CMB）と呼びます。

CMB は1965年に A.Penzias と R.Wilson が高感度マイクロ波アンテナの設置をしていた際、（説明のつかない雑音として）偶然発見されました。CMB はビッグバン宇宙論の観測的証拠であるだけでなく、我々の宇宙の歴史や物理的性質や組成を知る上で非常に重要な観測データとなります。２人は後にこの発見によってノーベル物理学賞を受賞しています。

CMB の著しい特徴として、非常に等方的であることが挙げられます。宇宙のありとあらゆる方向から同じように観測され、そのスペクトルは約3Kの非常に低い温度の黒体放射としてよく近似できます。一方でわずかながら非等方性を持っており、この非等方性こそが初期宇宙の密度揺らぎを反映していると考えられています。このわずかな密度揺らぎから重力不安定によって銀河団・銀河が生まれ、星が生まれ、惑星が生まれたのですから、CMBは現在の宇宙への進化を記述する最も古い地図とも見なせるわけです。しかしながら、この非等方性は天の川銀河中の星間ガスやダストのシグナルに比べて非常に小さいため、高精度観測とともに天の川銀河からのコンタミネーションをいかに上手く除去できるかが CMB観測の鍵となります。

CMB 観測プロジェクト

CMB を精度良く観測するために、これまで以下の３つの衛星ミッションが NASA および ESA によって行われました。

COBE (1989-1996), NASA
CMBの非等方性が初めて詳細に観測されたプロジェクト
WMAP (2001-2010), NASA
より高精度・高解像度で対局的な CMB 揺らぎの測定を目指したプロジェクト
PLANCK (2009-2013), ESA/NASA
WMAPよりさらに小さいスケールで局所的に CMBを測定したプロジェクト

下図は、それぞれの観測から得られた宇宙全天の温度揺らぎ分布を示しています。時代が進むにつれてより高解像度の測定が実現できていることがわかると思います。

Wilkinson Microwave Anisotropy Probe: WMAP 衛星観測データ

では、実際に観測された CMBデータをいじって遊んでみましょう。今回は、解像度とデータ容量を考慮して中間の WMAP データを用います。データは以下のサイトで公開されており、研究・教育（・遊び）目的で好きにダウンロードして使うことができます。
https://lambda.gsfc.nasa.gov/product/

WMAP は全天における温度場分布を、以下の表のような複数の異なる周波数帯で観測しています。これらの周波数帯で観測されるデータを総合的に用いることによって、天の川銀河等からのコンタミネーションを上手く除去してゆきます。

バンド	K	Ka	Q	V	W
周波数 [GHz]	23	33	41	61	94

それでは、実際に以下のコマンドで各周波数帯の温度擾乱データをまとめてダウンロードしてみましょう。

wget https://lambda.gsfc.nasa.gov/data/map/dr5/skymaps/9yr/raw/wmap_band_imap_r9_nineyear_v5.tar.gz

それぞれのデータが25MBで５バンド分あるので、全部で126MBとなるので注意。ダウンロードしたファイルは展開して、ディレクトリ名を data などと変更しておきます。

準備: Healpy のインストール

全てのCMB データはHEALPix という全天表示に特化した特殊な座標システム（球座標は空間非一様なので等方的なCMB解析には適さない）で表されているため、まずはこの座標システム上で様々な計算・描画ができるpython パッケージをインストールします。以下のコマンドで簡単にインストールできます。

conda install -c conda-forge healpy

実際にプロット: 各周波数帯での温度擾乱分布の違い

それでは、実際に以下の python スクリプトを用いて、試しに K-バンドとW-バンドのデータをプロットしてみます。

plot_KW.py

import healpy as hp
import matplotlib.pyplot as plt 

#read the K- W-band maps
wmap_K = hp.fitsfunc.read_map('data/wmap_band_imap_r9_9yr_K_v5.fits')
wmap_W = hp.fitsfunc.read_map('data/wmap_band_imap_r9_9yr_W_v5.fits')

#plot the data
fig = plt.figure(figsize=(13,4))
plt.rcParams["font.size"] = 13

hp.mollview(wmap_K, cmap='turbo', min=0, max=2, title='K-band', sub=[1,2,1], cbar=False)
hp.mollview(wmap_W, cmap='turbo', min=0, max=2, title='W-band', sub=[1,2,2], cbar=False)

plt.show()

出力は以下の通り。ただし、共にカラーマップは0-2 mK で表示してあります。

おぉ。だいぶ様子が違いますね。
まず、K-バンドでは、赤道面（黄道面）を走る天の川銀河からの放射が非常に強く見えています。これは、天の川銀河に大量に存在する星間物質内で強く発生するシンクロトロン放射・制動放射が原因と考えられます。一方、W-バンドでは、星間物質由来のコンタミネーションは抑えられているものの、今後は逆にダスト由来の熱放射が混じってきてしまっています。
このように、他の Ka, Q, V-バンドでもそれぞれ特定の要因のコンタミネーションに感度を持っており、それぞれ微妙に異なる全天温度分布が得られます。そこで、これらのデータを上手く線形結合させることで、天の川銀河からのノイズ除去を試みます。

内部線形結合（Internal Linear Combination）: ILC法

ここでは最も単純な内部線形結合法（通称 ILC法）によるノイズ除去を試してみます。
まず各バンドにおける温度擾乱分布 $\Delta T_{\nu}$（ただし $\nu=$K,Ka,Q,V,W）はバンドによらない CMB本体からの放射と、周波数依存のノイズからなるとします。ILC重み係数を $w_{\nu}$として、最終的に重ね合わされたILCデータの全分散（total variation / TV）が最小となるように係数を決定します。この際、$\Sigma_{\nu}w_{\nu}=1$ を制約条件として、ラグランジュの未定条数法を用いることで$w_{\nu}$が以下のように求まります。より詳しいことを知りたい方は、Eriksen et al. (2004) などを参照。

バンド	K	Ka	Q	V	W
ILC 係数 $w_{\nu}$	0.028	-0.283	-0.370	1.666	-0.042

これを見ると、ILCマップへの大部分の寄与はV-バンドから来ていて、それの余分を主にKa-バンドとQ-バンドによって相殺していており、K-バンドとW-バンドはごくわずかな寄与に留まっていることが分かります。

それでは、実際にこれらの算出された係数を用いて各周波数マップを重ね合わせてILCマップを構築してみましょう。以下のpythonスクリプトを用います。ただし、線形結合するにあたって、各マップの空間解像度を同じにする必要があるため、hp.sphtfunc.smoothing関数を用いてガウス平均化しました。ここでは、簡単のため、平均化パラメータは既知としましたが、実際には各周波数におけるビーム関数から算出する必要があります（より詳しく知りたい方はこちらからビーム関数をダウンロード可能）。

plot_ILC.py

#read data
wmap_K = hp.fitsfunc.read_map('data/wmap_band_imap_r9_9yr_K_v5.fits')
wmap_Ka = hp.fitsfunc.read_map('data/wmap_band_imap_r9_9yr_Ka_v5.fits')
wmap_Q = hp.fitsfunc.read_map('data/wmap_band_imap_r9_9yr_Q_v5.fits')
wmap_V = hp.fitsfunc.read_map('data/wmap_band_imap_r9_9yr_V_v5.fits')
wmap_W = hp.fitsfunc.read_map('data/wmap_band_imap_r9_9yr_W_v5.fits')

#smoothing factors
sigma_K = 0.0
sigma_Ka = 0.23273851815047383
sigma_Q = 0.28770995290055545
sigma_V = 0.3291958289037056
sigma_W = 0.34569507986473813

#gaussian smoothing
wmap_K_smooth = hp.sphtfunc.smoothing(wmap_K, sigma=sigma_K/180.0*np.pi)
wmap_Ka_smooth = hp.sphtfunc.smoothing(wmap_Ka, sigma=sigma_Ka/180.0*np.pi)
wmap_Q_smooth = hp.sphtfunc.smoothing(wmap_Q, sigma=sigma_Q/180.0*np.pi)
wmap_V_smooth = hp.sphtfunc.smoothing(wmap_V, sigma=sigma_V/180.0*np.pi)
wmap_W_smooth = hp.sphtfunc.smoothing(wmap_W, sigma=sigma_W/180.0*np.pi)

#computed ILC coefficients
w_K = 0.02797237
w_Ka = -0.28258839
w_Q = -0.36961155
w_V = 1.66602016
w_W = -0.04179258

#construct the ILC map
wmap_ILC = w_K*wmap_K_smooth+w_Ka*wmap_Ka_smooth+w_Q*wmap_Q_smooth+w_V*wmap_V_smooth+w_W*wmap_W_smooth

#Official result
wmap_official = hp.fitsfunc.read_map('data/wmap_ilc_9yr_v5.fits')

#---- fig ----#
fig = plt.figure(figsize=(13,4))
plt.rcParams["font.size"] = 13

hp.mollview(wmap_ILC, cmap='turbo', min=-0.15, max=0.15, title='Our ILC map', sub=[1,2,1], cbar=False)
hp.mollview(wmap_official, cmap='turbo', min=-0.15, max=0.15, title='Official ILC map', sub=[1,2,2], cbar=False)

plt.show()

こうして得られたILCマップを下図（左）に示しました。なお、カラーマップは $\pm 0.15$ mKとしてあります。背景放射の温度が約3Kであったことを考えると、擾乱はざっくり約$10^{-4}$の大きさということになります（ただしスケール依存であることに注意）。一応比較のため、NASA が出している「公式の」ILCマップも下図（右）に示しました。これも以下のリンクからダウンロードすることができます。

wget https://lambda.gsfc.nasa.gov/data/map/dr5/dfp/ilc/wmap_ilc_9yr_v5.fits

どうでしょう、なかなかいい感じなのではないでしょうか。
そもそもILC法は場当たり的で「物理的背景を持たない」最も単純な線形ノイズ除去手法であったことを思い出せば、たった５つの周波数帯データを用いただけでここまで綺麗に天の川銀河由来のノイズが除去できてしまうのですから、驚くべきことだと思います。やはり公式データに比べると、黄道面に残っているノイズが目立ち少し見劣りしますが、公式のILCではより洗練された手法・より複数の波長帯データを使用しているので、仕方ありません。

CMBのパワースペクトル

それでは、得られたILCマップを球面調和関数展開して、パワースペクトルを計算してみましょう。
実は、ノイズの影響を極力減らすために既知の光源やノイズ源を全天マップから上手く隠してくれるマスクデータなるものも公式に提供されています。せっかくなので、これもダウンロードして、マスクの効果も調べてみましょう。

wget https://lambda.gsfc.nasa.gov/data/map/dr5/ancillary/masks/wmap_temperature_kq85_analysis_mask_r9_9yr_v5.fits

マスクの概形は下図の通りです。マスクを噛ませることで、ILCマップの中でも特にクリーンな箇所のみに注目できるようになります。

パワースペクトルの計算には、hp.sphtfunc.anafast関数を用います。ここでは、CMBの温度擾乱も等方的に分布していると仮定し、球面調和スペクトルを位数 $m$ について足し合わせ次数 $l$ のみの１次元関数として計算しています。

plot_powerspec.py

#read data
wmap_ILC = hp.fitsfunc.read_map('data/wmap_ilc_9yr_v5.fits')
mask = hp.fitsfunc.read_map('data/wmap_temperature_kq85_analysis_mask_r9_9yr_v5.fits')

#compute power spectra
lmax = 2000
l = np.arange(lmax+1)

#without mask
Cl = hp.sphtfunc.anafast(wmap_ILC,lmax=lmax)
Pl = l*(l+1)*Cl/(2.0*np.pi)

#with mask
Cl_msk = hp.sphtfunc.anafast(wmap_ILC*mask,lmax=lmax)
Pl_msk = l*(l+1)*Cl_msk/(2.0*np.pi)

#---- fig ----#
fig = plt.figure(figsize=(6.5,5))
ax = fig.add_subplot(111)
ax.plot(l,Pl,'-r',linewidth=2,label='without mask')
ax.plot(l,Pl_msk,'-b',linewidth=2,label='with mask')
ax.set_xlim(0,1500)
ax.set_ylim(6.0e-9,4.0e-3)
ax.set_yscale('log')
ax.set_xlabel(r'Spherical harmonic degree $l$')
ax.set_ylabel(r'Power [mK$^{2}$]')
ax.legend()
plt.show()

こうして求めたパワースペクトルが下図になります。赤・青はそれぞれマスク無し・ありに対応しています。$l$ が小さいほど大きな空間スケール、$l$ が大きいほど小さな空間スケールを表しています。

マスク無しの場合だと、$l\approx 120$付近の第一ピークしか再現されていません。なるほど、いくらILC法で天の川銀河からのコンタミネーションを頑張って除去してもやはりノイズの影響は凄まじいらしいです。
一方で、マスク有りの場合だと、$l\approx 550$ 付近の第２ピーク、$l\approx 830$ 付近の第３ピーク、$l\approx 1130$ 付近の第４ピークまでちゃんと分解できていることが確かめられます。公式のパワースペクトルと比較してみても、定性的にはピークを良く再現できてると言えると思います。ILC法恐るべし。

宇宙論パラメータの決定

CMB 観測では、この温度擾乱パワースペクトルを正確に計算することが、非常に重要な手続きとなっています。というのも、このスペクトル形状は理論上様々な宇宙論パラメータに敏感であることが知られており、スペクトルを基に未知の宇宙論パラメータ決定することが可能になるからです。
まず、上で見たようなスペクトルの $l$ が比較的大きい箇所の複数の山や谷はビッグバン初期宇宙を伝播していた音波のスケールを表していると考えられています。音波の情報はその流体の密度を反映しているため、初期宇宙の水素・ヘリウムといった組成に関する情報（宇宙論の言葉でいえばバリオン密度）が得られます。また、より大スケール（小さい $l$）第一ピークのパワー形状からは初期宇宙の曲率を求めることができます。

これらのCMB の詳細な解析や理論モデルとの比較の結果、現在では「標準宇宙モデル」と呼ばれるものが確立しています。このモデルによれば宇宙はほぼ平坦で、我々が良く知る元素（バリオン）は全体の数%しか存在せず、残る大部分は暗黒物質と暗黒エネルギーが占めていると考えられているそうです（これ以上は専門外なので立ち入りません）。

まとめ

いかがだったでしょうか。今回は NASAが公開している WMAP衛星による CMBデータをダウンロードしてきて、ILC法によるノイズ除去・さらにはパワースペクトル計算をして遊んでみました。

「宇宙論」「ビッグバン」と聞くと、ロマンがあって誰しも一度は興味を持ったことがあると思います。しかしながら、ちゃんと理解しようとすると大学・大学院で物理を本格的に勉強しなければならず、ハードルが高いと思います。今回の記事は、そんな状況を少しでも橋渡しできればと思い、執筆しました。

実は、今の時代、誰でも NASAのWebページから研究者が使うような生の観測データがダウンロードでき、手持ちのノートパソコンで手軽にデータ解析が行え、物理学史に残るような重要な研究を再現しながらその感動を追体験することが可能なんです。本当に凄いことだとは思いませんか？