Edited at

Astropyで宇宙論的な計算を行う


Introduction

このページは、内容としてはAstropyのCosmological Calculationsのページを和訳したものですが、ところどころで私なりのアレンジを加えています。

Astropyモジュールを使って、宇宙論的な計算を行う方法をご紹介します。


Installation

PyPIに登録されているので

pip install astropy

でインストールするか、Anacondaがインストールされていれば

conda install -c astropy sphinx-astropy

でインストールすることができます(Anacondaの場合は、すでにインストールされている場合もあります)。


初めてのastropy.cosmology


ハッブル定数

以下のようにすることで、WMAP9yrから推定された、redshift z=0でのハッブル定数を求めることができます。

>>> from astropy.cosmology import WMAP9

>>> WMAP9.H(0)
<Quantity 69.32 km / (Mpc s)>


固有距離

z=3における、1ミリ秒角の固有距離を、kpc単位で算出します。

>>> WMAP9.kpc_proper_per_arcmin(3)

<Quantity 472.97709236 kpc / arcmin>


共動距離

リストなどのシーケンスを渡すことで、一気に共動距離を求めることができます。以下はz=0, 1, ..., 9までの共動距離を求めています。

>>> WMAP9.comoving_distance(range(10))

<Quantity [ 0. , 3363.0706321 , 5291.73040581, 6503.90188771,
7344.44470431, 7968.80246409, 8455.45279875, 8848.30522052,
9173.9669946 , 9449.58596714] Mpc>


Flat Lambda CDM modelとニュートリノ種の実効数

平坦なΛ-CDMモデルでのニュートリノ種の実効数を求めることができます。ここではハッブル定数 H0, 物質の密度パラメータ Om0, CMB温度 Tcmb0を指定しています。

>>> from astropy.cosmology import FlatLambdaCDM

>>> cosmo = FlatLambdaCDM(H0=70, Om0=0.3, Tcmb0=2.725)
>>> cosmo.Neff
3.04

また、このときcosmo変数をprintしてみると

>>> print cosmo

FlatLambdaCDM(H0=70 km / (Mpc s), Om0=0.3, Tcmb0=2.725 K, Neff=3.04, m_nu=[0. 0. 0.] eV, Ob0=None)

という、オブジェクトになっていることがわかります。


.unitメソッド

>>> from astropy.cosmology import WMAP9

>>> H0 = WMAP9.H(0)
>>> H0.value, H0.unit
(69.32, Unit("km / (Mpc s)"))

とすると、WMAP 9yrにおけるハッブル定数の値と、その単位を表示することができます。


astropy.cosmologyを使ってみよう


import astropy.units as u

>>> from astropy.cosmology import FlatLambdaCDM

>>> import astropy.units as u
>>> cosmo = FlatLambdaCDM(H0=70 * u.km / u.s / u.Mpc, Tcmb0=2.725 * u.K, Om0=0.3)

u.km, u.s, u.Mpc, u.Kなどとして、上述のように単位を明記することを、astropy.cosmologyドキュメントでは推奨しています。


光度距離

以下のようにして光度距離を求めることができます。以下はz=4までの距離です。

>>> cosmo.luminosity_distance(4)

<Quantity 35842.353288 Mpc>


宇宙年齢

以下のようにして宇宙年齢を求めることができます。以下はz=0での宇宙年齢を求める例です。

>>> cosmo.age(0)

<Quantity 13.46170148 Gyr>

以下のように、リストのようなシーケンスを渡し、.valueメソッドを使うことで、リストの形で宇宙年齢を取得することも可能です。

>>> cosmo.age(range(10)).value

array([13.46170148, 5.74698037, 3.22268178, 2.10909267, 1.51293349,
1.15196553, 0.91440079, 0.74841291, 0.62710773, 0.53530494])


FLRW クラス

FLRWクラスには一様等方宇宙モデル(Friedmann-Lemaitre-Robertson-Walker cosmology)が記述されています。例えば、FlatLambdaCDM, LambdaCDM, wCDMなどが計算できます。


WMAP 7yrの観測結果から宇宙論パラメータたちを算出する

astropy.cosmologyには様々な宇宙論パラメータが揃っています。以下はWMAP 7yrの結果からz=0での臨界密度を算出する方法です。

>>> from astropy.cosmology import WMAP7

>>> WMAP7.critical_density(0)
<Quantity 9.30966846e-30 g / cm3>

以下では物質密度パラメータ、ダークエネルギー密度パラメータを算出しています。

>>> WMAP7.Om([0, 1, 2])

array([0.272 , 0.74898523, 0.90905237])

>>> WMAP7.Ode([0, 1, 2])

array([0.72791572, 0.25055061, 0.0901026 ])

WMAP 7yrでは光子とニュートリノを考慮しているため、これらを足し合わせても1にはなりません。


バリオン密度バラメータ

以下のように、z=0でのバリオンの密度パラメータOb0を設定することも可能です。

>>> from astropy.cosmology import FlatLambdaCDM

>>> FlatLambdaCDM(H0=70, Om0=0.3, Ob0=0.05)
FlatLambdaCDM(H0=70 km / (Mpc s), Om0=0.3, Tcmb0=0 K, Neff=3.04, m_nu=None, Ob0=0.05)

この場合、このクラスからはz=0におけるダークマター密度パラメータしか得ることができません。redshift毎の変化を見たい場合には、ダークマターの密度パラメータOdmとバリオンの密度パラメータObを使う必要があります。もしz=0でのバリオン密度パラメータOb0が指定されていない場合には、Odm関数は以下のようにValueErrorを返します。

>>> from astropy.cosmology import FlatLambdaCDM

>>> cosmo = FlatLambdaCDM(H0=70, Om0=0.3)
>>> cosmo.Odm(1)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/anaconda3/lib/python3.7/site-packages/astropy/cosmology/core.py", line 574, in Odm
raise ValueError("Baryonic density not set for this cosmology, "
ValueError: Baryonic density not set for this cosmology, unclear meaning of dark matter density


name変数

Cosmologicalインスタンスはname変数を持っており、この変数に宇宙論モデルを記述しておくことができます。

>>> from astropy.cosmology import FlatwCDM

>>> cosmo = FlatwCDM(name='SNLS3+WMAP7', H0=71.58, Om0=0.262, w0=-1.016)
>>> cosmo
FlatwCDM(name="SNLS3+WMAP7", H0=71.6 km / (Mpc s), Om0=0.262, w0=-1.02, Tcmb0=0 K, Neff=3.04, m_nu=None, Ob0=None)


サポートしているダークエネルギーのモデル

上のname変数の部分で紹介したスクリプトはダークエネルギーのモデル計算一例です。他の異なるダークエネルギーのモデルもサポートされています。詳細はSpecifying a dark energy modelを参照してください。


宇宙論パラメータの変更

宇宙論パラメータはイミュータブル(変更可能な変数)です。例えば、WMAP 9yrの結果のうち、z=0での物質密度パラメータOm0のみ変更したいとしましょう。その場合、cloneメソッドを使って新しいインスタンスを作成する方法が便利です。試しに

>>> from astropy.cosmology import WMAP9

>>> newcosmo = WMAP9.clone(name='WMAP9 modified', Om0=0.3141)

としてnewcosmoインスタンスを作成し、元のWMAP9とnewcosmoでのハッブルパラメータ、物質密度パラメータ、ダークエネルギー密度パラメータを確認してみましょう。

>>> WMAP9.H0

<Quantity 69.32 km / (Mpc s)>
>>> newcosmo.H0
<Quantity 69.32 km / (Mpc s)>
>>> WMAP9.Om0
0.2865
>>> newcosmo.Om0
0.3141
>>> WMAP9.Ode0
0.7134130719051658
>>> newcosmo.Ode0
0.6858130719051657


z_at_value

z_at_valueに、用いる宇宙論モデルの.ageメソッドと年齢を引数にすることで、その年齢におけるredshiftを算出することができます。以下の例ではPlanck 2013を宇宙論モデルに指定したときの、2Gyrの宇宙のredshiftを算出しています。

>>> from astropy.cosmology import Planck13, z_at_value

>>> z_at_value(Planck13.age, 2*u.Gyr)
3.198120965652982

試しにPlanck13.ageを表示してみます。

>>> Planck13.age

<bound method FLRW.age of FlatLambdaCDM(name="Planck13", H0=67.8 km / (Mpc s), Om0=0.307, Tcmb0=2.725 K, Neff=3.05, m_nu=[0. 0. 0.06] eV, Ob0=0.0483)>


すでに組み込まれている宇宙論モデルたち

astropy.cosmologyに組み込まれている宇宙論モデルは、cosmology.parameters.availableで確認することができます。それらを以下に表としてまとめました。

Name
Source
H0
Om
Flat

WMAP5
Komatsu et al. 2009
70.2
0.277
Yes

WMAP7
Komatsu et al. 2011
70.4
0.272
Yes

WMAP9
Hinshaw et al. 2013
69.3
0.287
Yes

Planck13
Planck Collab 2013, Paper XVI
67.8
0.307
Yes

Planck15
Planck Collab 2015, Papaer XIII
67.7
0.307
Yes

モデルの詳細は.__doc__で確認することができます。

>>> from astropy.cosmology import WMAP7

>>> WMAP7.__doc__
'WMAP7 instance of FlatLambdaCDM cosmology\n\n(from Komatsu et al. 2011, ApJS, 192, 18, doi: 10.1088/0067-0049/192/2/18. Table 1 (WMAP + BAO + H0 ML).)'


ダークエネルギーのモデル

リンクを参照


光子とニュートリノ

cosmologyクラスは光子のエネルギー密度・ニュートリノのエネルギー密度を計算することができます。デフォルトでは質量なしの仮定を行なっています。

z=0でのCMB温度T_cmb0, ニュートリノ実効数Neff, 各種のニュートリノ質量m_nuを設定することで、これらの計算を行うことができます。

以下はWMAP 7yrによって求められた、z=0での光子の密度パラメータOgamma0, ニュートリノの密度パラメータOnu0を求めています。

>>> from astropy.cosmology import WMAP7

>>> WMAP7.Ogamma0
4.985868989189756e-05
>>> WMAP7.Onu0
3.442275090310745e-05

Ogamma, Onuにシーケンスを渡すことで、一度に複数のredshiftにおける光子の密度パラメータ・ニュートリノ密度パラメータを求めることも可能です。

>>> WMAP7.Ogamma(range(10))

array([4.98586899e-05, 2.74583989e-04, 4.99898824e-04, 7.02950869e-04,
8.95950291e-04, 1.08437305e-03, 1.27046001e-03, 1.45521852e-03,
1.63914762e-03, 1.82251551e-03])
>>> z = [0, 1, 2]
>>> WMAP7.Onu(z)
array([3.44227509e-05, 1.89574501e-04, 3.45133270e-04])

もし宇宙論モデルから光子とニュートリノを除外したければ、Tcmb0=0とすることでそれが可能です。

>>> from astropy.cosmology import FlatLambdaCDM

>>> import astropy.units as u
>>> cos = FlatLambdaCDM(70.4*u.km/u.s/u.Mpc, 0.272, Tcmb0=0.0*u.K)
>>> cos.Ogamma0
0.0
>>> cos.Onu0
0.0

ニュートリノのみを無視したモデルを構築したければ、Neff=0とします。

>>> cos = FlatLambdaCDM(70.4, 0.272, Tcmb0=2.725, Neff=0)

>>> z = [0, 1, 2]
>>> cos.Ogamma(z)
array([4.98586899e-05, 2.74632798e-04, 5.00069284e-04])
>>> cos.Onu(z)
array([0., 0., 0.])

Neffでニュートリノ実効数を仮定した場合、各々のニュートリノ質量をインプットとして与える必要があります。もし一つしか値を設定しなかった場合には、全ての種類のニュートリノ質量が一律でその値になります。

>>> from astropy.cosmology import FlatLambdaCDM

>>> import astropy.units as u
>>> H0 = 70.4 * u.km / u.s / u.Mpc
>>> m_nu = 0 * u.eV
>>> cos = FlatLambdaCDM(H0, 0.272, Tcmb0=2.725, m_nu=m_nu)
>>> cos.has_massive_nu
False
>>> cos.m_nu
<Quantity [0., 0., 0.] eV>
>>> m_nu = [0.0, 0.05, 0.1] * u.eV
>>> cos = FlatLambdaCDM(H0, 0.272, Tcmb0=2.725, m_nu=m_nu)
>>> cos.has_massive_nu
True
>>> cos.m_nu
<Quantity [0. , 0.05, 0.1 ] eV>
>>> z = [0, 1, 15]
>>> cos.Onu(z)
array([0.00327 , 0.00896814, 0.01257904])
>>> cos.Onu(1) * cos.critical_density(1)
<Quantity 2.44437837e-31 g / cm3>

全てではありませんが、一応この記事はここまで。