[CAS-Lab] ブラックホールをあなたの手の上に　一般相対性理論４

物理学ノート
前編

本稿は下記物理学ノート/一般相対性理論４で数式検証、グラフ描画に使用したMathematicaのプログラムの公開（個人での利用を許可）と簡単な解説となる。尚、一般相対性理論３の数式処理対象はあまり無かったので割愛した。

一般相対性理論4 [ブラックホール編]
本稿は一般相対性理論3 [重力方程式編]から続いての解説となる。重力方程式を解き、シュバルツシルト解を得てブラックホールがどのように記述されるかを解説る。引き続き数式処理 (CAS) を使ってブラックホールの構造等を視覚化を行い、数式の示すものを実感できるように心がける。
academians.net

リーマン幾何学で記述された、曲がった空間でのベクトルの動きや接続、曲率などの数式は難解で、その意味するところを直観的に理解することが難しいと思う。本稿に沿って、自分で処理を確認しながら実行することにより、その数式の意味を実感できるようになると思う。

サンプルプログラム実行前の準備

1) Mathematicaの環境構築

Mathematicaの環境が無い場合はこちらを参考に環境を整備してほしい。

尚、4.1節はMaxima版も準備したのでこちらも利用可能だ。

2) サンプルプログラム使用パッケージ（描画関数や座標変換関数等）のセーブ

以下ボタンで［.nb］ファイルをダウンロードしMathematicaで開いて実行できる。実行後いくつかのパッケージ（xxx.pac) がWindowsの場合デフォルトで”ドキュメント”　ディレクトリにセーブされるはずだ。

(*Mathematica*)

(*---------------基本パッケージ-----------------*)
　　：
　関数、変数定義　省略
　　：
(*---------------変換パッケージ-----------------*)
　　：
　関数、変数定義　省略
　　：
(*---------------描画パッケージ-----------------*)
　　：
　関数、変数定義　省略
　　：
(*---------ベクトル平行移動関連の共通関数--------*)
　　：
　関数、変数定義　省略
　　：
Save[(*cnd*)"udVectTranslate.pac",
 {
 セーブする関数、変数名のリスト  
 }];

上記ではソースコードの転記と解説を省略しました。ダウンロードした実物で内容は確認してください。

サンプルプログラムの実行

以下のボタンでサンプルプログラムをダウンロードしMathematicaで開いて実行できます。

最初にグローバル変数等の初期化と上記でセーブしたパッケージ の内、”udPlotPac1.pac” 、 “udElemntPac1.pac” と”udTransform1.pac”を読み込ます。これにより以降のプログラムで使用可能となります。

(*初期化、ユーザー定義関数パッケージの読み込み*)
(*Mathematica*)
ClearAll["Global`*"]

Get["udPlotPac1.pac"];
Get["udElemntPac1.pac"];
Get["udTransform1.pac"];
Get["udVectTranslate.pac"];

prON = False;

サンプルプログラムの実行

以上で作業準備は終わった。以降『一般相対性理論４』の数式検証とグラフ化に関するプログラムになる。各章節事にリンクを貼ってあるので内容を読み比べながらプログラムを実行してほしい。尚、内容はほぼ同じだが数式処理ソフトの差異で数式や画像の構成に多少の違いがあることをご承知ねがいたい。

4.1 シュバルツシルト解

ClearAll[c, t, r, theta, phi, swzgxx]
swzX = {x0, r, theta, phi};
swzgxx = dm[{g00[r], g11[r], -r^2, -r^2*Sin[theta]^2}];
% // mf

swzC3s = C3s[swzgxx, swzX];
PrintNonZero[swzC3s, "C3s_"];

swzR4 = R4[swzC3s, swzX];
PrintNonZero[swzR4, "R4_"];
swzGxx = Gxx[swzR4, swzgxx];
swzGxxUL = iv[swzgxx].swzGxx // Factor;
PrintNonZero[swzGxxUL, "Gul_"]
DSolve[swzGxxUL[[1, 1]] == 0, g11[r], r][[1, 1]] 
swzg11F = DSolve[swzGxxUL[[1, 1]] == 0, g11[r], r][[1, 1]] /. E^C[1] -> a // ExpandAll // Simplify
g00eq = DSolve[swzGxxUL[[1, 1]] == swzGxxUL[[2, 2]], g00[r], r][[1,1]]
swzg00F = (g00eq /. swzg11F // ExpandAll // Simplify) /. C[1] -> b
swzgxxF = swzgxx /. swzg00F /. swzg11F // MatrixForm

ニュートンの重力理論近似

ntXc = {x0, x1, x2, x3}
ntXf = {X0, X1, X2, X3}

rulexX = {X0 -> x0, X1 -> x1 - (h - (1/2)*g*x0^2/c^2), X2 -> x2,  X3 -> x3}

ntXcf = ntXf /. rulexX
Jh = Table[D[ntXcf[[ix]], ntXc[[iy]]], {ix, 4}, {iy, 4}];
eta = dm[{1, -1, -1, -1}];
ntgxx = tp[Jh].eta.Jh // MatrixForm
ntPhi = g*(x1 - h)
Phieq = Phi == ntPhi /. x1 -> h - (1/2)*g*x0^2/c^2 
x0x0eq = Solve[Phieq /. x0^2 -> x0x0, x0x0][[1, 1]] /. x0x0 -> x0^2
ntg00F = ntgxx[[1, 1, 1]] /. x0x0eq /. Phi -> -G*M/r

ntg00F == swzgxxF[[1, 1, 1]] // ExpandAll
aeq = Solve[(ntg00F == swzgxxF[[1, 1, 1]] /. b -> -1 // ExpandAll ),a](*[[1]]*)// Simplify
swzXFr = swzX /. x0 -> c*t

swzgxxFr = swzgxxF /. b -> -1 /. aeq[[1]] // ExpandAll

シュバルツシルト解の形状

const = {M -> 2.0 10^33, G -> 6.7 10^-14, c -> 3.0 10^8};
redshift = Solve[swzgxxFr[[1, 1, 1]] == 0 /. const, r][[1, 1]]
horizon = Solve[1/swzgxxFr[[1, 2, 2]] == 0 /. const, r][[1, 1]]
grRH = SphericalPlot3D[
   Evaluate[{r /. redshift, 1.01 r /. horizon}], {\[Theta], 0,Pi}, {\[Phi], 0, Pi}];
grCntP = Graphics3D[{Black, Sphere[{0, 0, 0}, 70]}];
Show[grRH, grCntP, ViewPoint -> {10, -5, 2}]

カー解の形状

Clear[m,r,a,theta,phi];

X={t,r,theta,phi};
gxxk = {{(a^2 - 2 m r + r^2 - a^2 Sin[theta]^2)/(r^2 + 
       a^2 Cos[theta]^2), 0, 
    0, (2 a m r Sin[theta]^2)/(r^2 + 
       a^2 Cos[theta]^2)}, {0, -((r^2 + a^2 Cos[theta]^2)/(a^2 - 
         2 m r + r^2)), 0, 0}, {0, 0, -r^2 - a^2 Cos[theta]^2, 
    0}, {(2 a m r Sin[theta]^2)/(r^2 + a^2 Cos[theta]^2), 0, 
    0, (-(a^2 + r^2)^2 Sin[theta]^2 + 
       a^2 (a^2 + r (-2 m + r)) Sin[theta]^4)/(r^2 + 
       a^2 Cos[theta]^2)}};
% // mf
redshift = Solve[gxxk[[1, 1]] == 0, r] // TrigReduce // ExpandAll
horizon = Solve[1/gxxk[[2, 2]] == 0, r] // Simplify
grRS = SphericalPlot3D[Evaluate[r /. redshift /. m -> 1 /. a -> .995], {theta, 0, Pi}, {phi, 0, Pi}, AspectRatio -> Automatic(*,Boxed\[Rule]False,Axes\[Rule]False*)];
grHZ = SphericalPlot3D[Evaluate[r /. horizon /. m -> 1 /. a -> .995], {theta, 0, Pi}, {phi, 0, Pi}, AspectRatio -> Automatic(*,Boxed\[Rule]False,Axes\[Rule]False*)];

Show[grRS, grHZ, PlotRange -> All, ViewPoint -> {10, -5, 2}]

4.2 質点運動の解析

前節で得られた計量を使いアインシュタインの確立した重力理論で質点がブラックホール近傍でどのように移動するかを実験する。本当に時空の曲率が4元移動している物体に作用し、あたかも重力が働いているかの様に見える現象を再現してみよう。質点の移動は測地線を進むベクトルを利用する。理論的な裏付けが取れた段階で修正される可能性が大きいことをご承知おき願いたい。

実験準備

(*実験専用の描画関数の設定*)
pFig420a[U123_, Chs123_, figC_, iPu_, iVu_, smpl_] := Module[{grgeo},
  (*測地線の軌跡を計算しPrに設定*)
  {Pr, Vr, Ar} = 
   Part[VectTrans[U123, Chs123, iPu, iVu, {0, 0, 0}, True, True, 
     n], {1, 2, 3}];
  (*Do[Print["P= ",Pr[[ix]],"  V= ",Vr[[ix]],"  A= ",Ar[[ix]]],{ix,1,
  3}];*)
  
  (*円柱座標からXYZ座標への変換*)
  pulist3d = 
   Table[{Pr[[ix, 2]] Cos[Pr[[ix, 3]]], Pr[[ix, 2]] Sin[Pr[[ix, 3]]], 
     Pr[[ix, 1]]}, {ix, 1, Length[Pr], smpl}];
  (*シュバルツシルト半径の描画*)
  grgeo = ListPointPlot3D[pulist3d, PlotStyle -> {figC}]
  ]
setIntV4[Be_, iPu_, rot_] := 
 Module[{eV4, pri}, 
  eV4 = {U2X[{1, 0}, Be][[1]], 
     dm[{1, 1/iPu[[2]]}].rt[rot].{U2X[{1, 0}, Be][[2]], 0}} // Flatten;
  Print["Inital iPu:{t,r,\[Phi]}=", iPu, ", 4V:{t,r,\[Phi]}=", eV4, 
   ", {t,\!\(\*SqrtBox[\(\*SuperscriptBox[\(x\), \(2\)] + \
\*SuperscriptBox[\(y\), \(2\)]\)]\)}=", {eV4[[1]], 
    Sqrt[eV4[[2]]^2 + (iPu[[2]]*eV4[[3]])^2]}];
  eV4
  ]
trp2xyz[trpP_] := {trpP[[2]]*Cos[trpP[[3]]], trpP[[2]]*Sin[trpP[[3]]],
   trpP[[1]]}
pIntV4[iPu_, eV4_, ds_, asc_: 1] :=
 Graphics3D[{Blue, Thick, Arrowheads[0.02], 
   Arrow[{trp2xyz[iPu], 
     trp2xyz[iPu] + asc (trp2xyz[iPu + ds eV4] - trp2xyz[iPu])/ds}]}], PlotStyle -> {figC}]
]
(*＝＝＝ブラックホール近傍での質点移動軌跡（測地線）挙動観測実験＝＝＝*)
U123 = {t, r, phi}

(*シュバルツシルト計量の設定*)
swg123a = 
  Table[swzgxxFr[[ix, iy]], {ix, {1, 2, 4}}, {iy, {1, 2, 4}}] /. 
      theta -> Pi/2 /. M -> c^2 /. G -> 1 /. c -> 1;
% // mf
(*クリストッフェル三指標の設定*)
Chs123a = C3s[swg123a, U123] // Simplify
(*ブラックホールのシュバルツシルト半径の描画情報を生成*)
horizona = Solve[1/swg123a[[2, 2]] == 0][[1]]
grHZa = SphericalPlot3D[
   Evaluate[r /. horizona], {theta, 0, 2 Pi}, {phi, 0, Pi}, 
   AspectRatio -> Automatic];
grHZc = Graphics3D[
   Cylinder[{{0, 0, 0}, {0, 0, 1000}}, Evaluate[r /. horizona]]];

prON = False;

シュバルツシルト計量を確認。シュバルツシルト半径は２だ。

\[\small \left( \begin{array}{ccc} 1-\frac{2}{r} & 0 & 0 \\ 0 & \frac{r}{2-r} & 0 \\ 0 & 0 & -r^2 \\ \end{array} \right) \]

実験１　ブラックホール脇を横切り移動

初期設定と設定した4元速度ベクトルの確認

n = 1000; ds = 1/n;
(* 初期位置　t=0,r=10,phi=0 *)
iPu = {0, 10, 0};
(* 初期4元速度ベクトルの設定
　Beta：0.6　
　空間ベクトル：ｘｙ平面でブラックホールに向かって約0.26Pi右方向） *)
eV4 = setIntV4[0.6, iPu, 0.74 Pi];
grIntV4 = pIntV4[iPu, eV4, ds, 5];
grStrt = ListPointPlot3D[{trp2xyz[iPu]}, PlotStyle -> Blue];
Show[grStrt, grIntV4, grHZc, AxesOrigin -> {0, 0, 0}, 
 PlotRange -> 10 {{-1, 1}, {-1, 1}, 1 {-0.1, 1}}, 
 ViewPoint -> {0, 0, 1000}]

シミュレーションの実行

iVu = 50 eV4 ds // N;
(*質点の軌跡（測地線）の計算と描画情報セーブ*)
grgeo2a = pFig420a[U123, Chs123a, Red, iPu, iVu, 10];
(*描画*)
Show[grgeo2a, grHZc, grStrt, grIntV4, AxesOrigin -> {0, 0, 0}, 
 PlotRange -> 10 {{-1, 1}, {-1, 1}, 5 {-0.1, 1}}, 
 ViewPoint -> {0, 0, 1000}]

空間X-Y面にて、右から侵入しブラックホールに引き寄せられて進行方向が変わってしまった例となる。この例でも重力はうまく再現できている様に思われる。

実験２　ブラックホールを一周して離脱

eV4 = setIntV4[0.6, iPu, 0.788 Pi];
iVu = 50 eV4 ds // N
(*質点の軌跡（測地線）の計算と描画情報セーブ*)
grgeo2a = pFig420a[U123, Chs123a, Red, iPu, iVu, 10];
(*描画*)
Show[grgeo2a, grHZc, AxesOrigin -> {0, 0, 0}, 
 PlotRange -> 10 {{-1, 1}, {-1, 1}, 10 {-0.1, 1}}, 
 ViewPoint -> {0, 0, 1000}]

出力結果

空間X-Y面にて、右から侵入しブラックホールに近接たまま一周し左上方向に離脱した。ブラックホールを一周する間、プロット点の間隔が伸びており、侵入時より高速で運動ている事がわかる。

実験３　ブラックホールを円軌道で周回する

(*実験円軌道*)
U123 = {t, r, phi}
iPu = {0, 5, 0 Pi}(*t,r,phi*)

(*シュバルツシルト計量の設定*)
swg123a = 
  Table[swzgxxFr[[ix, iy]], {ix, {1, 2, 4}}, {iy, {1, 2, 4}}] /. 
      theta -> Pi/2 /. M -> c^2 /. G -> 1 /. c -> 1;
% // mf
(*クリストッフェル三指標の設定*)
Chs123a = C3s[swg123a, U123] // Simplify
(*円軌道の測地線方程式 \[CapitalGamma].V.V=0 で r'=0 \
の条件で角速度：theta'を求めて初期速度（光速度比：Be)を算出*)
res = Solve[Chs123a.{t', 0, phi'}.{t', 0, phi'} == 0 // Simplify, 
   phi'];
Print["res=", res]
Be = (iPu[[2]] phi' /. res /. t' -> 1 /. r -> iPu[[2]] // N)[[2]];
(*結果は Be=1/Sqrt[r] *)
Print["Be=", Be]
(*1周期座標時間を計算し通過目標の青い時空点を作成*)
2 Pi iPu[[2]]^(3/2) // N
cyc1 = trp2xyz[{2 Pi iPu[[2]]^(3/2), iPu[[2]], iPu[[3]]}] // Simplify
grCyc1 = ListPointPlot3D[{Evaluate[cyc1 // N]}, PlotStyle -> Green];

n = 5000; ds = 1/n;
(* 初期4元速度ベクトルの設定 *)
eV4 = setIntV4[Be, iPu, 0.5 Pi];
Print["eV4=", eV4, ": ", 1 (trp2xyz[iPu + ds eV4] - trp2xyz[iPu])/ds]
grIntV4 = pIntV4[iPu, eV4, ds, 4];
grStrt = ListPointPlot3D[{trp2xyz[iPu]}, PlotStyle -> Blue];
(*Graphics3D[{Blue,Thick,Arrowheads[0.02],Arrow[{trp2xyz[iPu],trp2xyz[\
iPu]+5(trp2xyz[iPu+ds eV4]-trp2xyz[iPu])/ds}]}];*)
Show[grStrt(*grCyc*), grIntV4, grHZc, AxesOrigin -> {0, 0, 0}, 
 PlotRange -> 10 {{-1, 1}, {-1, 1}, 1 {-0.1, 1}}, 
 ViewPoint -> {-100, -70, 200}]
iVu = 200 eV4 ds // N;
(*質点の軌跡（測地線）の計算と描画情報セーブ*)
grgeo3a = pFig420a[U123, Chs123a, Red, iPu, iVu, 20];

Show[grStrt, grgeo3a, grHZc, grCyc1, AxesOrigin -> {0, 0, 0}, 
 PlotRange -> 10 {{-1, 1}, {-1, 1}, 30 {-0.1, 1}}, 
 ViewPoint -> {-100, -70, 200}]

出力結果

設定値に非常に敏感で、少しのズレですぐに、外に広がるかブラックホールに落下する。

実験４　静止質点のブラックホールへの落下

n = 5000; ds = 1/n;
(* 初期位置　t=0,r=5,phi=0 *)
iPu = {0, 10, 0};
(* 初期4元速度ベクトルの設定
　Beta：0　（静止）
　空間ベクトル：ｘｙ平面でブラックホール方向 *)
eV4 = setIntV4[0, iPu, 0 Pi]
grIntV4 = pIntV4[iPu, eV4, ds, 4];
grStrt = ListPointPlot3D[{{iPu[[2]], 0, 0}}, PlotStyle -> Blue];
Show[grStrt, grHZc, grIntV4, AxesOrigin -> {0, 0, 0}, 
 PlotRange -> 11 {{-1, 1}, {-1, 1}, 1 {-0.1, 1}}, 
 ViewPoint -> {100, 0, 0}]
iVu = 50 eV4 ds // N
(*質点の軌跡（測地線）の計算と描画情報セーブ*)
grgeo4a = pFig420a[U123, Chs123a, Red, iPu, iVu, 1];
pulistRZ = Table[Part[pulist3d[[ix]], {1, 3}], {ix, Length[pulist3d]}];
(*描画*)
Show[grgeo4a, grHZc, AxesOrigin -> {0, 0, 0}, 
 PlotRange -> 11 {{-1, 1}, {-1, 1}, 7 {-0.1, 1}}, 
 ViewPoint -> {100, 0, 0}]

出力結果

r=10 静止。速度ベクトルは真上に1。

最初は r=10 静止していた(つまり4元速度の時間成分１、空間成分０）質点が、重力により徐々にブラックホールに引き寄せられていく。

x-z面の軌跡情報を理論値との比較のため生成、保存

sz = Length[pulistRZ]
pulistRZc = {}; pulistRZs = {}; ix = 1;
While[(ix <= sz) && (pulistRZ[[ix, 2]] < 70), 
 pulistRZc = Append[pulistRZc, pulistRZ[[ix]]]; 
 pulistRZs = 
  Append[pulistRZs, {pulistRZ[[ix, 
     1]], (1 - (pulistRZ[[ix + 1, 1]] - 
          pulistRZ[[ix, 1]])/(pulistRZ[[ix + 1, 2]] - 
           pulistRZ[[ix, 2]])^0.5 ix)}]; ix++]
gr1 = Show[ListPlot[{pulistRZc(*,pulistRZs*)}, PlotStyle -> {Red}], 
  PlotRange -> {{0, 12}, {0, 60}}, AxesOrigin -> {2, 0}]

座標時間の数値解析値

理論式をeqCTimeに設定し、グラフを生成（青の点線）。上記シミュレーション結果（赤）と重ねて描画。

ClearAll[m, r, r0, t]
eqCTime = (Sqrt[-2 m + 
     r0] (4 Sqrt[2] m^(3/2)
        ArcTan[(Sqrt[2] Sqrt[m] Sqrt[r0 - r[t]])/(
        Sqrt[2 m - r0] Sqrt[r[t]])] Sqrt[r0 - r[t]] - 
      Sqrt[2 m - r0] r0 Sqrt[r[t]] + Sqrt[2 m - r0] r[t]^(3/2) + 
      Sqrt[2 m - r0] Sqrt[
       r0] (4 m + r0) ArcSin[Sqrt[r[t]]/Sqrt[r0]] Sqrt[
       1 - r[t]/r0]))/(Sqrt[4 m - 2 r0] Sqrt[m (r0 - r[t])])
Ttime = (eqCTime /. r0 -> 10 /. r[t] -> r) // Simplify // PowerExpand
c1 = Abs[Ttime /. r -> 9.999999 /. m -> 1]
plT = Plot[
   Evaluate[Ttime /. {m -> 1, r0 -> 10}] + c1, {r, 2.001, 9.999}, 
   PlotStyle -> {Dashed, Blue}, PlotRange -> {{0, 12}, {0, 60}}, 
   AxesOrigin -> {2, 0}];
Show[plT]
Show[gr1, plT]

座標時間の理論値

座標時間の理論値vs数値解析値

見た目上はピッタリ重なっており、シミュレーションは正しそうだ。

おわりに

以上で一旦終わりとする。これから先は宇宙論で本格的なスタディが必要そうだ。もし、時間が出来て軽く出来て面白い実験が見つかったら追記しよう。

一般相対性理論4 [ブラックホール編]
本稿は一般相対性理論3 [重力方程式編]から続いての解説となる。重力方程式を解き、シュバルツシルト解を得てブラックホールがどのように記述されるかを解説る。引き続き数式処理 (CAS) を使ってブラックホールの構造等を視覚化を行い、数式の示すものを実感できるように心がける。
academians.net
Post Views: 2

コメント