[size=150][b]このワークシートは[url=https://www.geogebra.org/m/twxxx3yq]Math by Code[/url]の一部です。[br][/b][/size][size=150][size=100][br]前回までにアインシュタイン方程式に到達しました。[br]この方程式は、[b]非線形[/b]で[b]座標系[/b]が任意で[b]解釈[/b]の自由度もあるという、[br][color=#0000ff][b]ゆるく複雑さの高いもの[/b][/color]です。[br]条件をつけて解釈することで、宇宙にまつわる様々な事象を予測する道具にすることができます。[br][br]今回は、ブラックホールの存在予測につながったシュワルツシルト解について学ぼう。[/size][/size]

[br][b][color=#0000ff]条件０・アインシュタインの方程式はRij −1/2gij R＝Tij。[br]条件１・天体の外部が真空、つまり、Tij=0。[br]条件２・物質分布は球対称、つまり、球座標系。物質は原点に集中している。[br][/color][/b]とします。[br][br]条件０，１から曲率はゼロとなり、[b][color=#0000ff]R=Rij＝0[/color][/b][br]条件２から、球面座標系の[b]一般座標をxi=(c, r, t, p)とする。cはctの略、rは半径、tはθ、pはφ[/b]。[br]共変計量テンソル [b]g_ij=diag(A(r), -B(r), -r[sup]2[/sup], -r[sup]2[/sup] sin[sup]2[/sup]t )[/b][br]としましょう。

[size=150][b]＜リッチテンソルを求める＞[br][/b][b][br][/b][size=100]リッチテンソルを求めるには、[br]リーマンテンソルのRi_jklの２つの軸（1番，3番）を縮約します。[br][/size][b]つまり、[/b][color=#0000ff]R[sub]ij[/sub]=Σ[sub]m[/sub] R[sup]m[/sup]imj[/color][b]を計算しましょう。[br][br][/b][/size][color=#9900ff][b][u][size=150]課題：リッチテンソルを求めるにはどうしたらよいでしょうか。[br][/size][/u][/b][/color][br]前もって、リーマンテンソルR[sup]i[/sup]jklの辞書を作っておきましょう。[br]さらにその前に、[br]計算の高速化のためにg_mnから[b]g^mn=g_mn.inv()[/b]を求めておきます。[br]そして、クリ１とクリ２を同じループで作りながら、これもgamma_ijkとして辞書化しておきます。[br][br]そうすれば、リーマンテンソルR[sup]i[/sup]jklの辞書作りは、[br][br] [b]diff(gamma[i][j][l], xs[k])-diff(gamma[i][j][k], xs[l])[/b]に[br]ループをまわして、[br]term3 = 0[br] for m in range(dim):[br][b] term3 += gamma[i][k][m] * gamma[m][j][l] - gamma[i][l][m] * gamma[m][j][k][br][/b]この合計をさらに加えればできます。[br][br]#=====================================[br]import sympy as sp[br]from sympy import symbols, sin, diff, Matrix, Function, simplify[br][br]# 変数と関数の定義[br][b][color=#0000ff]c, r, t, p = symbols('c r t p')[br]xs = [c, r, t, p][br]dim = len(xs)[br][br]A = Function('A')(r)[br]B = Function('B')(r)[br][br]# 計量テンソル g_mn[br]M = sp.diag(A, -B, -r**2, -r**2 * sin(t)**2)[br]M_inv = M.inv()[br]# 辞書を使って gamma[i][j][k] = Γ^i_jk の形で保存[br][/color][/b]gamma = {}[br][br]print("計算中: クリストッフェル記号...")[br]for i in range(dim):[br] gamma[i] = {}[br] for j in range(dim):[br] gamma[i][j] = {}[br] for k in range(dim):[br] # 第1種クリストッフェル記号の計算[br] res2 = 0[br] for m in range(dim):[br][color=#0000ff] # 第1種 Γ_{m,jk} = 1/2 * (g_{mj,k} + g_{mk,j} - g_{jk,m})[br][/color] gamma1_mjk = 0.5 * (diff(M[m, j], xs[k]) + diff(M[m, k], xs[j]) - diff(M[j, k], xs[m]))[br][color=#0000ff] # 第2種 Γ^i_jk = g^im * Γ_{m,jk}[br][/color] res2 += M_inv[i, m] * gamma1_mjk[br] gamma[i][j][k] = simplify(res2)[br][br]# リーマンテンソルの辞書を初期化・格納[br]print("計算中: リーマンテンソル R^i_jkl ...")[br]riemann = {}[br]for i in range(dim):[br] riemann[i] = {}[br] for j in range(dim):[br] riemann[i][j] = {}[br] for k in range(dim):[br] riemann[i][j][k] = {}[br] for l in range(dim):[br][color=#0000ff] # R^i_{jkl} の計算[br] term1 = diff(gamma[i][j][l], xs[k])[br] term2 = diff(gamma[i][j][k], xs[l])[br] term3 = 0[br] for m in range(dim):[br] term3 += gamma[i][k][m] * gamma[m][j][l] - gamma[i][l][m] * gamma[m][j][k][br] val = simplify(term1 - term2 + term3)[br] riemann[i][j][k][l] = val [br][/color] # if val != 0:[br] # print(f"R^{xs[i]}_{xs[j]}{xs[k]}{xs[l]} = {val}")[br][br]# リッチテンソルの計算[br]print("\n計算中: リッチテンソル R_ij ...")[br]ricci = {}[br]for i in range(dim):[br] ricci[i] = {}[br] for j in range(dim):[br][color=#0000ff] # R_ij = Σ R^m_imj[br][/color] res_ricci = 0[br] for m in range(dim):[br][color=#0000ff] res_ricci += riemann[m][i][m][j][br][/color] [br] final_val = simplify(res_ricci)[br] ricci[i][j] = final_val[br] [br] if final_val != 0:[br] print(f"R_{xs[i]}{xs[j]} = {final_val}")[br]#=========================================[br][b][OUT][br]計算中: クリストッフェル記号...[br]計算中: リーマンテンソル R^i_jkl ...[br][br]計算中: リッチテンソル R_ij ...[br]R_cc = 0.5*Derivative(A(r), (r, 2))/B(r) - 0.25*Derivative(A(r), r)*Derivative(B(r), r)/B(r)**2 - 0.25*Derivative(A(r), r)**2/(A(r)*B(r)) + 1.0*Derivative(A(r), r)/(r*B(r))[br]R_rr = -0.5*Derivative(A(r), (r, 2))/A(r) + 0.25*Derivative(A(r), r)*Derivative(B(r), r)/(A(r)*B(r)) + 0.25*Derivative(A(r), r)**2/A(r)**2 + 1.0*Derivative(B(r), r)/(r*B(r))[br]R_tt = 0.5*r*Derivative(B(r), r)/B(r)**2 - 0.5*r*Derivative(A(r), r)/(A(r)*B(r)) + 1.0 - 1.0/B(r)[br]R_pp = (0.5*r*A(r)*Derivative(B(r), r) - 0.5*r*B(r)*Derivative(A(r), r) + 1.0*(B(r) - 1)*A(r)*B(r))*sin(t)**2/(A(r)*B(r)**2)[br][/b][br]

計算結果を見やすくまとめて＝０としましょう。偏微分はすべてｒに関するものだけです。[br][math]R_{cc}=\frac{A''}{2B}-\frac{A'B'}{4B^2}-\frac{A'^2}{4AB}+\frac{A'}{rB}=0[/math][br][math]R_r_r=-\frac{A''}{2A}+\frac{A'B'}{4AB}+\frac{A'^2}{4A^2}+\frac{B'}{rB}=0[/math] [br][math]R_t_t=\frac{rB'}{2B^2}-\frac{rA'}{2AB}+1.0-\frac{1.0}{B}=0[/math][br][math]R_{pp}=\frac{\left(0.5r*A*B'-0.5r*B*A'+(B-1)*A*B\right)*sin^2(t)}{A*B^2}=0[/math] [br]最初の2式から、[br][math]\frac{R_{cc}}{A}+\frac{R_r_r}{B}=\frac{A''}{2AB}-\frac{A'B'}{4AB^2}-\frac{A'^2}{4A^2B}+\frac{A'}{rAB}-\frac{A''}{2AB}+\frac{A'B'}{4AB^2}+\frac{A'^2}{4A^2B}+\frac{B'}{rB^2}=\frac{A'}{rAB}+\frac{B'}{rB^2}=\frac{1}{rAB^2}\left(A'B+AB'\right)=0[/math]　[br]積の微分公式を思いだそう。(AB)'=0となるから、AB=定数。[br]時空が先々平坦ならA=1,B=1と仮定してよい。AB=1となるから、B=1/A。[br]これを第３式に代入しよう。[br][math]R_t_t=-\frac{rA'}{2}-\frac{rA'}{2}+1.0-1.0A=-rA'+1-A=0[/math] となる。[br]'r=１なので、(rA)'-1=0とかけるから、rで積分すると、rA-r=C（積分定数）となり、A=1+C/r。[br]ニュートン力学の意味で天体の質量をMとして、[br]たとえば、天体を含む半径ｒの領域の体積積分をしてガウスの定理を使うと、[br]面積分をしてCを-２GM/c^2と求めることができるから、[br]A=[math]1-\frac{2GM}{c^2r}[/math] （B=1/A）となる。[br]これによって、[br]共変計量テンソル [b]g_ij=diag([/b][math]1-\frac{2GM}{c^2r}[/math][b] , [/b][math]-\left(1-\frac{2GM}{c^2r}\right)^{-1}[/math][b] , -r[/b][sup]2[/sup][b], -r[/b][sup]2[/sup][b] sin[/b][sup]2[/sup][b]θ )　[br][/b]が決まるね。[br][br][b][size=150]＜シュワルツシルト半径とブラックホール＞[br][/size][/b][br]A＝０のとき半径ｒ＝[math]\frac{2GM}{c^2}[/math] となる。[br]これを[b]シュワルツシルト半径[/b]と呼ぶ。略してrs[br][br]A=0のとき、g_cc＝０となり、g_rrは逆数で発散する。[br]重力天体の半径rがrsより大きく、真空領域を含まないなら問題ない。[br]重力天体の半径rがrsより小さいとか、大きくても重力分布がrsの内側に集中しているものが、[br]ブラックホールだ。[br][br]ブラックホールの近辺を想像してみよう。[br]ｒがrsに近いときどうなるだろうか。[br]g00の1メモリが無限に小さくなる。進むために[b]時間が膨大にかかり、[/b]遠くから見て止まって見える。[br]grrの1メモリは負の無限大に発散する。[b]半径r方向に空間が無限に引き延ばされ[/b]る。[br]rsは[b]事象の境界線（地平線）[/b]という名前があり、吸い込まれたら見えなく観測できない。[br]しかし、吸い込まれる粒子から見れば連続的に移動しているだけだ。急激だけれどね。

[b][size=150]＜シュワルツシルト半径の具体化＞[/size][/b][br]G＝6.7×10[sup]-11[/sup][m^3/ (kg*s^2)][br]c=3×10[sup]8[/sup][m/s] から、[br]1/c[sup]2[/sup]=1.111*10^-15[br]とするとき、シュワルツシルト半径rsを求めてみよう。rは実半径です。[br][br]太陽の場合[br][br]M＝2×10[sup]30[/sup][kg]、r=6.96×10[sup]8[/sup][m][br]rs=(2G*1/c^2)*M=2*6.7*1.111[sup](-11-15)[/sup]*2×10[sup]30[/sup]=(1.489*10[sup]-27[/sup])*(2*10[sup]30[/sup])=2978[m] ＝約3㎞[br]r/rs=6.96*10**8/2978=230000[br]今の太陽の半径が23万分の１に中心に固まれば、ブラックホール並みになる。[br][br]地球の場合[br][br]M＝6×10[sup]24[/sup][kg]、r=6.38×10[sup]6[/sup][m][br]rs=(2G*1/c^2)*M=(1.489*10[sup]-27[/sup])*(6*10[sup]24[/sup])=8.9*10[sup]-3[/sup][m] ＝約9ミリ[br]3.38*10**6/0.0089=380000000[br]今の地球の半径が38000万分の１に中心に固まれば、ブラックホール並みになる。