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>
全てではありませんが、一応この記事はここまで。