Pythonによる膨張宇宙と宇宙背景放射軌道の計算

１． はじめに

宇宙は約138億年前にインフレーションに続くビッグバン膨張により生成されたとされている．図1は，参考文献(1)で紹介されている宇宙の歴史である．

　　　　　　　　　　図1　宇宙の歴史

宇宙年齢38万年は宇宙の晴れ上がり時期とされ，この時期の光は地球に宇宙背景放射として届き，宇宙形成解明のための手がかりとして盛んに研究されている．
　本報告では，標準宇宙モデルによる宇宙膨張計算方法を改良し，宇宙背景放射の光(電波)の伝搬軌跡をpythonにより計算するプログラムを作成したので報告する．

２． 問題設定

1)標準宇宙モデルによる膨張宇宙の定式化とpythonによる計算
2)宇宙背景放射の伝搬 計算手順とpythonによる計算
3)解析システム
　　RaspberryPi4＋Python

3． 標準宇宙モデルによる膨張宇宙の定式化とpythonによる計算

3.1 目的

　標準宇宙モデルは変数としてスケールファクタ:a を使っているが，宇宙背景放射の伝搬を計算するときは変数として時間:t を使う方が便利である．ここでは，aをtに変換し，膨張宇宙を変数として時間:tに対して計算，作図する方法を説明する．

3.2 定式化

参考文献(2)をまとめた結果を下記に示す．
　
上式から分かるように膨張宇宙の半径は，変数としてスケールファクタ:a (a=0～1.0)に対して計算されている．
　時間:H0tとスケールファクタ:aは下式により変換出来る．
　
したがって，
　時間:H0t　⇔　スケールファクタ:a (a=0～1.0) 　⇔　膨張宇宙の半径
　は，一対一で対応しているので，時間:H0tに対して膨張宇宙の半径を作図すれば時間軸に対する図を作成できる．
　今後の計算に際して，パラメータは下記を使用した．
　

3.3 pythonプログラム

　作成したプログラムは下記である．適宜コメントを入れておく．

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

Mpc=3260000
c=300000
H0=67.15
c_H0=Mpc*c/H0*10**-8
print("c_H0=",c_H0)
#
# スケールファクタ:a (a=0～1.0) 　⇔　膨張宇宙の半径
#
fig = plt.figure(figsize=(6, 18))
ax = fig.add_subplot(111)

ax.set_xlim(0, 500)
ax.set_ylim(0, 150)
ax.grid(which = "major", axis = "x", color = "green", alpha = 0.8,
        linestyle = "--", linewidth = 1)
ax.grid(which = "major", axis = "y", color = "green", alpha = 0.8,
        linestyle = "--", linewidth = 1)
#　
y = lambda x: 1/np.sqrt(0.317/x+0.683*x**2)
y1 = lambda x: 1/np.sqrt(0.317*x+0.683*x**4)
#
N=500
zp=0
s = np.linspace(zp, 1.0, N)
H0t= np.linspace(0.0001, 1.0, N)
LH= np.linspace(0.0001, 1.0, N)
#
for i in range(N):
    if i==0:
        H0t[i]=0.0038/c_H0
        LH[i]=0.5/c_H0
    else:
        s0=s[i]
        integ0 = integrate.quad(y, zp, s0)
# 年齢H0t[0]を38万年とする．
        H0t[i]=integ0[0]+0.0038/c_H0
        integ0=integrate.quad(y1, zp, s0)
# 年齢H0t[0]を38万年の宇宙半径を0.5億光年とする．
        LH[i]=integ0[0]+0.5/c_H0
        if i<10: print("sf=",i, s[i], "H0t=",H0t[i]*c_H0, "LH=",LH[i]*c_H0)

for i in range(N):
        s0=s[i]
        integ0 = integrate.quad(y, 0.001, s0)
# 年齢H0t[0]を38万年とする．
        H0t[i]=integ0[0]+0.0038/c_H0
        integ0=integrate.quad(y1, 0.001, s0)
# 年齢H0t[0]を38万年の宇宙半径を0.5億光年とする．
        LH[i]=integ0[0]+0.5/c_H0
#        if i<10: print("sf=",i, s[i], "H0t=",H0t[i]*c_H0, "LH=",LH[i]*c_H0)
#
plt.plot(LH*c_H0, s*c_H0, color = "green")
plt.ylabel("Scale factor: handred million year")
plt.xlabel("Space: handred million year *c")
#
# 時間:H0t　⇔　スケールファクタ:a (a=0～1.0) 
#
fig = plt.figure(figsize=(6, 18))
ax = fig.add_subplot(111)
# 
ax.set_xlim(0, 150)
ax.set_ylim(0, 150)
#
ax.grid(which = "major", axis = "x", color = "green", alpha = 0.8,
        linestyle = "--", linewidth = 1)
ax.grid(which = "major", axis = "y", color = "green", alpha = 0.8,
        linestyle = "--", linewidth = 1)
plt.ylabel("Scale Factor: handred million year")
plt.xlabel("Time H0t: handred million year")
plt.plot(H0t*c_H0, s*c_H0, color = "green")
#
# 時間:H0t　⇔　膨張宇宙の半径
#
fig = plt.figure(figsize=(6, 18))
ax = fig.add_subplot(111)
# 
ax.set_xlim(0, 500)
ax.set_ylim(0, 150)
#
ax.grid(which = "major", axis = "x", color = "green", alpha = 0.8,
        linestyle = "--", linewidth = 1)
ax.grid(which = "major", axis = "y", color = "green", alpha = 0.8,
        linestyle = "--", linewidth = 1)
plt.plot(LH*c_H0, H0t*c_H0, color = "green")
plt.ylabel("Time H0t: handred million year")
plt.xlabel("Space: handred million year *c")
# 
plt.show()

　図2 スケールファクタ:a (a=0～1.0) 　⇔　膨張宇宙の半径

　図3 時間:H0t　⇔　スケールファクタ:a (a=0～1.0)

　　　図4 時間:H0t　⇔　膨張宇宙の半径

4. 宇宙背景放射の伝搬 計算手順とpythonによる計算

4.1 目的

　宇宙年齢38万年の宇宙の光の一部が宇宙膨張 (赤方偏移は約1000)に逆らいながら138億年かけて現在地球に届いている．これが宇宙背景放射と言われている．宇宙年齢38万年の宇宙の大きさは現在の宇宙の大きさ約500億光年の1/1000の0.5億光年であり，宇宙膨張がなければ0.5億年で地球に届くはずであるが，宇宙膨張のため光は膨張に流されて，地球に到達するのに長い時間が必要になる．ここでは標準宇宙モデルを用いて宇宙背景放射の伝達経路を計算する．

4.2 計算手順

　計算の基本的なアルゴリズムは参考文献(3)を参考とした．
　図4に計算手順のイメージを示す．
　膨張開始初期は，ratio[i]≒1のため膨張速度 > 光の速度のため光の軌跡も膨張するが,膨張速度が低下し始め，ratio[i]<1となると膨張速度 < 光の速度になり光の軌跡は収束する．
　
　　　図4 計算手順のイメージ

4.3 pythonプログラム

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

Mpc=3260000
c=300000
H0=67.15
c_H0=Mpc*c/H0*10**-8
#print("c_H0=",c_H0)
#
fig = plt.figure(figsize=(6, 18))
ax = fig.add_subplot(111)
# 
ax.set_xlim(0, 500)
ax.set_ylim(0, 150)
#
ax.grid(which = "major", axis = "x", color = "green", alpha = 0.8,
        linestyle = "--", linewidth = 1)
ax.grid(which = "major", axis = "y", color = "green", alpha = 0.8,
        linestyle = "--", linewidth = 1)
#
y = lambda x: 1/np.sqrt(0.317/x+0.683*x**2)
y1 = lambda x: 1/np.sqrt(0.317*x+0.683*x**4)

N=500
zp=0.0
s = np.linspace(zp, 1.0, N)
H0t= np.linspace(0.0001, 1.0, N)
LH= np.linspace(0.0001, 1.0, N)
star0= np.linspace(0.0001, 1.0, N)
star1= np.linspace(0.0001, 1.0, N)
v= np.linspace(0.0001, 1.0, N)
ratio=np.linspace(0.0001, 1.0, N)

integ0=integrate.quad(y1, 0.0, 1)
print(c_H0*integ0[0])
#print(t)

for i in range(N):
    if i==0:
        H0t[i]=0.0038/c_H0
        LH[i]=0.5/c_H0
    else:
        s0=s[i]
        integ0 = integrate.quad(y, zp, s0)
        H0t[i]=integ0[0]+0.0038/c_H0
        integ0=integrate.quad(y1, zp, s0)
        LH[i]=integ0[0]+0.5/c_H0
        if i<10: print("sf=",i, s[i], "H0t=",H0t[i]*c_H0, "LH=",LH[i]*c_H0)
i=0
print("s[]=",i,s[i])
print("H0t[]=",i,H0t[i]*c_H0)
print("LH[]=",i,LH[i]*c_H0)
#
plt.plot(LH*c_H0, H0t*c_H0, color = "green")
plt.title("Space Expansion")
plt.ylabel("Scale factor: handred million year")
plt.xlabel("Space: handred million year *c")

ps=0.435
for i in range(500):
        t0=H0t[i]
        s0=s[i]
#        print("s0=",i,s0)
        if i==0:
            ratio[i]=1
            star0[i]=ps*LH[i]
            star1[i]=ps*LH[i]
            print("age=",i, t0, t0*c_H0,"r=",star0[i]*c_H0,star1[i]*c_H0)                
        else:
            star0[i]=ps*LH[i]
#            print("star0=",star0[i])
            dt=H0t[i]-H0t[i-1]
#            print("dt=",dt)
            v[i]=(LH[i]-LH[i-1])/dt*ps #
            vv=y1(s0)*ps#
            if i==10: print("v=",v[i],vv)           
            if i==1:
                print("ratio[i-1]",ratio[i-1])
                print("star[i-1]=",i, star0[i-1]*c_H0,star1[i-1]*c_H0)
            star1[i]=star1[i-1]+v[i]*ratio[i-1]*dt-1.0*dt
#            print("star1=",star1[i])
            ratio[i]=star1[i]/star0[i]
#            print("ratio=",ratio[i])
plt.plot(star1*c_H0, H0t*c_H0, color = "red")
plt.plot(star0*c_H0, H0t*c_H0, color = "blue")
plt.show()

　計算結果を図5に示す．ps=0.47は，現在（宇宙年齢138億年）で地球に到達するように調整した結果である．すなわち，現在観測している背景放射は宇宙年齢38万年のとき0.5*0.47=0.235億年離れた地点の光ということになる．

　　　　　　　　　　図5 宇宙背景放射軌道計算結果

5. おわりに

　138億年前の光がどのように我々の地球に届くのかを標準宇宙膨張モデルを用いて計算した．本結果を使えば遠方の銀河，星からの軌道計算が可能である．図6に129億年前の恒星エアレンデル(参考文献(4))の軌道計算結果を参考に示す．

　　　　　　図6 129億年前の恒星エアレンデルの軌道計算結果

参考文献
(1) ビックバンモデルを正しく理解する（ひの市民大学資料）
(2) 標準的宇宙モデルは一様等方宇宙に対応する次の線素： ds2
(3) 宇宙晴れ上がりの光は超光速で膨張する宇宙空間をどのようにして地球に到達するか
(4) 史上最も遠い恒星を観測、今は地球から２８０億光年の距離