Hubble (1929) のデータ
Hubble (1929) の Table 1 のデータから Hubble 定数を推定する。
以下のような内容の hubble.dat を用意する。
1列目 2列目 3列目
name r(Mpc) v(km/s)
"S.Mag." 0.032 +170
"L.Mag." 0.034 +290
"6822" 0.214 -130
"598" 0.263 -70
"221" 0.275 -185
"224" 0.275 -220
"5457" 0.45 +200
"4736" 0.5 +290
"5194" 0.5 +270
"4449" 0.63 +200
"4214" 0.8 +300
"3031" 0.9 -30
"3627" 0.9 +650
"4826" 0.9 +150
"5236" 0.9 +500
"1068" 1.0 +920
"5055" 1.1 +450
"7331" 1.1 +500
"4258" 1.4 +500
"4151" 1.7 +960
"4382" 2.0 +500
"4472" 2.0 +850
"4486" 2.0 +800
"4649" 2.0 +1090
銀河までの距離と後退速度のグラフ
In [1]:
%gnuplot inline svg size 640,480 fixed enhanced font ',16'
In [2]:
set xrange [-0.4:2.4]
set yrange [-400:1400]
set xtics 1
set ytics 500
set zeroaxis
set grid
set xlabel "r (Mpc)"
set ylabel "v (km/s)"
set title "Hubble (1929) Table 1."
plot "hubble.dat" using 2:3 pt 7 ps 0.6 notitle
最小二乗法によるフィット
$v = H_0 r$ を仮定して最小二乗法により $H_0$ を推定する。
In [4]:
f(x) = H0 * x
fit f(x) "hubble.dat" using 2:3 via H0
Out[4]:
In [5]:
print sprintf("H0 = %4.1f", H0)
Out[5]:
フィットした直線もグラフに描く
In [6]:
plot "hubble.dat" using 2:3 pt 7 ps 0.6 notitle, \
[0:2.4] f(x) lw 2 title sprintf("H0 = %4.1f", H0)
ハッブル年齢（$H_0$ の逆数）を求める
$\Omega_{\Lambda} = 0$ の場合は，宇宙年齢の上限が $\displaystyle \frac{1}{H_0}$ となる。また，$\Omega_{\rm m} = 1, \Omega_{\Lambda} = 0$ の場合の宇宙年齢は $\displaystyle \frac{2}{3}\frac{1}{H_0}$ となる。
In [8]:
# Mpc を km で
Mpc = 3.0856775814127E19
# km/(s Mpc)
Hinv = 1/H0 * Mpc
# 1/H0 は秒なので「億年」に換算する
okunen = 60*60*24*365.25*1.E8
print sprintf("%3.1f 億年", Hinv/okunen)
Out[8]:
演習問題：現在のデータから $H_0$ を推定する
Hubble (1929) のデータから得られたハッブル定数の値，およびその逆数から得られる宇宙年齢の上限の値は，現在我々が知っている地球年齢（約46億年）よりも短いだけでなく，当時（20世紀前半）の知見でも地球より宇宙の年齢が短くなりそうな危うい状況であった。
その原因は，距離の推定の誤差であった。
現在，我々は以下のリンクからそれぞれの銀河の距離（距離はハッブルの法則以外の手段で求められていると仮定しよう）と赤方偏移を確認できる。これらの値を用いてハッブル定数の推定を行い，ハッブル年齢がいくらになるか調べよ。
In [9]:
reset
set xrange [-1:21]
set yrange [-400:1400]
set ytics 500
set zeroaxis
set grid
set xlabel "r (Mpc)"
set ylabel "v (km/s)"
set title "現在のデータを用いたハッブルダイアグラム"
plot "rvnow.dat" using 2:3 pt 7 ps 0.6 notitle
In [10]:
f(x) = H0 * x
fit f(x) "rvnow.dat" using 2:3 via H0
print sprintf("H0 = %4.1f", H0)
Out[10]:
In [11]:
set key top left
plot "rvnow.dat" using 2:3 pt 7 ps 0.6 notitle, \
[0:21] f(x) lw 2 title sprintf("H0 = %4.1f", H0)
