Return to 宇宙論パラメータと宇宙年齢

補足:Matplotlib でスケール因子のグラフを描く

導出については,以下のページ:

必要なモジュールの import

In [1]:
import numpy as np
import matplotlib.pyplot as plt

# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
In [2]:
# デフォルト設定
config = {
    'axes.labelsize': 'x-large',
    'mathtext.fontset': 'cm'
}
plt.rcParams.update(config)

$t = 0$ からの $a(t)$ の立ち上がりを揃えたグラフの例

異なる $\Omega_{\rm m}$ の場合のスケール因子の時間変化のグラフを,$t=0$ でのスケール因子 $a(t)$ の傾きを揃えて描く場合。同じようにビッグバンで始まった宇宙の膨張が,$\Omega_{\rm m}$ の値によってその後の膨張の仕方に違いがあらわれ,ある場合には途中で膨張が止まって収縮に転じたり,ある場合には永久に膨張が続いたりするんだなぁ… ということを理解するのに適切なグラフ。

$\Omega_{\Lambda} = 0, \Omega_{\rm m} > 1$ すなわち $k > 0$ の場合

$\Omega_{\rm m} \rightarrow \Omega$ として

\begin{eqnarray}
a_1(u, \Omega) = \frac{a}{a_0}
&=& \frac{\Omega}{2 (\Omega -1)}
\left(1 -\cos\left(\sqrt{\Omega-1} u\right)\right)\\
t_1(u, \Omega) = H_0 t &=& \frac{\Omega}{2 (\Omega -1) }
\left(u -\frac{\sin\left(\sqrt{\Omega-1} u\right)}{\sqrt{\Omega-1}}\right)
\end{eqnarray}

となりそうだが,$|u| \ll 1$ での振るまいが後述の $\Omega_{\rm m} = 1$ の場合のように($\Omega$ の値によらずに)

\begin{eqnarray}
a_1 &\simeq& \frac{u^2}{4} \\
t_1 &\simeq& \frac{u^3}{12}
\end{eqnarray}

となるためには,以下のように縦横等倍で縮尺を変えてやればよい。

\begin{eqnarray}
a_1(u, \Omega) \equiv \frac{a}{a_0} \times \Omega^{-1}
&=& \frac{1}{2 (\Omega -1)}
\left(1 -\cos\left(\sqrt{\Omega-1} u\right)\right)\\
t_1(u, \Omega) \equiv H_0 t \times \Omega^{-1}&=& \frac{1}{2 (\Omega -1) }
\left(u -\frac{\sin\left(\sqrt{\Omega-1} u\right)}{\sqrt{\Omega-1}}\right)
\end{eqnarray}

In [3]:
def a1(u, Omega):
    return 1/(2*(Omega-1))*(1-np.cos(np.sqrt(Omega-1)*u))
def t1(u, Omega):
    return 1/(2*(Omega-1))*(u-np.sin(np.sqrt(Omega-1)*u)/np.sqrt(Omega-1))

$\Omega_{\Lambda} = 0, \Omega_{\rm m} < 1$ すなわち $k < 0$ の場合

同様に

\begin{eqnarray}
a_2(u, \Omega) \equiv \frac{a}{a_0}\times \Omega^{-1}
&=& \frac{1}{2 (1-\Omega)}
\left(\cosh\left(\sqrt{1-\Omega} u\right) -1\right)
\\
t_2(u, \Omega) \equiv H_0 t \times \Omega^{-1}&=& \frac{1}{2 (1 -\Omega) }
\left(\frac{\sinh\left(\sqrt{1-\Omega} u\right)}{\sqrt{1-\Omega}}- u\right)
\end{eqnarray}

In [4]:
def a2(u, Omega):
    return 1/(2*(1-Omega))*(np.cosh(np.sqrt(1-Omega)*u)-1)
def t2(u, Omega):
    return 1/(2*(1-Omega))*(np.sinh(np.sqrt(1-Omega)*u)/np.sqrt(1-Omega)-u)

$\Omega_{\Lambda} = 0, \Omega_{\rm m} = 1$ すなわち $k = 0$ の場合

\begin{eqnarray}
a_3(u) \equiv \lim_{\Omega\rightarrow 1} a_1(u, \Omega) &=& \frac{u^2}{4} \\
t_3(u) \equiv \lim_{\Omega\rightarrow 1} t_1(u, \Omega) &=& \frac{u^3}{12}
\end{eqnarray}

In [5]:
def a3(u):
    return u**2/4
def t3(u):
    return u**3/12
In [6]:
def a(u, Omega):
    if Omega > 1:
        return a1(u, Omega)
    elif Omega == 1:
        return a3(u)
    else:
        return a2(u, Omega)
    
def t(u, Omega):
    if Omega > 1:
        return t1(u, Omega)
    elif Omega == 1:
        return t3(u)
    else:
        return t2(u, Omega)
In [7]:
Omega1 = 1.1
Omega2 = 0.9
# 0 <= u <= u1 まで 
u1 = 2*np.pi/np.sqrt(Omega1 -1)
urange = np.linspace(0, u1, 100)
# 横軸表示範囲
tend = t1(u1, Omega1)*1.02

plt.plot(t(urange, Omega2), a(urange, Omega2), 
         color='blue');
plt.plot(t(urange, 1     ), a(urange, 1    ),
         color='black');
plt.plot(t(urange, Omega1), a(urange, Omega1),
         color='red');

# 横軸縦軸の表示範囲
plt.xlim(0, tend)
plt.ylim(0, 42)
# 座標軸の目盛を非表示に
plt.xticks([0])
plt.yticks([0])

# 数式フォント設定
plt.title(r"$\Omega_{\Lambda}=0$ の場合のスケール因子の時間発展")
plt.xlabel(r"$t$")
plt.ylabel(r"$a(t)$")

# 曲線の近くにラベルと数式フォント設定
plt.text(60, 38, r"$\Omega_{\rm{m}} < 1\ \ (k < 0)$", 
        {'color':'blue', 'size':'x-large'})
plt.text(60, 18, r"$\Omega_{\rm{m}} = 1\ \ (k = 0)$", 
        {'color':'black', 'size':'x-large'})
plt.text(60, 10.5, r"$\Omega_{\rm{m}} > 1\ \ (k > 0)$", 
        {'color':'red', 'size':'x-large'});

$t = t_0$ で $a(t_0)$ と $H_0 = \frac{\dot{a}}{a}|_{t_0}$ を揃えたグラフの例

異なる $\Omega_{\rm m}$ の場合のスケール因子の時間変化のグラフを,現在時刻 $t=t_0$ でのスケール因子 $a(t)$ の傾きを揃えて描く場合。現在時刻でのスケール因子の傾きを表すハッブル定数 $H_0$ の値が同じでも,時間を遡るとやがて $a(t)$ がゼロになる時刻すなわち宇宙年齢が異なるのだなぁ… ということを理解するのに便利なグラフ。

$\Omega_{\Lambda} = 0, \Omega_{\rm m} > 1$ すなわち $k > 0$ の場合

\begin{eqnarray}
\frac{a}{a_0} \equiv x
&=& \frac{\Omega_{\rm m}}{2 (\Omega_{\rm m} -1)}
\left(1 -\cos\left(\sqrt{k} \eta\right)\right)
\\
H_0 t &=& \frac{\Omega_{\rm m}}{2 (\Omega_{\rm m} -1)^{\frac{3}{2}} }
\left(\sqrt{k} \eta -\sin\left(\sqrt{k} \eta\right)\right)
\end{eqnarray}

これから
\begin{eqnarray}
\cos\sqrt{k} \eta &=& 1 – \frac{2 (\Omega_{\rm m} -1)}{\Omega_{\rm m}}x\\
\sqrt{k} \eta &=& \cos^{-1} \left(1 – \frac{2 (\Omega_{\rm m} -1)}{\Omega_{\rm m}}x \right)
\end{eqnarray}

$\Omega_{\rm m} \rightarrow \Omega$ として

\begin{eqnarray}
\therefore\ \ T_1(x, \Omega) &\equiv& H_0 t = \frac{\Omega}{2 (\Omega -1)^{\frac{3}{2}} }
\left(\sqrt{k} \eta -\sin \sqrt{k} \eta \right) \\
&=& \frac{\Omega}{2 (\Omega -1)^{\frac{3}{2}} }
\left(\cos^{-1} \left(1 – \frac{2 (\Omega -1)}{\Omega} x \right)
– \sqrt{1 – \left(1 – \frac{2 (\Omega -1)}{\Omega} x \right)^2} \right)
\end{eqnarray}

として,$H_0 t$ を $x$ で表す。

In [8]:
def T1(x, Omega):
    return (Omega/(2*(Omega-1)*np.sqrt(Omega-1)) * 
            (np.arccos(1 - 2*(Omega-1)/Omega * x) 
             - np.sqrt(1-(1-2*(Omega-1)/Omega * x)**2)))

$\Omega_{\Lambda} = 0, \Omega_{\rm m} < 1$ すなわち $k < 0$ の場合

\begin{eqnarray}
\frac{a}{a_0} \equiv x
&=&\frac{\Omega_{\rm m}}{2 (1-\Omega_{\rm m} )}
\left(\cosh \sqrt{k} \eta -1\right)a
\end{eqnarray}

これから
\begin{eqnarray}
\cosh\sqrt{k} \eta &=& 1 + \frac{2 (1-\Omega_{\rm m})}{\Omega_{\rm m}}x \\
\sqrt{k} \eta &=& \cosh^{-1} \left(1 + \frac{2 (1-\Omega_{\rm m})}{\Omega_{\rm m}}x\right)
\end{eqnarray}

$\Omega_{\rm m} \rightarrow \Omega$ として

\begin{eqnarray}
\therefore\ \
T_2(x, \Omega) &\equiv& H_0 t = \frac{\Omega}{2 (1-\Omega)^{\frac{3}{2}} }
\left(\sinh \sqrt{k} \eta – \sqrt{k} \eta \right) \\
&=& \frac{\Omega}{2 (1-\Omega)^{\frac{3}{2}} }
\left(
\sqrt{\left(1 + \frac{2 (1-\Omega)}{\Omega}x \right)^2 – 1} –
\cosh^{-1} \left(1 + \frac{2 (1-\Omega)}{\Omega}x\right) \right)
\end{eqnarray}

として,$H_0 t$ を $x$ で表す。

In [9]:
def T2(x, Omega):
    return (Omega/(2*(1-Omega)*np.sqrt(1-Omega)) * 
            (np.sqrt((1+2*(1-Omega)/Omega * x)**2 -1) 
             -np.arccosh(1 + 2*(1-Omega)/Omega * x)))

独り言:あぁそう。NumPy は $\cosh^{-1} x$ を np.arccosh(x) としているんだぁ。逆双曲線関数は arc と書くべきではないんだよなぁ。

$\Omega_{\Lambda} = 0, \Omega_{\rm m} = 1$ すなわち $k = 0$ の場合

\begin{eqnarray}
\frac{a}{a_0} = x &=& \left(\frac{3}{2}H_0 t \right)^{\frac{2}{3}}
\end{eqnarray}$$\therefore\ \ T_3(x) \equiv H_0 t = \frac{2}{3} x^{\frac{3}{2}}$$

In [10]:
def T3(x):
    return 2/3 * x*np.sqrt(x)
In [11]:
def T(x, Omega):
    if Omega > 1:
        return T1(x, Omega)
    elif Omega == 1:
        return T3(x)
    else:
        return T2(x, Omega)    
In [12]:
Omega1 = 2
Omega2 = 0.1

x = np.linspace(0, 3, 100)
x1 = np.linspace(0, 2, 100)

plt.plot(T(x, Omega2)-T(1, Omega2), x, color='blue',
         label=r"$\Omega_{\rm m} < 1\ \ (k < 0)$")
plt.plot(T(x, 1     )-T(1, 1     ), x, color='black',
         label=r"$\Omega_{\rm m} = 1\ \ (k = 0)$")
plt.plot(T(x1, Omega1)-T(1, Omega1), x1, color='red',
         label=r"$\Omega_{\rm m} > 1\ \ (k > 0)$")

# 横軸縦軸の表示範囲
plt.xlim(-1, 1.6)
plt.ylim(0, 3)
# グリッド(格子線)の表示
plt.grid(which="major", linestyle = 'dotted')

# 数式フォント設定
plt.title(r"$\Omega_{\Lambda}=0$ の場合のスケール因子の時間発展")
plt.xlabel(r"$H_0 (t -t_0)$")
plt.ylabel(r"$a(t)/a_0$")
plt.legend(prop={"size":"large"})

# 現在 t=t0 
plt.axhline(1, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6)

# 副目盛も表示
plt.minorticks_on()
# 副目盛には grid をつけない
plt.grid(False, which="minor");
plt.legend(prop={"size":"large"});

$k = 0, \ \Omega_{\Lambda} > 0$ の場合の追加

宇宙定数がある場合のスケール因子の時間変化についても追加しておく。

\begin{eqnarray}
\frac{a}{a_0} = x &=&
\left\{\sqrt{\frac{\Omega_{\rm m}}{1-\Omega_{\rm m}}}
\sinh\left(\frac{3\sqrt{1-\Omega_{\rm m}}}{2} H_0 t\right)\right\}^{\frac{2}{3}} \end{eqnarray}

より $\Omega_{\rm m} \rightarrow \Omega$ として

$$T_4(x, \Omega) \equiv H_0 t = \frac{2}{3\sqrt{1-\Omega}}
\sinh^{-1} \left( \sqrt{\frac{1-\Omega}{\Omega}} x^{\frac{3}{2}}\right)$$

In [13]:
def T4(x, Omega):
    return 2/(3*np.sqrt(1-Omega)) * np.arcsinh(np.sqrt((1-Omega)/Omega) * x*np.sqrt(x))

独り言:あぁそう。NumPy は $\sinh^{-1} x$ を np.arcsinh(x) としているんだぁ。逆双曲線関数は arc と書くべきではないんだよなぁ。

また,$t = 0$ からの $a(t)$ の立ち上がりを揃えたグラフにするには,以下のように定義すればよいであろう。

In [14]:
def a4(u, Omega):
    return (np.sqrt(1/(1-Omega))*np.sinh(3*np.sqrt(1-Omega)*u**3/12/2))**(2/3)
def t4(u, Omega):
    return u**3/12

$t = 0$ からの $a(t)$ の立ち上がりを揃えたグラフの例

In [15]:
Omega1 = 1.1
Omega2 = 0.7
# 0 <= u <= u1 まで 
u1 = 2*np.pi/np.sqrt(Omega1 -1)
urange = np.linspace(0, u1, 100)
# 横軸表示範囲
tend = t1(u1, Omega1)*1.1

plt.plot(t4(urange, 0.999), a4(urange, 0.999), 
         color='purple');
plt.plot(t(urange, Omega2), a(urange, Omega2), 
         color='blue');
plt.plot(t(urange, 1     ), a(urange, 1    ),
         color='black');
plt.plot(t(urange, Omega1), a(urange, Omega1),
         color='red');

# 横軸縦軸の表示範囲
plt.xlim(0, tend)
plt.ylim(0, 120)
# 座標軸の目盛を非表示に
plt.xticks([0])
plt.yticks([0])

# 数式フォント設定
plt.title(r"スケール因子の時間発展")
plt.xlabel(r"$t$")
plt.ylabel(r"$a(t)$")

# 曲線の近くにラベルと数式フォント設定
plt.text(80, 75, r"$k=0, \Omega_{\Lambda} >0$", 
        {'color':'purple', 'size':'large'})
plt.text(80, 62, r"$k<0, \Omega_{\Lambda} =0$", 
        {'color':'blue', 'size':'large'})
plt.text(80, 30, r"$k=0, \Omega_{\Lambda} =0$", 
        {'color':'black', 'size':'large'})
plt.text(80, 9, r"$k>0, \Omega_{\Lambda} =0$", 
        {'color':'red', 'size':'large'});

$t = t_0$ で $a(t_0)$ と $H_0 = \frac{\dot{a}}{a}|_{t_0}$ を揃えたグラフの例

In [16]:
Omega1 = 2
Omega2 = 0.3

x = np.linspace(0, 3, 100)
x1 = np.linspace(0, 2, 100)

plt.plot(T4(x, Omega2)-T4(1, Omega2), x, color='purple',
         label=r"$\Omega_{\rm m} =%.1f, \ \Omega_{\Lambda} =%.1f$" % (Omega2, 1-Omega2))
plt.plot(T(x, Omega2)-T(1, Omega2), x, color='blue',
         label=r"$\Omega_{\rm m} =%.1f, \ \Omega_{\Lambda} =%.1f$" % (Omega2, 0))
plt.plot(T(x, 1     )-T(1, 1     ), x, color='black',
         label=r"$\Omega_{\rm m} =%.1f, \ \Omega_{\Lambda} =%.1f$" % (1, 0))
plt.plot(T(x1, Omega1)-T(1, Omega1), x1, color='red',
         label=r"$\Omega_{\rm m} =%.1f, \ \Omega_{\Lambda} =%.1f$" % (Omega1, 0))

# 横軸縦軸の表示範囲
plt.xlim(-1, 1.6)
plt.ylim(0, 3)
# グリッド(格子線)の表示
plt.grid(which="major", linestyle = 'dotted')

# 数式フォント設定
plt.title(r"スケール因子の時間発展")
plt.xlabel(r"$H_0 (t -t_0)$")
plt.ylabel(r"$a(t)/a_0$")
plt.legend(prop={"size":"large"})

# 現在 t=t0 
plt.axhline(1, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6)

# 副目盛も表示
plt.minorticks_on()
# 副目盛には grid をつけない
plt.grid(False, which="minor");
plt.legend(prop={"size":"large"});