1. No category

卒業論文
exotic matter の存在しない領域をもつ
ワームホールのシャドー
氏名：森 友紀
名古屋大学理学研究科
素粒子宇宙物理学専攻
宇宙論研究室
2014 年 3 月 14 日
概要
巨視的なワームホールは、誕生したとしても自己重力によってすぐに崩壊してしまう。この崩壊を支え、
ワームホールを維持する為に必要な物質を exotic matter と呼び、通常の物質では考えられない負のエネル
ギーを持っている。Morris ＝ Thorne が提唱した通常のワームホールでは、この exotic matter がワームホー
ルスロート部分の周りを満たしており、必ず exotic matter を通過しないとワームホールを通り抜けることが
できない。exotic matter を通過することは非常に不安定で、すぐに崩壊したり、ブラックホールになったり
するのである。しかし、Edward Teo が、回転するより一般的なワームホールを考え、計量に非球対称な成分
を持たせることによって、exotic matter を通過することなくワームホールを通り抜けられるということが分
かっている。
本研究では Teo の非球対称なワームホールモデルを用いて、光子の軌道を計算し、振る舞いを見て、実際に
観測者からどのように見えるかを議論した。また、exotic matter に出会わずに通過できる軌道があることが
分かった。
目次
目次
目次
第 1 章 導入
3
1.1
メトリック . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.2
測地線方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.3
反変ベクトル・共変ベクトル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.4
共変微分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
1.5
曲率 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
1.6
完全流体の場合のエネルギー・運動量テンソル . . . . . . . . . . . . . . . . . . . . . . . . . .
7
1.7
アインシュタイン方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
1.8
ハミルトンヤコビ方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
第 2 章 ブラックホール
11
シュバルツシルトブラックホール . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.1.1 事象の地平面 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.1.2 特異点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
カーブラックホール . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.2.1 事象の地平面 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.2.2 特異点リング . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.2.3 静止限界 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.2.4 エルゴ領域 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.1
2.2
第 3 章 ブラックホール周辺の光子の軌道計算
13
3.1
数値計算方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
3.2
ブラックホール周辺の光子の軌道計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
3.2.1 カー解で測地線方程式を解く . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
3.2.2 初期条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
3.2.3 不安定軌道 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
3.3
数値計算の結果
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
第 4 章 ワームホール
20
23
4.1
ワームホールと exotic matter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
4.2
Morris=Thorne モデルのワームホール . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
4.3
Edward Teo モデルのワームホール . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
第 5 章 ワームホール周辺の光子の軌道計算
26
数値計算方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
5.1.1 Edward Teo モデルのワームホール解で測地線方程式を解く . . . . . . . . . . . . . . . . .
26
5.1.2 初期条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
5.1.3 不安定軌道 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
5.1
1
目次
5.2
目次
数値計算の結果
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
第 6 章 まとめ
33
38
2
第1章
第1章
導入
導入
1.1 メトリック
一般相対性理論は、時間と空間によって作られる４次元時空上に展開される。時間座標と空間座標をまと
めて、
x0 = ct, x1 = x, x2 = y, x3 = z
(1)
と書ける。xµ (µ = 0, 1, 2, 3) とし、これを４元位置ベクトルと呼ぶ。平坦な４次元ミンコフスキー空間におけ
る４元位置ベクトルの微小変化を xµ *xµ + dxµ とすると、微小距離を
ds2 = −(dx0 )2 + (dx1 )2 + (dx2 )2 + (dx3 )2
(2)
と定義し、この微小距離を不変間隔と呼ぶ。また、これを ηµν という行列を用いると、
ds2 = ηµν dxµ dxν
(3)
と書ける。この、ηµν をメトリックと呼ぶ。このときは空間の座標系に直交座標を用いたが、より一般的な座
標系を用いると、不変距離は
ds2 = gµν dxµ dxν
(4)
と書ける。また、ds2 > 0 のとき空間的、ds2 < 0 のとき時間的、ds2 = 0 が光の場合（ヌル）として、観測す
るものを光子か質量をもった粒子かを決めている。
1.2 測地線方程式
測地線とは、重力場が存在する場合に、質点や光が重力の影響を受けてどのような経路を取るのかを表す式
である。２点間の距離が空間的である場合、不変距離は、
ds =
√
gµν dxµ dxν
(5)
を積分して得られる。まず、2 点 P,Q の隔たりが空間的である場合について考える。パラメーターを u を用
いると、２点を結ぶ曲線の長さ s は、
∫
∫
Q
s=
Q
ds =
P
P
ds
du =
du
∫
Q
√
gµν
P
dxµ dxν
du
du du
(6)
である。測地線とは、この２点間の長さ s が最小になるような曲線である。したがって、この距離が最短とな
るには最小作用の原理を用いればよい。x˙ = dx/du とすると、ラグラジアンは、
L=
√
gµν x˙µ x˙ν
(7)
と書けるので、最短距離になる条件は、L がオイラー・ラグランジュ方程式を満たすことである。オイラー・
ラグランジュ方程式は、
∂L
d ∂L
−
=0
du ∂xµ
∂xµ
3
(8)
1.2
測地線方程式
第1章
導入
であり、これを計算する。まず第１項は、
√
∂L
∂
gµν x˙µ x˙ν )
=
(
∂xµ
∂ x˙λ
(
)
1
∂gµν µ ν
∂ x˙µ κ
∂ x˙ν
µ
˙
˙
˙
˙
= √
x x + gµν
x + gµν x
˙λ
∂ x˙λ
∂ x˙λ
2 gαβ x˙α x˙β ∂ x
=√
gλν x˙ν
(9)
gαβ x˙α x˙β
であるので、したがって、
d
du
(
∂L
∂ x˙λ
)
(
1
gλν
˙ x˙ν (gαβ x˙α x˙β ) + gλν x¨ν (gαβ x˙α x˙β )
(gαβ x˙α x˙β )3/2
)
1
˙ x˙α x˙β + gαβ x¨α x˙β + gαβ x˙α x¨β )
− gλν x˙ν (gαβ
2
]
[
1
1
=
gλσ x¨σ + g σα (∂ρ gαν + ∂ν gαρ )x˙ρ x˙ν
2
(gαβ x˙α x˙β )1/2
]
[
1
1 κσ
˙
κ
β
γ
ν
α
¨
˙
˙
˙
−
gλν x × 2gακ x x + g (∂γ gσβ + ∂β gσγ − ∂σ gβγ )x x
2
2(gαβ x˙α x˙β )3/2
=
(10)
となる。ここで、∂λ gµν = ∂gµν /∂xλ とした。次に、第２項は、
√
∂L
∂
=
(
gµν x˙µ x˙ν )
∂xλ
∂xλ
1
= √
(∂λ gµν )x˙µ x˙ν
2 gαβ x˙α x˙β
1
1
=√
gλσ × g σα (∂α gρν )x˙ρ x˙ν
2
gαβ x˙α x˙β
(11)
である。したがって、式 (8) は、
[
]
d ∂L
∂L
1
1 σα
σ
ν
ρ
¨
˙
˙
− µ =
gλσ x + g (∂ρ gαν + ∂ν gαρ − ∂α gρν )x x
du ∂xµ
∂x
2
(gαβ x˙α x˙β )1/2
[
]
1
1
−
gλν x˙ν × 2gακ x˙α x¨κ + g κσ (∂γ gσβ + ∂β gσγ − ∂σ gβγ )x˙β x˙γ
2
2(gαβ x˙α x˙β )3/2
=0
(12)
である。これが成立するには、
x¨µ + Γµνλ x˙ν x˙λ = 0
(13)
1
Γµνλ ≡ g µκ (∂λ gκν + ∂ν gλκ − ∂κ g λν )
2
(14)
を満たせばよい。ここで、
µ
と定義し、Γνλ をクリストフェル記号と呼ぶ。また式 (14) を測地線方程式と呼ぶ。今、２点間の隔たりが空
間的で微分は不変間隔 s であるので、測地線方程式は、
ν
λ
d2 xν
µ dx dx
+
Γ
=0
νλ
ds2
ds ds
4
(15)
1.3
反変ベクトル・共変ベクトル
第1章
導入
である。また、２点間の隔たりが時間的である場合は固定時間 τ での微分に置き換わり、
d2 xν
dxν dxλ
+ Γµνλ
=0
2
dτ
dτ dτ
(16)
と書ける。さらに、ヌルの場合は、アフィンパラメーター λ を用いて、
d2 xν
dxν dxλ
=0
+ Γµνλ
2
dλ
dλ dλ
(17)
P µ Pµ = gµν P µ P ν = 0
(18)
と書ける。ただしここで、
を同時に満足しなければならない。これを null 条件と呼ぶ。
1.3 反変ベクトル・共変ベクトル
0
座標 xµ について xµ *x µ (x0 , x1 , x2 , x3 ) と変換することを一般座標変換という。すべての物理量が、一般
座標変換に対して同じ変換性を持っている訳ではないが、それぞれ、スカラーやベクトルという量は一般座標
変換に対する変換性を持っている。同じ点で、一般座標変換によって値が変化しない量をスカラーと呼ぶ。あ
0
る点の座標が２つの座標系によって xµ , x µ とかけるとすると、
0
φ0 (x µ ) = φ(xµ )
(19)
となる φ がスカラーである。座標上のどの点でもこの関係が成り立っている場合をスカラー場と呼ぶ。
次にベクトルを定義する。一般座標変換のもとで
∂x0µ ν
A (x)
∂xν
Aµ (x)*A0µ (x0 ) =
(20)
のように変換するベクトル Aµ (x) を反変ベクトルといい、
Bµ (x)*Bµ0 (x0 ) =
∂xν
Bν (x)
∂x0µ
(21)
のように変換するベクトル Bµ (x) を共変ベクトルという。
1.4 共変微分
まず、スカラー場 φ の場合について考える。位置 xµ でのスカラー場の大きさを φ(x) とし、この位置から
微小に離れた位置 xµ + dxµ でのスカラー場の大きさは φ(x + dx) とする。 このとき、φ(x)*φ(x + dx) と
変化させると、2 つのスカラー量の差 δφ は、
δφ = φ(x + dx) −φ(x)
(22)
のようになる。したがって、スカラー場の微分を、
(∂µ φ)hµ ≡ lim
ε→0
φ(x + εh) −φ(x)
ε
(23)
とすると、スカラー量の差は、
δφ = (∂µ φ)dxµ + O(dxµ )
5
(24)
1.5
曲率
第1章
導入
のようにかける。
次に、ベクトル場の微分について考える。V µ (x)*V µ (x + εh) と平行移動を εhλ に沿って行うとき、平行
移動したものは、
V µ (x) −ε hν Γµνκ V κ (x)
(25)
とかける。したがってベクトル場の微分は、
V µ (x + εh) − [V µ (x) − εhν Γµνκ V κ (x)]
= hµ (∂ν V µ + Γµνκ V κ )
ε→0
ε
= hµ ∇ν V µ
lim
(26)
である。これを改めて書き直すと、
∇ν V µ ≡∂ν V µ + Γµνκ V κ
(27)
である。これを反変ベクトル V µ の共変微分と呼ぶ。
また、共変ベクトル Vµ の共変微分について考える。共変ベクトルの平行移動は、
Vµ (x) + εhν Γκνµ Vκ (x)
(28)
と書けるので、共変ベクトルの共変微分は、
∇ν Vµ ≡∂ν V µ − Γκνµ Vκ
(29)
と定義できる。
1.5 曲率
スカラー φ の 2 階共変微分について考えると、
∇ν ∇µ φ = ∂ν (∇µ φ) −Γκνµ ∇κ φ
= ∂µ ∂ν φ−Γκνµ ∂κ φ
(30)
である。これより、添字 µ, ν を入れ替えても不変であり、
∇ν ∇µ φ = ∇µ ∇ν φ
(31)
であることが分かる。
次にベクトル aµ の２階共変微分について考えると、
∇ν ∇κ aµ = ∂ν (∇κ aν ) − Γσνκ (∇σ aµ ) − Γσνµ (∇κ aσ )
= ∂ν (∂κ aµ − Γλκµ aλ ) − Γσνκ (∂σ aµ − Γλσν aλ ) − Γσνµ (∂κ aσ − Γλκσ aλ )
= ∂ν ∂κ aµ − (Γλκµ ∂ν ∂λ + Γσνµ ∂κ aσ − Γσνκ (∂σ aν + Γλσµ aλ ) − aλ (∂ν Γλκµ − Γσνµ Γλκσ )
(32)
である。これより、ベクトルの 2 階共変微分は微分する順序によって結果が変わるということが分かる。つま
り、ある点 P から点 Q へベクトルを平行移動する際に、κ 方向に移動してから ν 方向に移動した結果と ν 方
向に移動してから κ 方向に移動した結果では異なる。その差は、
∇ν ∇κ aµ − ∇κ ∇ν aµ = aλ (∂κ Γλσµ − ∂ν Γλκµ + Γσνµ Γλκσ − Γσκµ Γλνσ )
λ
≡ aλ Rµκν
(33)
6
1.6
完全流体の場合のエネルギー・運動量テンソル
第1章
導入
である。改めて書き直すと、
µ
Rνλκ
= ∂λ Γµνκ − ∂κ Γµνλ + Γµσλ Γσνκ − Γµσκ Γσνλ
(34)
であり、これをリーマン曲率テンソルと呼ぶ。リーマン曲率テンソルこそ時空の曲がりを表す指標となる。
リーマン曲率テンソルの２つの足の縮約をとることで、
κ
Rµν ≡Rµκν
= g κσ Rσµκν
(35)
となり、これをリッチテンソルと呼ぶ。リーマン曲率テンソルにメトリックを作用させることで、共変なリー
σ
マン曲率テンソルは Rµνλκ = gµσ Rνλκ
と表され、
∇λ Rµνκσ + ∇σ Rµνλκ + ∇κ Rµνσλ = 0
(36)
という共変微分の関係式も成立する。これをビアンキ恒等式と呼ぶ。ビアンキ恒等式の縮約をとると、
(
)
1 µν
µν
∇ν R − g R = 0
2
(37)
である。また、Gµν = Rµν − 21 g µν R として、これをアインシュタインテンソルと呼ぶ。
1.6 完全流体の場合のエネルギー・運動量テンソル
Newton 力学の場合、ポテンシャル φ を考え、質量がない真空中では、
52 φ = 0
(38)
と書ける。このような形の式 をラプラス方程式という。質量が密度 ρ(x) で分布している場合の重力場の方程
式はラプラス方程式より、
52 φ(x) = 4πGρ(x)
(39)
となり、このような方程式をポアソン方程式という。一般相対論において、このポアソン方程式に相当する式
を求める。まず、ポアソン方程式の右辺について考える。密度も圧力と組み合わされて、４元化することが可
能である。これは、流体の場合、xν 一定面を横切る４元運動量 pµ の流れとして表され、この T µν を成分と
するテンソルをエネルギー・運動量テンソルと呼ぶ。まず、この流体に対して局所ローレンツ系を考える。ま
た、流体の単位当たりに含まれる粒子数を n = n(xµ )、粒子１つ当たりの質量を m とする。このとき単位体
積当たりの質量、つまり質量密度 ρ = ρ(xµ ) は ρ = nm と表される。今、この流体の一部分に対して局所ロー
レンツ系をとったことにより、考えている部分での流体の４元速度は、
uµ =
dxµ
dt xµ
=
= γ x˙µ
dτ
dτ dt
(40)
と表される。ここで、γ = dt/dτ とした。流体が光速に比べて小さいとすると、γ ' 1 であるから、u0 'c, ui 'ni
とかける。この条件の元では、T 00 は、時間一定面を横切る運動量の０成分の流れである。時間一定面を横切
る運動量の０成分は p0 = mc、流れは nu0 = nc であるから、p0 ×nc = mnc2 = ρc2 となる。したがって、
T 00 はエネルギー密度を表している。
次に T i0 は、時間一定面を横切る運動量の i 成分の流れなので、pi ×nu0 = mv i ×nc = cρv i となり、これは
単位体積当たりの運動量密度に光速度 c をかけたものである。また T 0i は、xi 一定面を横切る運動量 p0 の流
7
1.7
アインシュタイン方程式
第1章
導入
れなので、p0 ×nui = mc×nv i = cρv i となり、T i0 と等しい。物理的な意味はエネルギー流速である。
そして T ij は、xj 面を横切る運動量 pi の流れなので、pi ×nuj = mv i ×nv j = ρv i v j である。i = j のとき、
流体を構成する粒子は v i だけ移動する。単位体積当たりの粒子数は n だったので、nv i は単位時間に xi 一定
面の単位面積を通過した粒子の総数を表す。よって、pi ×nv i は単位時間に単位面積当たりの運動量の総変化
量を表している、したがってこれは、圧力である。一方、i6=j の場合は応力を表している。
圧力を持つ完全流体の場合、局所ローレンツ静止系 T¯µν = diag(ρc2 , p, p, p) となる。テンソルの線形性を考
えると、T¯ 00 = ρc2 の部分は圧力０の流体と同じく、ρuµ uν によって表すことができる。一方、圧力の部分
は、局所ローレンツ静止系の場合に、ui = 0 となってしまうので、４元速度を使って表すことができない。
局所ローレンツ静止系の場合には、g µν = η µν = diag(−1, 1, 1, 1) である。pg µν は局所ローレンツ静止系で、
(−p, p, p, p) となるので、００成分を除けばうまくいっている。T¯µν = ρuµ uν + pη µν としたいが、余分な項
−p が現れてしまう。しかし、この局所ローレンツ静止系では、uµ uν は００成分にしか影響を与えないので、
係数は自由にやりくりして良い。新たに付け加えるべきは +p = ku0 u0 より k = p/c2 とすればよい。した
がって、
(
)
p
T¯µν = ρ + 2 uν uµ + pg µν
c
(41)
を得る。次に、T µν についてエネルギーと運動量の保存について考える。局所ローレンツ系をとると、エネル
ギー保存則は、∂µ T 0µ = 0 である。また、運動量保存則は、∂µ T iµ = 0 と表される。以上の局所ローレンツ系
での保存則を、一般の座標系に拡張すると、
∇µ T µν = 0
(42)
と書ける。
1.7 アインシュタイン方程式
式 (42) とアインシュタインテンソルを関係づけて、重力場の方程式を
Gµν = κT µν
(43)
とし、κ を求める。重力場が弱いときは、メトリックのミンコフスキー時空の ηµν からずれはわずかであると
考えられる。したがって、ミンコフスキー時空からのずれを hµν とすると、メトリックは、
gµν = ηµν + hµν
(44)
と表される。このメトリックは、時間的にゆっくりと変化しているので、近似的に ∂0 gµν = 0 つまり、
∂0 hµν = 0 である。弱い重力場中での質点の運動は、測地線方程式
d2 xµ
dxν dxλ
+ Γµνλ
=0
2
dτ
dτ dτ
(45)
によって記述できる。それがニュートン運動方程式と一致するということから、メトリックのずれ ηµν と重力
ポテンシャルの関係が求まる。それを重力場の方程式に代入することによってポアソン方程式と比較し κ を
決める。クリストフェル記号は、近似を用いると、
Γµνλ =
1
1 νκ
g (∂λ gκν + ∂ν gκλ − ∂κ gλν ) = η µκ (∂λ hκν + ∂ν hκλ − ∂κ hλν )
2
2
8
(46)
1.7
アインシュタイン方程式
第1章
導入
と書ける。したがって、
1 0κ
η (∂0 hκ0 + ∂0 hκ0 − ∂κ h00 )
2
1
= − ∂0 h00 = 0
2
Γ000 =
(47)
である。次に、
1 iκ
η (∂0 hκ0 + ∂0 hκ0 + ∂0 h00 )
2
1
= − ∂ i h00
2
Γi00 =
(48)
である。以上より、測地線方程式は、
d2 x0
=0
dτ 2
(49)
d2 xi
1 i
dx0 dx0
−
∂
h
=0
00
dτ 2
2
dτ dτ
(50)
と表される。 ここで、x0 = cτ であり、dx0 /dτ = c なので、
d2 xi
1
= c2 ∂ i h00
dτ 2
2
(51)
d2 xi
= −∂ i φ
dt2
(52)
である。これとニュートンの運動方程式
を比較すると、
h00 = −
2φ
c2
(53)
であればよい。ここまでに導出した式を重力場の方程式に代入すると、
1 2
1
∇ φ = κρc2
c2
2
(54)
となる。これとポアソン方程式を比較することで、
κ=
8πG
c4
(55)
を得ることができる。したがって、重力場の基礎方程式は、
8πG
1
Rµν − Rg µν = 4 T µν
2
c
であることが分かった。これをアインシュタイン方程式と呼ぶ。
9
(56)
1.8
ハミルトンヤコビ方程式
第1章
導入
1.8 ハミルトンヤコビ方程式
質点系の自由度を n、一般化座標を q1 , ...qn 、ラグランジアンを L とする。このとき、
pi =
H=
∂L
∂ q˙i
n
∑
(57)
pi q˙i − L
(58)
i=1
と書け、pi を qi に共役な一般運動量といい、H をハミルトニアンという。このとき、
∂H
dqi
=
dt
∂pi
dpi
∂H
=−
dt
∂qi
(59)
(60)
が成り立ち、これをハミルトンの正準方程式という。q1 , ..., qn , p1 , ..., pn を変数とする関数 Qi , Pi を考えると、
Qi = Qi (q1 , ..., qn , p1 , ..., pn )
Pi = Pi (q1 , ..., qn , p1 , ..., pn )
(61)
(62)
である。このとき、(q, p) * (Q, P ) を正準変換という。正準変換によって変換された後のハミルトニアン K
が K = 0 になるような特別な正準変換を考える。ハミルトニアンがゼロなので、新しい正凖変数は時間変化
しないので、
∂K
Q˙ i =
=0
∂Pi
∂K
P˙i = −
=0
∂Qi
(63)
(64)
である。この正準変換の母関数を S(q, P, t) とすると、
H(q, p, t) +
と書ける。ここで、正準運動量 p は pi =
∂S
∂t
∂S
=0
∂t
(65)
と表されるので、
)
∂S
∂S
,t +
=0
H q,
∂q
∂t
(
と書かれる。これをハミルトンヤコビ方程式と呼ぶ。
10
(66)
第2章
第2章
ブラックホール
ブラックホール
2.1 シュバルツシルトブラックホール
回転・電荷なしのブラックホールモデルである。シュバルツシルト計量とは、球対称真空、かつ静的で時間
変化が無いという仮定の下でアインシュタイン方程式を解いたときの解であり、
(
)
2GM
1
dr2 + r2 (dθ2 + sin2 dφ2 )
ds = − 1 −
(cdt)2 +
2
rc
1 − 2GM
2
rc
2
(67)
と書ける。ここで、G は万有引力定数、M はブラッックホールの質量、c は光の速度である。
2.1.1 事象の地平面
シュバルツシルト計量の grr 成分が発散する場所を事象の地平面と呼び、この事象の地平面の半径を特に、
シュバルツシルト半径と呼ばれ、
rg =
2GM
c2
(68)
と表される。この半径の内側では、光でさえも脱出できなくなるほど時空の曲がりが大きくなっている。
2.1.2 特異点
中心（r = 0）は特異点となっていて、ブラックホールの強い重力場をつくる物質が、そこに集中している
場所である。事象の地平面の内部ではいったんその内側にはいると、いくら外向きに加速しようが、 絶えず
内向きにどんどん落ちていって、その多くの軌道は最後は特異点に到達する。
2.2 カーブラックホール
回転あり、電荷なしのブラックホールモデルである。カー計量とは、φ 軸対称真空、φ 方向に回転し、かつ
静的で時間変化が無いという仮定の下でアインシュタイン方程式を解いたときの解であり、
(
)
2M r
4aM r sin2 θ
ρ2
2M ra sin2 θ
ds2 = − 1 − 2
(cdt)2 −
dtdφ + dr2 + ρ2 dθ2 + (r2 + a2 +
)sin2 θdφ2
2
ρ
ρ
∆
ρ2
(69)
∆ = r2 − 2M r + a2 , ρ2 = r2 + a2 cos2 θ
(70)
と書ける。ここで、a は回転の大きさ、M はブラッックホールの質量、c は光の速度である。
2.2.1 事象の地平面
カー計量の grr 成分が発散する場所、つまり ∆ = 0、さらに、
√
r± = M ± M 2 − a2
(71)
によって与えられる。r+ は外部地平面、r− は内部地平面と呼ばれ、r+ を事象の地平面とする。a = 0 のと
き、シュバルツ シルト半径に一致する。
11
2.2
カーブラックホール
第2章
ブラックホール
2.2.2 特異点リング
ρ2 = 0 となるところ、つまり r = 0, θ = π/2 は特異点で、半径 r = 0 のリング状の特異点になっている。
これは、適当な座標のもとでは有限な大きさのリングになっている。
2.2.3 静止限界
メトリックの時間と空間部分の符号が入れ替わる境界を静止限界と呼ぶ。事象の地平面と回転軸（θ = 0, π)
で接している。静止限界より外側では、原理的には逆噴射などによって回転方向の力に逆らうことができる
が、静止限界よる内側では、必ずカーブラックホールの回転によって引きずられてしまうのである。つまり、
観測者が決して静止できなくなる限界である。
2.2.4 エルゴ領域
事象の地平面と静止限界の間の領域をエルゴ領域と呼び、gtt > 0、つまり、
M+
√
√
M 2 − a2 < r < M + M 2 − a2 cos2 θ
(72)
の領域で与えられる。観測者が赤道上面上 (θ = π/2) にいる場合は、a が大きくなれば、事象の地平面の半径
は小さくなるのでエルゴ領域は大きくなり、反対に、a が小さくなれば、事象の地平面の半径は大きくなるの
でエルゴ領域は小さくなる。a = 0 のとき、事象の地平面と静止限界は一致し、エルゴ領域はなくなる。この
領域内では、
『慣性の引きずり』と呼ばれる現象が見られ、ブラックホールの回転と同じ方向に引きずられる。
エルゴ領域内では、空間に存在する光子が光速以上で引きずられるので、エルゴ領域の外で静止した観測者か
ら見て静止することはできない。
12
第3章
第3章
ブラックホール周辺の光子の軌道計算
ブラックホール周辺の光子の軌道計算
3.1 数値計算方法
今回はより計算の精度を上げる為に、４次のルンゲクッタ法を用いた。初期値を (ti , yi , fi )、変化分をそれ
ぞれ h, k, l として、変化後の位置は (ti + h, yi + k, fi + l) と表される。
dy
= f (t, y)
dt
df
= g (t, f (t, y))
dt
(73)
(74)
という連立常微分方程式の解を、４次のルンゲクッタ法によって求める。
K0 = hf (ti , yi )
L0 = hg(ti , fi )
h
K0
, yi +
)
2
2
h
L0
= hg(ti + , fi +
)
2
2
h
K1
= hf (ti + , yi +
)
2
2
L1
h
)
= hf (ti + , yi +
2
2
= hf (ti + h, yi + K2 )
K1 = hf (ti +
L1
K2
L2
K3
(75)
L3 = hg(ti + h, fi + L2 )
これらの K, L を用いると、
1
(K0 + 2K1 + 2K2 + K3 )
6
1
L = (L0 + 2L1 + 2L2 + L3 )
6
K=
(76)
(77)
となる。
3.2 ブラックホール周辺の光子の軌道計算
3.2.1 カー解で測地線方程式を解く
カー計量のラグラジアンは、
2L = −(1 −
ρ2 2
2M ra sin2 θ
2M r 2 4aM r sin2 θ
2
2
2
2
)dt
−
dtdφ
+
dr
+
ρ
dθ
+
(r
+
a
+
)sin2 θdφ2
ρ2
ρ2
∆
ρ2
∆ = r2 − 2M r + a2 , ρ2 = r2 + a2 cos2 θ
(78)
(79)
であり、メトリックは、



gµν = 

−(1 −
2M r
ρ2 )
0
0
2
− 2aM ρr2sin
θ
0
ρ2
∆
0
0
2
− 2aM ρr2sin
0
0
ρ2
0
13

θ




0
0
(r2 + a2 +
2M ra sin2 θ
)sin2
ρ2
θ
(80)
ブラックホール周辺の光子の軌道計算
3.2
第3章



g µν = 

2
− ρΣ2 ∆
0
0
2M ra
− ρ2 ∆
0
∆
ρ2
0
0

1
ρ2
ra
− 2M
ρ2 ∆
0
0
0
∆−a sin θ
ρ2 ∆ sin2 θ
0
0
2
ブラックホール周辺の光子の軌道計算




(81)
2
である。だだし、Σ2 = (r2 + a2 )2 − a2 ∆ sin2 θ とする。これより、ハミルトンヤコビの方程式に代入し、ア
フィンパラメーター λ を導入すると、
2
∂S
∂S ∂S
= g µν µ ν
∂λ
∂x ∂x
( )2
( )2
( )2
( )2
4M ra ∂S ∂S
∆ − a2 sin2 θ ∂S
∆ ∂S
1 ∂S
∂S
Σ2
+ 2
− 2
−
−
= 2
ρ ∆ ∂t
ρ ∆ ∂t ∂φ
∂φ
ρ2 ∂r
ρ2 ∂θ
ρ ∆ sin2 θ
[
]2
[
]2
( )2
( )2
1
∂S
1
∂S
∆ ∂S
1 ∂S
∂S
2
2
2 2 ∂S
= 2
(r + a )
−a
− 2 2
+
− 2
− 2
(a sin θ)
ρ ∆
∂t
∂φ
∂t
∂φ
ρ
∂r
ρ
∂θ
ρ sin θ
(82)
である。ここで、δ1 を粒子の質量として、S = 12 δ1 λ − Et + Lz φ + Sr (r) + Sθ (θ) とすると、
]2
]2
1 [
1 [ 2
−(a sin2 θ)E + Lz − ∆
(r + a2 )E − aLz −
ρ δ1 =
∆
sin2 θ
(
2
(
∆
dSr
dr
)2
dSr
dr
)2
(
−
1 2
[(r + a2 )E − aLz ]2 + (Lz − aE)2 + δ1 r2
∆
(
)2
dSθ
+
+ (L2z csc2 θ) cos2 θ + δ1 a2 cos2 θ = 0
dθ
dSθ
dθ
)2
(83)
+
(84)
となる。L を分離定数とすると、
(
∆2
dSr
dr
)2
= [(r2 + a2 )E − aLz ]2 − ∆[(Lz − aE)2 + δ1 r2 + L]≡R
(
dSθ 2
dθ
(85)
)
= L − (L2z csc2 θ − a2 E 2 + δ1 a2 ) cos2 θ≡Θ
(86)
と表される。したがって、
1
S = δ1 λ − Et + Lz φ +
2
∫
r
√
∫
R
dr +
∆
θ
√
Θdθ
(87)
と改めて書き直すことができる。
∂S
1
=
∂L
2
∫
∫
1 θ 1 ∂Θ
1 ∂R
√
√
dr +
dθ = 0
2
∆ R ∂L
Θ ∂L
∫ r
∫ θ
1
1
√ dr =
√ dθ
∴
R
Θ
r
(88)
(89)
なので、
dr
=
dθ
√
14
R
Θ
(90)
3.2
ブラックホール周辺の光子の軌道計算
第3章
ブラックホール周辺の光子の軌道計算
と求まる。
∫
1 ∂R
1 θ 1 ∂Θ
√
√
dr +
dθ = 0
2
∆ R ∂δ1
Θ ∂δ1
∫ r 2
∫ θ
r
cos2 θ
√ dr +
∴λ=
a2 √ dθ
R
Θ
∂S
1
1
= λ+
∂δ1
2
2
∫
r
(91)
(92)
なので、
(
ρ4
(
4
ρ
dr
dλ
dθ
dλ
)2
=R
(93)
=Θ
(94)
)2
と求まる。
∫
∫
1 ∂R
1 θ 1 ∂Θ
√
√
dr +
dθ = 0
2
∆ R ∂E
Θ ∂E
∫ r
r2
√ r[r2 E − a(Lz − E)]dr
∴ t = Eλ + 2M
∆ R
∂S
1
= −t +
∂E
2
r
(95)
(96)
なので、
(
2
ρ
dt
dλ
)
=
)
1 ( 2
Σ E − 2aM rLz
∆
(97)
と求まる。
∂S
1
= −φ +
∂Lz
2
∫
r
∴φ=a
∫
r
1 ∂R
1
√
dr +
2
∆ R ∂Lz
r2
√ [(r2 + a2 )E − aLz ]dr +
∆ R
∫
∫
θ
θ
1 ∂Θ
√
dθ = 0
Θ ∂Lz
1
a2 √ (Lz csc2 θ − aE)dθ
Θ
(98)
(99)
なので、
(
ρ2
dφ
dλ
)
=
1
∆
(
2aM rE + (ρ2 − 2M r)Lz
1
sin2 θ
)
cos2 θ
(100)
と求まる。以上より、測地線方程式が求まったので、これを数値計算で解く。
3.2.2 初期条件
観測者の周辺では局所ローレンツ系であるということから初期条件を求める。測地線方程式より、４元運動
量は、
Pt =
1 (
ρ2 ∆
Σ2 E − 2aM rLz
)
√
R
r
P = 2
ρ
√
Θ
Pθ = 2
ρ
(
)
1
1
2aM rE + (ρ2 − 2M r)Lz 2
cos2 θ
Pφ = 2
ρ ∆
sin θ
15
(101)
(102)
(103)
(104)
3.2
ブラックホール周辺の光子の軌道計算
第3章
ブラックホール周辺の光子の軌道計算
である。したがって、Pµ = gµν P ν より、
Pt = −E
√
R
Pr =
∆
√
Pθ = Θ
Pφ = Lz
(105)
(106)
(107)
(108)
と書ける。曲がった時空と局所ローレンツ系とを結びつける四脚場ベクトルは、
eµ(t) =
eµ(r)
=
eµ(θ) =
eµ(φ) =
[ √
[
[
ra
√2M
2
2
Σ2
ρ2 ∆
√
0
ρ ∆Σ
∆
ρ2
0
]
0
0
]
0
]
√1 2 0
ρ
]
√
ρ2
0
2
2
Σ sin θ
0 0
[
0 0
(109)
である。したがって、局所ローレンツ系における４元運動量は、
√
P (t) = −P(t) = −eµ(t) Pµ =
√
P
(r)
= P(r) =
eµ(r) Pµ
√
P (θ) = P(θ) =
1
=
∆
Θ
ρ2
√
P (φ) = P(φ) = Lz
Σ2
2M ra
(E −
Lz )
2
ρ ∆
Σ2
R
ρ2
(110)
(111)
(112)
Σ2
ρ2
sin2 θ
(113)
と書ける。いま、null 測地線方程式に対して δ1 = 0 であり、ξ = Lz /E, η = Q/E 2 という新たなパラメー
ターを導入し、アフィンパラメーターを λnew = Eλ とすると、
R = [(r2 + a2 ) − aξ]2 − ∆[(ξ − a)2 + η]
ξ2
Θ = η − ( 2 − a2 ) cos2 θ
sin θ
(114)
(115)
と書き直すことができる。
さらに、無限遠にいる観測者が見る像の見かけの位置を与える α と β の天球座標に直す。θ0 をブラック
ホールの回転軸と観測者の視線の間の角度とすると、
(
rP (φ)
P(t)
)
ξ
sin θ0
r*∞
√(
)
( (θ) )
ξ2
rP
=
η + a2 cos2 θ0 −
β=
P(t) r*∞
sin2 θ0
α=
=
(116)
(117)
となる。ここで r*∞ のとき、
P (t) = P(t) = 1
16
(118)
3.2
ブラックホール周辺の光子の軌道計算
第3章
ブラックホール周辺の光子の軌道計算
とすると、
√
βP (t)
β
= = ρ2 P θ
r
r
√
αP (t)
α
Σ2
=
= =
sin θP φ
r
r
ρ2
P (θ) =
(119)
P (φ)
(120)
より、
(
θ
P =
(
Pφ =
β
√
r ρ2
α
r sin θ
)
=
)
r*∞
√
ρ2
Σ2
β
r2
(121)
=
r*∞
α
r2 sin θ
(122)
であることが分かる。null 条件 P µ Pµ = 0 より、
√
P
(r)
=− 1−
α2 + β 2
r2
(123)
であるので、以上より、初期条件における測地線方程式が求まった。
3.2.3 不安定軌道
ブラックホールのまわりの光子の軌道には、無限遠から無限遠にぬけていく軌道、ブラックホールに落ちて
いく軌道、そして、ブラックホールの周りをぐるぐると円を描く軌道と３種類存在することが分かっている。
このブラックホールの周りで円を描く軌道を不安定軌道と呼ぶ。測地線方程式が、
(
のように書けるとすると、V = 0,
dV
dr
dr
dλ
)2
+V =0
(124)
2
= 0, ddrV2 < 0 となるとき、つまり、ポテンシャル V が不安定点をとる
ときが不安定軌道になる条件である。この条件を解けば、不安定軌道になるときの衝突パラメーターを求める
ことができる。
具体的に、赤道上面での不安定軌道になるとき、すなわち β = 0 のときの衝突パラメーターを求める。いま、
式 (93) より、測地線方程式は、
(
2
ρ
dr
dλ
)2
= R = [(r2 + a2 ) − aξ]2 − ∆[(ξ − a)2 + η]
(125)
[
]
= r (2M − r)α2 − 4aM α + (r3 + a2 r + 2M a2 ) = −V
(126)
である。θ = π/2 を代入すると、
(
r2
dr
dλ
)2
と書ける。これより、ポテンシャル V と r のグラフは、
17
3.2
ブラックホール周辺の光子の軌道計算
図1
第3章
ブラックホール周辺の光子の軌道計算
a = 0, θ = π/2, β = 0 のときのポテンシャル
√
√
a = 0 のときシュバルツシルトブラックホールである。このとき、不安定軌道になるのは α = 3 3, −3 3
のときであることが分かる。
18
3.2
ブラックホール周辺の光子の軌道計算
図2
第3章
ブラックホール周辺の光子の軌道計算
a = 0.5, θ = π/2, β = 0 のときのポテンシャル
ブラックホールの回転の大きさは a = 0.5。このとき、不安定軌道になるのは α = 4.0, −6.0 のときである
ことが分かる。
また、今回は測地線方程式を求める際に変数分離をすることができたので、理論的に不安定軌道を求めるこ
とができ、
V = R = [(r2 + a2 ) − aξ]2 − ∆[(ξ − a)2 + η] = 0
(127)
∂V
∂R
=
= 4r3 + 2(a2 − ξ 2 − η)r + 2M [η + (ξ − a)2 ] = 0
∂r
∂r
(128)
a2 (r − M )ξ 2 − 2aM (r2 − a2 )ξ − (r2 + a2 )[r(r2 + a2 ) − M (3r2 − a2 )] = 0
(129)
より、
と書ける。これを解くと、
ξ=
M (r2 − a2 )±r∆
a(r − M )
(130)
ξ=
M (r2 − a2 ) + r∆
a(r − M )
(131)
と求まる。ただし、
のとき、
η=−
19
r4
a2
(132)
3.3
数値計算の結果
第3章
ブラックホール周辺の光子の軌道計算
であり、これは不適切である。したがって、
M (r2 − a2 ) − r∆
a(r − M )
(133)
r3 [4M ∆ − r(r − M )2 ]
a2 (r − M )2
(134)
M (r2 − a2 ) − r∆
a sin θ0 (r − M )
(135)
ξ=
η=
である。これと、式 (116)(117) より、
α=
[M (r2 − a2 ) − r∆]2
β = a cos θ0 −
a2 (r − M )2
2
2
2
(
cos θ0
sin θ0
)2
+
r3
[4a2 M − r(r − 3M )2 ]
a2 (r − M )2
(136)
と求まる。これより、理論的に不安定軌道になる衝突パラメーターを求めることができる。
3.3 数値計算の結果
今回の計算では、C = G = M = 1、初期条件 r = 200, φ = 0、として計算を行う。なお、計算コードは巻
末に付録として記載してある。
図3
a = 0, θ = π/2, β = 0 のときの軌道
青：事象の地平面 (半径 2)
20
3.3
数値計算の結果
第3章
ブラックホール周辺の光子の軌道計算
a = 0 のときシュバルツシルトブラックホールである。このとき、事象の地平面と静止限界は一致し、エル
√
ゴ領域は見られない。衝突パラメーターを変化させていくと、無限遠から無限遠に抜けていき、α = 3 3 のと
√
きに円軌道になる。それよりも内側になると、ブラックホールに落ち込む軌道になり、そして再び α = −3 3
のとき円軌道になり、また無限遠から無限遠に抜けていく。つまり、ブラックホールの回転が無いとき衝突パ
ラメーターは対称である。
図4
a = 0.5, θ = π/2, β = 0 のときの軌道
青：事象の地平面 (半径 1.866)、緑：エルゴ領域 (半径 2)
ブラックホールの回転の大きさは a = 0.5。ブラックホールの回転は反時計回り。衝突パラメーターを変化
させていくと、無限遠から無限遠に抜けていき、α = 4.0 付近で円軌道になる。それよりも内側になると、エ
ルゴ領域内でブラックホールの回転と同じ方向に引きずられたのちブラックホールに落ち込む。そして、再び
α = −6.0 付近で円軌道になり、また無限遠から無限遠に抜けていく。ここで、理論値と数値計算の結果に多
少のズレが生じているが、これは、実際には観測者は無限遠にいるとしたが、数値計算の初期条件で r = 200
と設定したためだと考えられる。したがって、r のを大きくすればするほど計算の誤差は少なくなっていくと
思われる。
α と β の変化させていき、不安定軌道になる値をプロットしていくと、ブラックホールのシルエットが描
ける。
21
3.3
数値計算の結果
第3章
図5
ブラックホール周辺の光子の軌道計算
ブラックホールシャドー
a はブラックホールの回転、th は回転軸と観測者の間の角度を表している。この結果は、参考文献 [9] の高
橋労太 (2005) の論文を正しく再現できている。
22
第4章
第4章
ワームホール
ワームホール
4.1 ワームホールと exotic matter
ブラックホールは必ず特異点に落ちていくので
通過することはできないが、ワームホールは誕生
した段階で進行方向に対して地平面を持たず、特
異点も持たないような構造になっているので通
過可能である。このワームホールのトンネル構
造ののどの部分をスロートと呼ぶ。また、巨視的
なワームホールは、誕生したとしても自己重力
によってすぐに崩壊してしまう。この、崩壊を支
え、ワームホールを維持する為に必要な物質を
exotic matter と呼ぶ。
図 6 ワームホールの概念図
通常の物質は正のエネルギー、正の質量をもっているが、exotic matter は「負のエネルギー」
、
「負の質量」
をもった仮想的な物質である。つまり、エネルギー・運動量テンソルのエネルギー密度の項、(t, t) 成分は負
となっている。ヌルベクトルを、
(
µ
k =
とすると、（ただし、e−µ = 1 −
b
r
1
ω
, −e−µ/2 , 0,
N
N
)
(137)
とおいた。
）
Tµν k µ k ν < 0
となる。アインシュタイン方程式 Rµν − 12 Rgµν =
(
8πG
c4 Tµν
(138)
より、
)
1
Rµν − Rgµν k µ k ν = Rµν k µ k ν < 0
2
(139)
である。つまり、exotic matter はスロートで null エネルギー条件
Rµν k µ k ν > 0
を破るということを示している。
23
(140)
4.2
Morris=Thorne モデルのワームホール
第4章
ワームホール
4.2 Morris=Thorne モデルのワームホール
回転のなしの単純化されたワームホールモデルである。Morris=Thorne のワームホールモデルの計量は、
(
b
ds = −N (cdt) + 1 −
r
2
2
)−1
dr2 + r2 K 2 [dθ2 + sin2 θdφ2 ]
2
(141)
N = K = 1, b = 1
(142)
Rµν k µ k ν < 0
(143)
と書ける。この計量を用いて計算すると、
となり、必ずヌルエネルギー条件を破る。つまり、このモデルでは、スロートのまわりが exotic matter で満
たされており、測地線は exotic matter を必ず通過しなければならない通過する際にゆらぎを与えてしまう
と、ゆらぎが増大してしまい、ブラックホールになってしまうような不安定なものである。
4.3 Edward Teo モデルのワームホール
Morris=Thorne のモデルより一般的な回転するワームホールモデルである。Teo のワームホールモデルの
計量は、
(
)−1
b
ds2 = −N 2 (cdt)2 + 1 −
dr2 + r2 K 2 [dθ2 + sin2 θ(dφ − ωdt)2 ]
r
N =K =1+c
(4a cos θ)2
,b = 1
r
(144)
(145)
と書ける。この計量を用いて計算すると、
( )2
( )
dµ 1 d(rK)
dω
sin2 θ
1
dµ
−
−
dr rK dr
dθ
2N 2
4(rK)2 dθ
(
)
(
)
1
d dµ
1
d dN
−
sin
θ
+
sin
θ
2(rK)2 sin θ dθ dθ
(rK)2 N sin θ dθ dθ
(
)
dµ 1 drK
1
d dN
= eµ
+
sin
θ
dr rK dr
(rK)2 N sin θ dθ dθ
b
2(4a)2
1
= 3 −c
(3 cos2 θ − 1)
r K
r (rK)2 N
Rµν k µ k ν = eµ
(146)
である。ここで、
−
b
<0
r3 K
(147)
であるので、式 (146) の第２項目については、角度によって null エネルギー条件を満たし得る。つまり、
Teo のワームホールモデルは、exotic matter が存在しない領域を有しているということが分かる。例えば、
θ = π/2 のときは、
Rµν k µ k ν = −
1
(1 − 2(4a)2 c) > 0
r3
24
(148)
4.3
Edward Teo モデルのワームホール
第4章
ワームホール
を解く。ここで、
(4a)2 c = C
(149)
C > 0.5
(150)
とすると、
のとき、exotic matter が存在しない領域があることが分かった。また、C = 16 のとき、exotic matter の存
在しない領域は、
Rµν k µ k ν = −
b
r3 K
−
32
1
(3 cos2 θ − 1) > 0
r (rK)2 N
(151)
を解くと、
56° < θ < 124°
であることが分かった。
25
(152)
第5章
第5章
ワームホール周辺の光子の軌道計算
ワームホール周辺の光子の軌道計算
5.1 数値計算方法
5.1.1 Edward Teo モデルのワームホール解で測地線方程式を解く
ワームホール計量のラグラジアンは、
(
)−1
b
2L = −N 2 dt2 + 1 −
dr2 + r2 K 2 [dθ2 + sin2 θ(dφ − ωdt)2 ]
r
N =K =1+c
(153)
(4a cos θ)2
2a
, b = 1, ω = 3
r
r
(154)
であり、メトリックは、

gµν
−N 2 + r2 K 2 sin2 θω 2

0
=

0
−r2 K 2 sin2 θω

g µν
− N12
 0
=
 0
− Nω2
0
)−1
(
1 − rb
0
0
( 0 b)
1− r
0
0

−r2 K 2 sin2 θω

0


0
2
2 2
r K sin θ
0
0
r2 K 2
0

− Nω2
0
0
0
0
1
r2 K 2
(155)



(156)
N 2 −r 2 K 2 sin2 θω 2
r 2 K 2 sin2 θN 2
0
である。これより、ハミルトンヤコビの方程式に代入し、アフィンパラメーター λ を導入すると、
2
∂S
∂S ∂S
= g µν µ ν
∂λ
∂x ∂x
( )2
( )2 (
) ( )2
( )2
1
∂S
2ω ∂S ∂S
N 2 − r2 K 2 sin2 θω 2 ∂S
b
∂S
1
∂S
= 2
+ 2
−
−
1
−
−
2
2
2
2
2
2
N
∂t
N ∂t ∂φ
∂φ
r
∂r
r K
∂θ
r K sin θN
(157)
と書ける。ここで、δ1 を粒子の質量として、S = 21 δ1 λ − Et + Lz φ + Sr (r) + Sθ (θ) とすると、
(
(
)(
)
)
N2 2
b
dSθ
dSr
2 2 2
2
r N K δ1 = r K (E − ωLz ) −
Lz − r N K 1 −
−N
r
dr
dθ
sin2 θ
2
2
2
2
2
2
(
r K
N
(
2
b
1−
r
2
)(
2
dSr
dr
)
)2
− (E − ωLz ) + N δ1
2
2
((
+N
2
dSθ
dθ
)2
L2z
+
sin2 θ
(158)
)
=0
(159)
となる。Q を分離定数とすると、
(
)(
)2
(
)
b
dSr
Q
2
2
N 1−
= (E − ωLz ) − N δ1 + 2 2 ≡R
r
dr
r K
2
(
dSθ
dθ
)2
=Q−
26
L2z
≡Θ
sin2 θ
(160)
(161)
5.1
数値計算方法
第5章
ワームホール周辺の光子の軌道計算
と表される。したがって、
1
S = δ1 λ − Et + Lz φ +
2
∫
r
√
R
(
) dr +
2
N 1 − rb
∫
θ
√
Θdθ
(162)
と改めて書き直すことができる。
1
∂S
=
∂Q
2
∫
∫
1 θ 1 ∂Θ
∂R
√
√
dr +
dθ = 0
(
)
2
Θ ∂Q
N 2 1 − rb R ∂Q
∫ r
∫ θ
N
1
1
√(
√ dθ
∴
dr
=
)
2
2
Θ
1 − rb R r K
r
1
(163)
(164)
なので、
r2 K 2
dr
=
dθ
N
∫
∂S
1
1
= λ+
∂δ1
2
2
r
∫
√(
)
1 − rb R
Θ
(165)
1
∂R
1
√
(
) ∂δ dr + 2
b
1
2
N 1− r R
r
∴λ=
√(
N
1−
b
r
)
∫
θ
1 ∂Θ
√
dθ = 0
Θ ∂δ1
dr
(166)
(167)
R
なので、
√(
)
1 − rb R
dr
=
dλ
√ N
dθ
Θ
= 2 2
dλ
r K
(168)
(169)
と求まる。
∂S
1
= −t +
∂E
2
∫
∫
r
∫
1
1
∂R
√
(
) ∂E dr + 2
b
N2 1 − r R
r
∴t=
θ
1 ∂Θ
√
dθ = 0
Θ ∂E
(170)
1
E − ωLz
√(
dr
)
N
b
1− r R
(171)
dt
E − ωLz
=
dλ
N2
(172)
なので、
と求まる。
∂S
1
= −φ +
∂Lz
2
∫
∴φ=
r
∫
r
√(
∂R
√
dr +
(
)
N 2 1 − rb R ∂Lz
1
∫
ω(E − ωLz )
dr +
)
N
1 − rb R
1
27
∫
θ
1 ∂Θ
√
dθ = 0
Θ ∂Lz
θ
1 (−2Lz )
√
dθ
Θ sin2 θ
(173)
(174)
5.1
数値計算方法
第5章
ワームホール周辺の光子の軌道計算
なので、
(
N
2
dφ
dλ
)
= ω(E − ωLz ) +
N 2 Lz
r2 K 2 sin2 θ
(175)
と求まる。以上より、測地線方程式が求まったので、これを数値計算で解く。
5.1.2 初期条件
観測者の周辺では局所ローレンツ系であるということから初期条件を求める。測地線方程式より、４元運動
量は、
Pt =
Pr =
E − ωLz
2
√(N
)
1 − rb R
(176)
√ N
Θ
Pθ = 2 2
r K
ω(E − ωLz )
Lz
φ
P =
+ 2 2 2
2
N
r K sin θ
(177)
(178)
(179)
である。したがって、Pµ = gµν P ν より、
Pt = −E
√
1
R
(
)
Pr =
N
1 − rb
√
Pθ = Θ
Pφ = Lz
(180)
(181)
(182)
(183)
と書ける。四脚場ベクトルは、
eµ(t) =
eµ(r) =
eµ(θ) =
eµ(φ) =
[
[
[
[
1
N
ω
N
0
√
1−
0
0
0 0
]
0 0
b
r
1
rK
0
0
]
0 0
]
1
rK sin θ
]
(184)
である。したがって、局所ローレンツ系における４元運動量は、
P (t) = −P(t) = −eµ(t) Pµ =
√
P (r) = P(r) = eµ(r) Pµ =
√
Θ
(θ)
P = P(θ) =
rK
Lz
P (φ) = P(φ) =
rK sin θ
28
E − ωLz
N
1−
N
b
r
(185)
(186)
(187)
(188)
5.1
数値計算方法
第5章
ワームホール周辺の光子の軌道計算
と書ける。いま、null 測地線方程式に対して δ1 = 0 であり、ξ = Lz /E/, η = Q/E 2 という新たなパラメー
ターを導入し、アフィンパラメーターを λnew = Eλ とすると、
R = (1 − ωξ)2 − N 2
Θ=η−
η
r2 K 2
ξ
sin2 θ
(189)
(190)
と書き直すことができる。
さらに、無限遠にいる観測者が見る像の見かけの位置を与える α と β の天球座標に直す。θ0 をワームホー
ルの回転軸と観測者の視線の間の角度とすると、
(
rP (φ)
P(t)
)
ξ
sin
θ0
r*∞
√(
( (θ) )
)
rP
ξ2
β=
=
η−
P(t) r*∞
sin2 θ0
α=
=
(191)
(192)
となる。ここで r*∞ のとき、
P (t) = P(t) = 1
(193)
とすると、
αP (t)
α
= = rK sin θP φ
r
r
βP (t)
β
=
= = rKP θ
r
r
P (φ) =
(194)
P (θ)
(195)
より、
(
)
α
α
= 2
2 K sin θ r*∞
r
r
sin
θ
(
)
β
β
Pθ =
= 2
r2 K r*∞
r
Pφ =
(196)
(197)
であることが分かる。null 条件 P µ Pµ = 0 より、
√
P
(r)
=− 1−
α2 + β 2
r2
(198)
であるので、以上より、初期条件における測地線方程式が求まった。
5.1.3 不安定軌道
ワームホールスロートのまわりの光子の軌道には、ワームホールスロートに接触せず無限遠から無限遠にぬ
けていく軌道、ワームホールスロートに接触して我々の宇宙から別の宇宙にぬけていく軌道、そして、ワーム
ホールスロートの周りをぐるぐると円を描く不安定軌道と３種類存在することが分かっている。不安定軌道に
なる条件は既に確認済みである。この条件を解けば、不安定軌道になるときの衝突パラメーターを求めること
ができる。具体的に、赤道上面での不安定軌道になるとき、すなわち β = 0 のときの衝突パラメーターを求め
る。いま、式 (168) より、
(
dr
dλ
)2
](
)
[
(1 − ωξ)2 − N 2 r2ηK 2 1 − rb
=
N2
29
(199)
5.1
数値計算方法
第5章
ワームホール周辺の光子の軌道計算
である。θ = π/2 を代入すると、
(
dr
dλ
)2
[(
](
)2
)
2a
α2
1
=
1− 3α − 2
1−
= −V
r
r
r
(200)
である。これより、ポテンシャル V と r のグラフは、
図7
a = 0, θ = π/2, β = 0 のときのポテンシャル
ワームホールの回転の大きさは a = 0。このとき、不安定軌道になるのは α = 1.0, −1.0 のときであること
が分かる。
30
5.1
数値計算方法
第5章
図8
ワームホール周辺の光子の軌道計算
a = 0.5, θ = π/2, β = 0 のときのポテンシャル
ワームホールの回転の大きさは a = 0.5。このとき、不安定軌道になるのは α = 0.43, −2.59 のときである
ことが分かる。
c = 0 のときは測地線方程式を求める際に変数分離をすることができたので、理論的に不安定軌道を求める
ことができる。式 (168) より、測地線方程式は、
(
N2
dr
dλ
)2
[
= (1 − ωξ)2 − N 2
(
)
η ]
b
1
−
= −V
r2 K 2
r
(201)
2a
r3
(202)
(
)
b
1−
≡g
r
(203)
N = K = 1, b = 1, ω =
である。ここで、
(1 − ωξ)2 − N 2
η
≡f
r2 K 2
とおくと、不安定軌道を求めるには、
V = fg = 0
(204)
∂V
= f 0 g + f g0 = 0
∂r
(205)
を満たせば良い。
31
5.1
数値計算方法
第5章
ワームホール周辺の光子の軌道計算
f = 0 のとき、
(
)2
2a
η
V = 1− 3ξ − 2 =0
r
r
(206)
(
)(
)
∂V
2a
2a
η
=2 1− 3ξ
3 4ξ =2 3 =0
∂r
r
r
r
(207)
と書ける。これを解くと、
r3 ±3r3
8a
(208)
r3
2a
(209)
η=0
(210)
ξ=
と求まる。ただし、
ξ=
のとき、
であり、これは不適切である、したがって、
r3
4a
9r2
η=
4
ξ=−
(211)
(212)
である。これと、式 (191)(192) より、
r3
α=−
√4a sin θ
9r2
r6
β=
−
2
4
16a sin2 θ
(213)
(214)
と求まる。g = 0 のとき、
V =1−
b
=0
r
(215)
∂V
1
= 2 (1 − 4aξ + 4a2 ξ 2 − η) = 0
∂r
r
(216)
η = 1 − 4aξ + 4a2 ξ 2
(217)
(1 − 4a2 sin2 θ)α2 + 4a sin θα + β 2 − 1 = 0
(218)
と書ける。これを解くと、
これと、式 (191)(192) より、
と求まる。これらより、理論的に不安定軌道になる衝突パラメータ−を求めることができる。
また、c 6= 0 のときは測地線方程式を求める際に変数分離をすることができないので、理論的に不安定軌道
を求めることはできない。
32
5.2
数値計算の結果
第5章
ワームホール周辺の光子の軌道計算
5.2 数値計算の結果
今回の計算では、C = G = M = 1、初期条件 r = 200, φ = 0、として計算を行う。なお、ブラックホール
の数値計算の際に用いた計算コードのメトリック等を変えている。
図9
c = 0, a = 0, θ = π/2, β = 0 のときの軌道
緑：ワームホールスロート (半径 1)
c = 0, a = 0 の Morris と Thorne のワームホールである。ワームホールスロートと静止限界は一致し、エ
ルゴ領域は見られない。衝突パラメーターを変化させていくと、無限遠から無限遠に抜けていき、α = 1.0 の
ときに円軌道になる。それよりも内側になると、ワームホールスロートに接触したのちまた無限遠にはね返さ
れるような軌道になる。そして再び α = −1.0 のとき円軌道になり、また無限遠から無限遠に抜けていく。
この軌道を、３次元ユークリッド空間に埋め込みを行うと、
33
5.2
数値計算の結果
第5章
図 10
ワームホール周辺の光子の軌道計算
埋め込み
３次元ユークリッド空間に埋め込みを行うと、光子の軌道はスロートを上から下に通過していることが分か
る。ワームホールスロートを介して、一方は我々の宇宙、もう一方は我々の宇宙とは殆ど因果関係を持たない
別の空間に隔てられている。つまり、我々の宇宙にいる観測者から飛ばした光子が、ひとたびワームホールス
ロートに達すると、我々とは別の宇宙へと抜けていくのである。ここで、３次元ユークリッド空間への埋め込
みについて説明する。c = 0, a = 0 より計量は、
(
)−1
b
ds = −dt + 1 −
dr2 + r2 (dθ2 + sin2 θdφ2 )
r
2
2
(219)
であり、t, θ, φ 一定にすると、
)−1
(
b
dr2
ds = 1 −
r
2
= dz 2 + dr2
( )2
∂z
=
dr2 + dr2
∂r
[( )
]
2
∂z
=
+ 1 dr2
∂r
(220)
したがって、
(
∂z
∂r
)2
=
34
b
r−b
(221)
5.2
数値計算の結果
第5章
√
z = ±2 b(r − b)
ワームホール周辺の光子の軌道計算
(222)
となる。これより、(r, θ) 平面を、式 (222) を用いて仮想的に z 軸を設定して (r, θ, z) と３次元化する。
図 11 c = 0, a = 0.5, θ = π/2, β = 0 のときの軌道
緑：事象の地平面とエルゴ領域 (半径 1.866)
この計算では、光子がワームホールスロート達すると計算を止めている。ワームホールの回転は反時計回
り。衝突パラメーターを変化させていくと、無限遠から無限遠に抜けていき、α = 0.43 付近で円軌道になる。
それよりも内側になると、エルゴ領域内でワームホールの回転と同じ方向に引きずられたのち別の宇宙へと抜
けていく。そして、再び α = −2.59 付近で円軌道になり、また無限遠から無限遠に抜けていく。ここで、理
論値と数値計算の結果に多少のズレが生じているが、これは、実際には観測者は無限遠にいるとしたが、数値
計算の初期条件で r = 200 と設定したためだと考えられる。したがって、r のを大きくすればするほど計算の
誤差は少なくなっていくと思われる。
α と β の変化させていき、不安定軌道になる値をプロットしていくと、ワームホールのシルエットが描ける。
35
5.2
数値計算の結果
第5章
ワームホール周辺の光子の軌道計算
図 12 c = 0 のワームホールシャドー
a はワームホールの回転、th は回転軸と観測者の間の角度を表している。この数値計算の結果は、式
(213)(214)(218) より理論的に求めたものと正しく一致している。また、c を変化させていくと、
36
5.2
数値計算の結果
第5章
図 13
ワームホール周辺の光子の軌道計算
c を変化させたときのワームホールシャドー
a = 0, θ = π/2、赤：c = 0、緑：c = 8、青：c = 16
c を変化させていくと、細長いシルエットになることが分かった。
37
第6章
第6章
まとめ
まとめ
本卒業研究では、はじめに一般相対論の基礎を学んだ。ここでは、測地線方程式やアインシュタイン方程式
について導出した。
次に、カー解で測地線方程式を解き、数値計算を行い、不安定軌道になる点をプロットすることで、Rohta
Takahashi(2005) のブラックホールシャドーを再現した。ここでは、数値計算と理論がよく一致し、計算コー
ドが正しいことが分かった。
そして、Morris=Thorne のワームホールモデルと Edward Teo のワームホールモデルについて比較し、
Morris=Thorne のモデルではスロートが exotic matter で満たされているが、Teo のモデルでは、exotic
matter が存在しない領域があるということが分かった。また、それぞれの θ に対する c の制限、c に対する
θ の制限を得、具体的に求めることができた。したがって、本研究では Edward Teo のモデルを採用し、測地
線方程式を解き、数値計算を行い、不安定軌道になる点をプロットすることで、ワームホールシャドーを描い
た。c = 0 のときは、理論曲線が描けて数値計算を理論がよく一致した。c 6= 0 のときは、c を変化させるにつ
れて細長いシルエットになっていくということが分かった。
38
第6章
まとめ
謝辞
本研究および卒業論文作成を行うにあたり、スタッフの皆様、研究室の先輩方から多くのご協力をいただい
まして誠にありがとうございました。
杉山教授には、１年生の基礎セミナーから始まり、４年生前期のセミナーで英文訳から宇宙論の基礎を丁寧
に指導していただきました。私の大学４年間を通じて、杉山教授は学習面から、生活面等様々な場面でたくさ
んの助言をしていただいたこと、大変感謝しています。また、松原准教授、市來助教、日影特任助教、門田特
任助教、新田さん、関口さんには、中間発表の際に様々な助言と質問などをしていただきました。自分の知識
の未熟さを改めて実感し、理解を深めていくことができました。
特に、直接指導していただきました新田さんには本当に感謝しております。研究に対する稚拙な質問をした
り、同じ質問を繰り返してしまったりと多大なる迷惑をかけてしまいましたが、お忙しい中、私の質問等に優
しく、丁寧に指導していただきました。新田さんのおかげに無事卒業論文を書き上げることができました。本
当にありがとうございました。大学院進学後も迷惑をかけると思いますが、よろしくお願い致します。
また、日影さんにも直接指導していただきました。最初に、研究のテーマに対するアイデアをいただき、卒業
論文作成への道筋をたてていただきました。本当にありがとうございました。
研究室の先輩方には常日頃から優しくしていただき、この 1 年間はとても充実したものになりました。プロ
グラミングや研究に関して自分の質問に答えていただき、時には優しく、時に厳しく指導していただく中で成
長できたと実感しています。また、研究室以外の場でも優しく接していただきました。本当に感謝致します。
稲垣さん、竹内さん、渋沢さん、卒業後もがんばってください。
そして、同期の皆にはセミナーや大学院試験の際に一緒になって議論したり、質問し合ったりとてもお世話
になりました。特に、照屋さんには公私共にお世話になりました。大学院に進んでもがんばっていきましょ
う。
最後に、心配しつつもいつも暖かく見守り、応援してくれた家族に心から感謝致します。
39
参考文献
参考文献
参考文献
参考文献
[1] 杉山 直:『相対性理論』講談社
[2] 須藤 靖:『一般相対論入門』日本評論社
[3] Herbert Goldstein, John Safko, Charles Poole：古典力学（上）吉岡書店
[4] Herbert Goldstein, John Safko, Charles Poole：古典力学（下）吉岡書店
[5] Stuart L. Shapiro, Saul A. Teukolsky:『Black Holes, White Dwarfs and Neutron Stars』Wiley-VCH
[6] S. Chandrasekhar:『The Mathematical Theory of Black Holes』Oxford University Press, USA
[7] Charles W. Misner, Kip S. Thorne, John Archibald Wheeler:『GRAVITATON』W H Freeman & Co
[8] 高橋労太:”銀河中心のブラックホールの影” , 天文月報 2006 年 1 月
[9] Rohta Takahashi:”Black Hole Shadows of Charged Spinning Black Holes” (2005) [arXiv:astroph/0505316v1]
[10] Edward Teo:”otating traversable wormholes” (1998) [arXiv:gr-qc/980398v2]
[11] Petya G.Nedkova,Vassil Tinchev,Stoytcho S. Yazadjiev :”The Shadow of a Rotating Traversable
Wormhole” (2013) [arXiv:gr-qc/1307.7647v1]
40
参考文献
参考文献
付録：カーブラックホール時空上の光子軌道を求めるプログラム
#include <stdio.h>
#include<math.h>
#define M 1.0 #ブラックホールの質量
#define pi 3.14159265359
#define alpha {} #天球座標
#define beta {} #天球座標
#define a {} #角運動量
/\ast X,Y,Z definition \ast/
double X(double r,double th)
{
return pow(r,2.0)+pow(a*cos(th),2.0);
}
double Y(double r)
{
return pow(r,2.0)-2*M*r+pow(a,2.0);
}
double Z(double r,double th)
{
return pow(pow(r,2.0)+pow(a,2.0),2.0)-pow(a*sin(th),2.0)*Y(r);
}
/* metric g_{\mu\nu} */
double g00(double r,double th)
{
return 1-2*M*r/X(r,th);
}
double g03(double r,double th)
{
return 2*a*M*r*pow(sin(th),2.0)/X(r,th);
}
double g11(double r,double th)
{
return -X(r,th)/Y(r);
}
double g22(double r,double th)
{
return -X(r,th);
}
double g33(double r,double th)
{
return -(pow(r,2.0)+pow(a,2.0)+2*pow(a,2.0)*M*r*pow(sin(th),2.0)/X(r,th))*pow(sin(th),2.0);
}
/* metric g^{\mu\nu} */
41
参考文献
double
{
return
}
double
{
return
}
double
{
return
}
double
{
return
}
double
{
return
}
参考文献
f00(double r,double th)
Z(r,th)/(X(r,th)*Y(r));
f03(double r,double th)
2*a*M*r/(X(r,th)*Y(r));
f11(double r,double th)
1/g11(r,th);
f22(double r,double th)
1/g22(r,th);
f33(double r,double th)
-(Y(r)-pow(a*sin(th),2.0))/(X(r,th)*Y(r)*pow(sin(th),2.0));
/* derivative with respect to r, g_{\mu\nu,r}*/
double i00(double r,double th)
{
return 2*M*(pow(r,2.0)-pow(a*cos(th),2.0))/pow(X(r,th),2.0);
}
double i03(double r,double th)
{
return 2*a*M*pow(sin(th),2.0)*(pow(a*cos(th),2.0)-pow(r,2.0))/pow(X(r,th),2.0);
}
double i11(double r,double th)
{
return 2*(X(r,th)*(r-M)-r*Y(r))/pow(Y(r),2.0);
}
double i22(double r,double th)
{
return -2*r;
}
double i33(double r,double th)
{
return -2*pow(sin(th),2.0)*(r+pow(a,2.0)*M*pow(sin(th),2.0)*(pow(a*cos(th),2.0)-pow(r,2.0))/X(r,th));
}
/* derivative with respect to theta, g_{\mu\nu,th}*/
double j00(double r,double th)
{
return -4*pow(a,2.0)*M*r*sin(th)*cos(th)/pow(X(r,th),2.0);
}
double j03(double r,double th)
42
参考文献
{
return
}
double
{
return
}
double
{
return
}
double
{
return
}
参考文献
4*a*M*r*(sin(th)*cos(th)*X(r,th)+pow(a,2.0)*pow(sin(th),3.0)*cos(th))/pow(X(r,th),2.0);
j11(double r,double th)
2*pow(a,2.0)*sin(th)*cos(th)/Y(r);
j22(double r,double th)
2*pow(a,2.0)*cos(th)*sin(th);
j33(double r,double th)
-2*sin(th)*cos(th)*(pow(r,2.0)+pow(a,2.0)+2*pow(a,2.0)*M*r*pow(sin(th),2.0)/X(r,th))-a*pow(sin(th),2.0)*j
/* Christoffel symbols */
double G001(double r,double th)
{
return f00(r,th)*i00(r,th)*0.5+f03(r,th)*i03(r,th)*0.5;
}
double G002(double r,double th)
{
return f00(r,th)*j00(r,th)*0.5+f03(r,th)*j03(r,th)*0.5;
}
double G013(double r,double th)
{
return f00(r,th)*i03(r,th)*0.5+f03(r,th)*i33(r,th)*0.5;
}
double G023(double r,double th)
{
return f00(r,th)*j03(r,th)*0.5+f03(r,th)*j33(r,th)*0.5;
}
double G100(double r,double th)
{
return -f11(r,th)*i00(r,th)*0.5;
}
double G103(double r,double th)
{
return -f11(r,th)*i03(r,th)*0.5;
}
double G111(double r,double th)
{
return f11(r,th)*i11(r,th)*0.5;
}
double G112(double r,double th)
{
return f11(r,th)*j11(r,th)*0.5;
}
43
参考文献
double
{
return
}
double
{
return
}
double
{
return
}
double
{
return
}
double
{
return
}
double
{
return
}
double
{
return
}
double
{
return
}
double
{
return
}
double
{
return
}
double
{
return
}
double
{
return
}
参考文献
G122(double r,double th)
-f11(r,th)*i22(r,th)*0.5;
G133(double r,double th)
-f11(r,th)*i33(r,th)*0.5;
G200(double r,double th)
-f22(r,th)*j00(r,th)*0.5;
G203(double r,double th)
-f22(r,th)*j03(r,th)*0.5;
G211(double r,double th)
-f22(r,th)*j11(r,th)*0.5;
G212(double r,double th)
f22(r,th)*i22(r,th)*0.5;
G222(double r,double th)
f22(r,th)*j22(r,th)*0.5;
G233(double r,double th)
-f22(r,th)*j33(r,th)*0.5;
G301(double r,double th)
f03(r,th)*i00(r,th)*0.5+f33(r,th)*i03(r,th)*0.5;
G302(double r,double th)
f03(r,th)*j00(r,th)*0.5+f33(r,th)*j03(r,th)*0.5;
G313(double r,double th)
f03(r,th)*i03(r,th)*0.5+f33(r,th)*i33(r,th)*0.5;
G323(double r,double th)
f03(r,th)*j03(r,th)*0.5+f33(r,th)*j33(r,th)*0.5;
44
参考文献
参考文献
/* geodesic equations */
double v0(double r,double th, double u0, double u1, double u2, double u3)
{
return -2*G001(r,th)*u0*u1-2*G002(r,th)*u0*u2-2*G013(r,th)*u1*u3-G023(r,th)*u2*u3;
}
double v1(double r,double th, double u0, double u1, double u2, double u3)
{
return -G100(r,th)*u0*u0-2*G103(r,th)*u0*u3-G111(r,th)*u1*u1-2*G112(r,th)*u1*u2-G122(r,th)*u2*u2-G133(r,th)*u3*u
}
double v2(double r,double th, double u0, double u1, double u2, double u3)
{
return -G200(r,th)*u0*u0-2*G203(r,th)*u0*u3-G211(r,th)*u1*u1-2*G212(r,th)*u1*u2-G222(r,th)*u2*u2-G233(r,th)*u3*u
}
double v3(double r,double th, double u0, double u1, double u2, double u3)
{
return -2*G301(r,th)*u0*u1-2*G302(r,th)*u0*u2-2*G313(r,th)*u1*u3-2*G323(r,th)*u2*u3;
}
int main(void)
{
double tmax,dt,t;
double x0,x1,x2,x3;
double u0,u1,u2,u3,l[4][4],k[4][4];
double A,B,C;
/* initial conditions */
x0=0;
x1=100.0;
x2=pi/2;
x3=0;
u0=1;
u1=-sqrt((1-(pow(alpha,2.0)+pow(beta,2.0))/pow(x1,2.0)));
u2=beta/pow(x1,2);
u3=alpha/(pow(x1,2)*sin(x2));
tmax=10000;
dt=0.005;
FILE *output;
output=fopen("test.txt","w");
/* 4th Runge Kutta Calculation*/
for(t=0;t<tmax;t+=dt) {
l[0][0]=v0(x1,x2,u0,u1,u2,u3)*dt;
k[0][0]=u0*dt;
l[1][0]=v1(x1,x2,u0,u1,u2,u3)*dt;
k[1][0]=u1*dt;
l[2][0]=v2(x1,x2,u0,u1,u2,u3)*dt;
45
参考文献
参考文献
k[2][0]=u2*dt;
l[3][0]=v3(x1,x2,u0,u1,u2,u3)*dt;
k[3][0]=u3*dt;
l[0][1]=v0(x1+k[1][0]/2,x2+k[2][0]/2,u0+l[0][0]/2,u1+l[1][0]/2,u2+l[2][0]/2,u3+l[3][0]/2)*dt;
k[0][1]=(u0+l[0][0]/2)*dt;
l[1][1]=v1(x1+k[1][0]/2,x2+k[2][0]/2,u0+l[0][0]/2,u1+l[1][0]/2,u2+l[2][0]/2,u3+l[3][0]/2)*dt;
k[1][1]=(u1+l[1][0]/2)*dt;
l[2][1]=v2(x1+k[1][0]/2,x2+k[2][0]/2,u0+l[0][0]/2,u1+l[1][0]/2,u2+l[2][0]/2,u3+l[3][0]/2)*dt;
k[2][1]=(u2+l[2][0]/2)*dt;
l[3][1]=v3(x1+k[1][0]/2,x2+k[2][0]/2,u0+l[0][0]/2,u1+l[1][0]/2,u2+l[2][0]/2,u3+l[3][0]/2)*dt;
k[3][1]=(u3+l[3][0]/2)*dt;
l[0][2]=v0(x1+k[1][1]/2,x2+k[2][1]/2,u0+l[0][1]/2,u1+l[1][1]/2,u2+l[2][1]/2,u3+l[3][1]/2)*dt;
k[0][2]=(u0+l[0][1]/2)*dt;
l[1][2]=v1(x1+k[1][1]/2,x2+k[2][1]/2,u0+l[0][1]/2,u1+l[1][1]/2,u2+l[2][1]/2,u3+l[3][1]/2)*dt;
k[1][2]=(u1+l[1][1]/2)*dt;
l[2][2]=v2(x1+k[1][1]/2,x2+k[2][1]/2,u0+l[0][1]/2,u1+l[1][1]/2,u2+l[2][1]/2,u3+l[3][1]/2)*dt;
k[2][2]=(u2+l[2][1]/2)*dt;
l[3][2]=v3(x1+k[1][1]/2,x2+k[2][1]/2,u0+l[0][1]/2,u1+l[1][1]/2,u2+l[2][1]/2,u3+l[3][1]/2)*dt;
k[3][2]=(u3+l[3][1]/2)*dt;
l[0][3]=v0(x1+k[1][2],x2+k[2][2],u0+l[0][2],u1+l[1][2],u2+l[2][2],u3+l[3][2])*dt;
k[0][3]=(u0+l[0][2])*dt;
l[1][3]=v1(x1+k[1][2],x2+k[2][2],u0+l[0][2],u1+l[1][2],u2+l[2][2],u3+l[3][2])*dt;
k[1][3]=(u1+l[1][2])*dt;
l[2][3]=v2(x1+k[1][2],x2+k[2][2],u0+l[0][2],u1+l[1][2],u2+l[2][2],u3+l[3][2])*dt;
k[2][3]=(u2+l[2][2])*dt;
l[3][3]=v3(x1+k[1][2],x2+k[2][2],u0+l[0][2],u1+l[1][2],u2+l[2][2],u3+l[3][2])*dt;
k[3][3]=(u3+l[3][2])*dt;
x0=x0+(k[0][0]+2*k[0][1]+2*k[0][2]+k[0][3])/6;
x1=x1+(k[1][0]+2*k[1][1]+2*k[1][2]+k[1][3])/6;
x2=x2+(k[2][0]+2*k[2][1]+2*k[2][2]+k[2][3])/6;
x3=x3+(k[3][0]+2*k[3][1]+2*k[3][2]+k[3][3])/6;
u0=u0+(l[0][0]+2*l[0][1]+2*l[0][2]+l[0][3])/6;
u1=u1+(l[1][0]+2*l[1][1]+2*l[1][2]+l[1][3])/6;
u2=u2+(l[2][0]+2*l[2][1]+2*l[2][2]+l[2][3])/6;
u3=u3+(l[3][0]+2*l[3][1]+2*l[3][2]+l[3][3])/6;
if(x1<1.0)return 0;
if(x1>100)return 0;
#光子の座標
A=x1*sin(x2)*cos(x3);
46
参考文献
参考文献
B=x1*sin(x2)*sin(x3);
C=x1*cos(x2);
fprintf(output,"%f %f %f\n",A,B,C);
/*printf("%f %f\n",x,y);*/
}
fclose(output);
return 0;
}
47