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

# Introduction

このページは、内容としては[AstropyのCosmological Calculationsのページ](http://docs.astropy.org/en/stable/cosmology/index.html#using-astropy-cosmology)を和訳したものですが、ところどころで私なりのアレンジを加えています。
Astropyモジュールを使って、宇宙論的な計算を行う方法をご紹介します。

# Installation

[PyPI](https://pypi.org)に登録されているので

```
pip install astropy
```

でインストールするか、[Anaconda](https://www.anaconda.com)がインストールされていれば

```
conda install -c astropy sphinx-astropy
```

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

# 初めてのastropy.cosmology

## ハッブル定数

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

```python
>>> from astropy.cosmology import WMAP9
>>> WMAP9.H(0)
<Quantity 69.32 km / (Mpc s)>
```

## 固有距離

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

```python
>>> WMAP9.kpc_proper_per_arcmin(3)
<Quantity 472.97709236 kpc / arcmin>
```

## 共動距離

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

```python
>>> 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を指定しています。 

```python
>>> from astropy.cosmology import FlatLambdaCDM
>>> cosmo = FlatLambdaCDM(H0=70, Om0=0.3, Tcmb0=2.725)
>>> cosmo.Neff
3.04
```

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

```python
>>> 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メソッド

```python
>>> 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

```python
>>> 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までの距離です。

```python
>>> cosmo.luminosity_distance(4)
<Quantity 35842.353288 Mpc>
```

## 宇宙年齢

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

```python
>>> cosmo.age(0)
<Quantity 13.46170148 Gyr>
```

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

```python
>>> 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クラス](http://docs.astropy.org/en/stable/api/astropy.cosmology.FLRW.html#astropy.cosmology.FLRW)には一様等方宇宙モデル(Friedmann-Lemaitre-Robertson-Walker cosmology)が記述されています。例えば、FlatLambdaCDM, LambdaCDM, wCDMなどが計算できます。

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

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

```python
>>> from astropy.cosmology import WMAP7
>>> WMAP7.critical_density(0)
<Quantity 9.30966846e-30 g / cm3>
```

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

```python
>>> WMAP7.Om([0, 1, 2])
array([0.272     , 0.74898523, 0.90905237])
```

```python
>>> WMAP7.Ode([0, 1, 2])
array([0.72791572, 0.25055061, 0.0901026 ])
```

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

## バリオン密度バラメータ

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

```python
>>> 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`を返します。

```python
>>> 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変数を持っており、この変数に宇宙論モデルを記述しておくことができます。

```python
>>> 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](http://docs.astropy.org/en/stable/cosmology/index.html#specifying-a-dark-energy-model)を参照してください。

## 宇宙論パラメータの変更

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

```python
>>> from astropy.cosmology import WMAP9
>>> newcosmo = WMAP9.clone(name='WMAP9 modified', Om0=0.3141)
```

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

```python
>>> 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を算出しています。

```python
>>> from astropy.cosmology import Planck13, z_at_value
>>> z_at_value(Planck13.age, 2*u.Gyr)
3.198120965652982
```

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

```python
>>> 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__`で確認することができます。

```python
>>> 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).)'
```

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

[リンク](http://docs.astropy.org/en/stable/cosmology/index.html#specifying-a-dark-energy-model)を参照

## 光子とニュートリノ

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

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

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

```python
>>> from astropy.cosmology import WMAP7
>>> WMAP7.Ogamma0
4.985868989189756e-05
>>> WMAP7.Onu0
3.442275090310745e-05
```

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

```python
>>> 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`とすることでそれが可能です。

```python
>>> 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`とします。

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

```python
>>> 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>
```

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