Python で Planck 2018 results から宇宙年齢と物質密度を計算する

Python を電卓として使って宇宙年齢や宇宙空間の物質密度を計算する例。Python でなくても，$\sqrt{x}$ や $\tanh^{-1} x$ ができる電卓があればいいです。

Planck 2018 の宇宙論パラメータ

Planck 2018 results の宇宙論パラメータは以下の通り。

\begin{eqnarray}
H_0 &=& (67.4 \pm 0.5) \ \mbox{km/s/Mpc} \\
\Omega_{\rm m} &=& 0.315 \pm 0.007 \\
\Omega_{\Lambda} &=& 1 – \Omega_{\rm m}
\end{eqnarray}

宇宙年齢

$\Omega_{\rm m} + \Omega_{\Lambda} = 1$ すなわち $k = 0$ の場合の宇宙年齢はここ に書いているように

$$
H_0 t_0 = \frac{2}{3(\sqrt{1-\Omega_{\rm m} })}\tanh^{-1} \sqrt{1-\Omega_{\rm m} } \quad \mbox{for}\ \ \Omega_{\rm m} < 1$$

より

$$t_0 = \frac{1}{H_0} \times \frac{2}{3(\sqrt{1-\Omega_{\rm m} })}\tanh^{-1} \sqrt{1-\Omega_{\rm m} }$$

必要なモジュールの import

In [1]:
from sympy import *
from sympy.abc import *
from sympy import I, pi, E

import numpy as np

パーセクの定義とハッブル定数

$$1 \,\mbox{pc} \equiv \frac{\mbox{au}}{\mbox{1 arcsec}} = \frac{648000}{\pi} \,\mbox{au}$$

In [2]:
au = 149597870700 * m
# Mpc だから 1E6 pc で，km だから 1E3 m。
Mpc = 648000/float(pi) * au * 1E6
km = 1E3 * m

H0 = 67.4 * km/(s * Mpc)
H0
Out[2]:
$\displaystyle \frac{2.1842852410855 \cdot 10^{-18}}{s}$

ハッブル定数の逆数を億年で表すと…

In [3]:
okunen = 60*60*24*365.25 * 1E8 * s
print('ハッブル年齢は %.1f 億年' % (1/H0/okunen))

ハッブル年齢は 145.1 億年

SymPy で計算する宇宙年齢

In [4]:
def t0(H, Omega):
    """
    H0 (km/(s Mpc)) と Omegam から宇宙年齢を億年で返す。
    """
    au = 149597870700 * m
    # Mpc だから 1E6 pc で，km だから 1E3 m。
    Mpc = 648000/float(pi) * au * 1E6
    km = 1E3 * m
    H0 = H * km/(s * Mpc)
    
    okunen = 60*60*24*365.25 * 1E8 * s

    # SymPy の関数を使って計算。
    t = 2/(3*H0*sqrt(1-Omega))*atanh(sqrt(1-Omega))/okunen
    print('宇宙年齢は %.1f 億年' % t)
In [5]:
H_0 = 67.4
Omega_m = 0.315
t0(H_0, Omega_m)

宇宙年齢は 138.0 億年

NumPy で計算する宇宙年齢

In [6]:
def T0(H, Omega):
    """
    H0 (km/(s Mpc)) と Omegam から宇宙年齢を億年で返す。
    """
    au = 149597870700 * m
    # Mpc だから 1E6 pc で，km だから 1E3 m。
    Mpc = 648000/float(np.pi) * au * 1E6
    km = 1E3 * m
    H0 = H * km/(s * Mpc)
    
    okunen = 60*60*24*365.25 * 1E8 * s
    
    # NumPy の関数を使って計算。
    t = 2/(3*H0*np.sqrt(1-Omega))*np.arctanh(np.sqrt(1-Omega))/okunen
    print('宇宙年齢は %.1f 億年' % t)

問題 1

誤差を考えてハッブルパラメータの最大値，最小値をそれぞれ $H_{\rm max}, \ H_{\rm min}$，密度パラメータの最大値，最小値をそれぞれ $\Omega_{\rm max}, \ \Omega_{\rm min}$ とする。

\begin{eqnarray}
H_{\rm max} &=& 67.4 + 0.5\\
H_{\rm min} &=& 67.4 – 0.5\\
\Omega_{\rm max} &=& 0.315 + 0.007\\
\Omega_{\rm min} &=& 0.315 – 0.007\\
\end{eqnarray}

ハッブルパラメータを大きくすると，宇宙年齢 $t_0$ の値は大きくなるのか小さくなるのか，また密度パラメータを大きくすると，宇宙年齢 $t_0$ の値は大きくなるのか小さくなるのか，このへんを理解した上で，$t_0$ の最大値，最小値を求め，$138 \pm ??$ 億年の形に表せ。

In [7]:
# 解答例：
# 宇宙年齢は H0 に反比例し，Omega_mの単調減少関数だから...
# 宇宙年齢が最大となるのは
H_0 = 67.4 - 0.5
Omega_m = 0.315 - 0.007
print('最大となる', end='')
T0(H_0, Omega_m)

# 宇宙年齢が最小となるのは
H_0 = 67.4 + 0.5
Omega_m = 0.315 + 0.007
print('最小となる', end='')
T0(H_0, Omega_m)

最大となる宇宙年齢は 139.9 億年
最小となる宇宙年齢は 136.1 億年

臨界密度

ここ に書いているように

$$
\rho_{\rm cr}\equiv \frac{3 H_0^2}{8\pi G}
$$

宇宙空間の物質密度

宇宙空間の物質密度 $\rho$ は
$$\rho = \rho_{\rm cr} \Omega_{\rm m}$$から求める。

万有引力定数 $G$ は

$$ G = 6.674 \times 10^{-11} \, \frac{\mbox{m}^3}{\mbox{kg}\cdot\mbox{s}^2}$$

また，$H_0$ は別ページにまとめているように $s^{-1}$ の次元を持つので，適切に単位を合わせて計算すれば確かに $\rho_{\rm cr}$ は
質量密度 $\mbox{kg}/\mbox{m}^3$ の次元をもつことがわかる。

In [8]:
var('kg')
G = 6.674E-11 * m**3/(kg * s**2)

au = 149597870700 * m
# Mpc だから 1E6 pc で，km だから 1E3 m。
Mpc = 648000/float(pi) * au * 1E6
km = 1E3 * m
H0 = 67.4 * km/(s * Mpc)
Omega_m = 0.315
In [9]:
rho = 3*H0**2/(8*np.pi*G) * Omega_m
rho
Out[9]:
$\displaystyle \frac{2.68797019689748 \cdot 10^{-27} kg}{m^{3}}$

問題 2

現在の宇宙空間の物質密度は次のうちのどれに最も近いか。

  1. $1\,\mbox{cm}^3$ に電子が数個程度
  2. $1\,\mbox{cm}^3$ に陽子が数個程度
  3. $1\,\mbox{m}^3$ に電子が数個程度
  4. $1\,\mbox{m}^3$ に陽子が数個程度
  5. $1\,\mbox{km}^3$ に電子が数個程度
  6. $1\,\mbox{km}^3$ に陽子が数個程度

参考

陽子，電子の質量を以下のリンク等から確認しておくこと。