[size=150][b]このワークシートは[url=https://www.geogebra.org/m/twxxx3yq]Math by Code[/url]の一部です。[br][/b][/size][size=150][size=100][br]アインシュタイン方程式は、[b]非線形[/b]で[b]座標系[/b]が任意で[b]解釈[/b]の自由度もあるという、[br][color=#0000ff][b]ゆるく複雑さの高いもの[/b][/color]です。[br]条件をつけて解釈することで、宇宙にまつわる様々な事象を予測する道具にすることができます。[br]前回は、ブラックホールの存在予測につながったシュワルツシルト解について学びました。[br][br][/size][/size]今回は、宇宙項があって、ふくらむアインシュタイン方程式の解をさぐろう。
[b]条件1・一様で等方的空間。[br][/b][color=#0000ff][b]条件2・宇宙流体は静止しているので、エネルギーテンソルT[/b][sup]ij[/sup][b]=diag(ρ,p,p,p)[br][/b][/color][b]条件3・アインシュタインの方程式は[br] R[sub]ij[/sub] ー1/2g[sub]ij[/sub] R([/b]ー[b]Λ[/b][color=#0000ff][b]g[sub]ij[/sub] )[/b][/color][b]=[color=#0000ff]8πGT[/color][color=#0000ff][sub]ij[/sub]/c[sup]4[/sup][/color][/b][b]。[br]条件4・宇宙は膨張・縮小するとしても一様に伸縮する。[br][/b]とします。[br][br][b]条件1[/b]から、球面座標系[b]をxi=(r, [/b]θ, φ)[b]とする。[/b][br]空間部分の計量テンソル は[b]g_ij=diag(g(r), r[sup]2[/sup], r[sup]2[/sup] sin[sup]2[/sup]t ),[br][/b][b]g^ij=diag(1/g(r),1/r[/b][sup]2[/sup][b], 1/r[/b][sup]2[/sup][b] sin[/b][sup]2[/sup]θ[b] )[/b]とおけるね。[br][br][color=#9900ff][b][u][size=150]課題:等方的空間のリッチテンソルをもとめよう。[br][/size][/u][/b][/color][br]シュワルツシルト半径を求めたときのように、リッチテンソルを求めるコードを実行すればよいね。import sympy as sp[br]from sympy import symbols, sin, diff, Matrix, Function, simplify[br][br]# 変数と関数の定義[br][b]r, t, p = symbols('r t p')[br]xs = [r, t, p][br]dim = len(xs)[br]g = Function('g')(r)[br][/b][br]# 計量テンソル g_mn[br]M = sp.diag(g, r**2, r**2 * sin(t)**2)[br]M_inv = M.inv()[br]# 辞書を使って gamma[i][j][k] = Γ^i_jk の形で保存[br]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] # 第1種 Γ_{m,jk} = 1/2 * (g_{mj,k} + g_{mk,j} - g_{jk,m})[br] gamma1_mjk = 0.5 * (diff(M[m, j], xs[k]) + diff(M[m, k], xs[j]) - diff(M[j, k], xs[m]))[br] # 第2種 Γ^i_jk = g^im * Γ_{m,jk}[br] 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] # 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] # 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] # R_ij = Σ R^m_imj[br] res_ricci = 0[br] for m in range(dim):[br] res_ricci += riemann[m][i][m][j][br] [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][OUT][br]計算中: クリストッフェル記号...[br]計算中: リーマンテンソル R^i_jkl ...[br][br]計算中: リッチテンソル R_ij ...[br]計算中: リッチテンソル R_ij ...[br]R_rr = 1.0*Derivative(g(r), r)/(r*g(r))[br]R_tt = 0.5*r*Derivative(g(r), r)/g(r)**2 + 1.0 - 1.0/g(r)[br]R_pp = (0.5*r*Derivative(g(r), r) + 1.0*(g(r) - 1)*g(r))*sin(t)**2/g(r)**2[br]となる。[br]つまり、[br][br][math]R_r_r=\frac{g'}{r*g}[/math][br][math]R_t_t=\frac{rg'}{2g^2}+1-\frac{1}{g}[/math] [br][math]R_{pp}=\frac{rg'}{2g^2}+(1-\frac{1}{g})sin^2(t)[/math][br]計量テンソルの対角成分以外は0なので、この3つのテンソルに対応する対角成分をかけた積和が[br]スカラー曲率になるね。[br]R=simplify(M_inv[0,0]*ricci[0][0] + M_inv[1,1]*ricci[1][1] + M_inv[2,2]*ricci[2][2])[br]print(R)[br]#================================================[br][OUT][br]2.0*Derivative(g(r), r)/(r*g(r)**2) + 2.0/r**2 - 2.0/(r**2*g(r))[br]これをさらにカンタンにすると、2.0/r**2{1+ rg'/g[sup]2[/sup]-1/g}[br]=[math]\frac{2}{r^2}\left\{1-\left(\frac{r}{g}\right)'\right\}[/math] =C(定数)とすると、曲率が一定になる。[br][b]条件1[/b]からこの曲率が定数ならば、一様で等方的な空間になる。[br]積分して、g=1/(1-C/6 * r^2+ const/r) [br]原点を任意にできるようにするためと、Cとrの次元にそろえるために次のようにおこう。[br]g=1/(1-kr[sup]2[/sup]) ただし、k=0,1,-1のどれかにする。[br]これで、[b]条件1[/b]の一様で等方的な空間の設定になった。[br][br]空間部分の計量テンソル は[br][b]g_ij=diag(1/(1-k r[sup]2[/sup]), r[sup]2[/sup], r[sup]2[/sup] sin[sup]2[/sup]t )[br]r基底のメモリに着目しよう。[br][/b]k=0ではメモリは常に1だから、[b]平坦空間[/b]になる。[br]k=-1ではメモリは必ず正だが、原点から遠くなるとメモリが巨大数分の1となり、[br]微細になる。だから、逆に距離は極度に大きくなり、[b]果てのない開いた宇宙[/b]になる。[br]k=1ではメモリrは1以内でしかありえないために、[b]閉じた球状宇宙[/b]となる。
[b][color=#0000ff]条件1から設定したkは宇宙の形を決めた。[br]次は、条件4・宇宙は膨張・縮小するとしても一様に伸縮する。[br]を考えてみよう。[br][/color][/b]条件5から、空間の計量の線素をdl[sup]2[/sup]とすると、dlをa倍する相似係数(宇宙のスケール)を考えて、[br]しかも、これが時間とともに変わる関数a(t)としたい。[br]その4次元の一般座標の線素は[br]ds[sup]2[/sup]=(cdt)[sup]2[/sup]-(dl)[sup]2[/sup]=c[sup]2[/sup]dt[sup]2[/sup]-a[sup]2[/sup](t)[[b]1/(1-k r[sup]2[/sup]) dr[sup]2[/sup]+r[sup]2[/sup]dθ[sup]2[/sup]+ r[sup]2[/sup] sin[sup]2[/sup]tdφ[sup]2[/sup]] [br][/b]とおける。[br]これが、[b][color=#0000ff]フリードマン・ルメートル・ロバートソン・ウォーカ計量(FLRW計量)[/color][/b]です。[br][size=150][b][u][color=#9900ff]課題:FLRW計量のリッチテンソルとスカラー曲率を求めるにはどうしたらよいですか。[br][/color][/u][/b][/size][br]変数名がかぶりやすいので、今まで1文字でtheta,phiをt,pとしていたものを略さずにかきます。[br]その代わり、時間の基底変数はtではなくctにしましょう。[br]自然単位系ではないので、光速はcをつかいますが、微分の関係で時間の計量を1にして、[br]ct自体を1つの変数として微分しておき、最後にtで微分することで、時間微分表示に直しましょう。[br]微分の連鎖率により、[br]# dot{a}=da(ct)/dt= da(ct)/dct * dct/dt=a'*c[br]# ddot{a}=d(a'*c/dt)=d(a'(ct)/dt*c)+a'*dc/dt=a''*c2[br]# だから、[br]# a'(ct) = (1/c) * dot{a}, a''(ct) = (1/c^2) * ddot{a}[br]# 表示を分かりやすくするため、a1=dot{a}, a2=ddot{a} として代入表示[br]というのがそのやり方です。[br][br]import sympy as sp[br]from sympy import symbols, sin, diff, Function, simplify[br][br]# 1. 変数と座標の定義[br]# 微分は文字通り (ct) で行う[br]ct, r, theta, phi = symbols('ct r theta phi')[br]c, k = symbols('c k')[br]xs = [ct, r, theta, phi][br]dim = len(xs)[br]a = Function('a')(ct)[br][br]# 2. 標準的な FLRW 計量テンソル (Signature: + - - -)をそのまま使う。[br]# ds^2 = 1(dct)^2 - a(ct)^2 [ dr^2/(1-kr^2) + r^2 dθ^2 + r^2 sin^2θ dφ^2 ][br]M = sp.diag(1, [br] -a**2 / (1 - k * r**2), [br] -a**2 * r**2, [br] -a**2 * r**2 * sin(theta)**2)[br]M_inv = M.inv()[br][br]# 3. クリストッフェル記号の計算[br]gamma = {}[br]for i in range(dim):[br] gamma[i] = {}[br] for j in range(dim):[br] gamma[i][j] = {}[br] for l in range(dim):[br] res = 0[br] for m in range(dim):[br] g1 = 0.5 * (diff(M[m, j], xs[l]) + diff(M[m, l], xs[j]) - diff(M[j, l], xs[m]))[br] res += M_inv[i, m] * g1[br] gamma[i][j][l] = simplify(res)[br][br]# 4. リッチテンソル R_ij の計算[br]print("非ゼロのricciテンソル表示:")[br]ricci = [][br]for i in range(dim):[br] row = [][br] for j in range(dim):[br] res_ricci = 0[br] for m in range(dim):[br] # R_ij = ∂_m Γ^m_ij - ∂_j Γ^m_im + Γ^m_mn Γ^n_ij - Γ^m_jn Γ^n_im[br] term1 = diff(gamma[m][i][j], xs[m])[br] term2 = diff(gamma[m][i][m], xs[j])[br] term3 = sum(gamma[m][m][n] * gamma[n][i][j] for n in range(dim))[br] term4 = sum(gamma[m][j][n] * gamma[n][i][m] for n in range(dim))[br] res_ricci += term1 - term2 + term3 - term4[br] row.append(simplify(res_ricci))[br] if res_ricci != 0:[br] simple_ricci = res_ricci.subs({diff(a, ct, 2): a2/c**2, diff(a, ct): a1/c})[br] print(f"R_{xs[i]}{xs[j]} = {simple_ricci}")[br] ricci.append(row)[br][br]# 5. スカラー曲率 R = g^ij R_ij の計算[br]R_val = 0[br]for i in range(dim):[br] R_val += M_inv[i, i] * ricci[i][i][br][br]R_final = simplify(R_val)[br][br]#print("--- 導出結果 ---")[br]#print(f"スカラー曲率 R (ct微分版): {R_final}")[br]# 6. 時間微分 (t) への書き換え[br]# dot{a}=da(ct)/dt= da(ct)/dct * dct/dt=a'*c[br]# ddot{a}=d(a'*c/dt)=d(a'(ct)/dt*c)+a'*dc/dt=a''*c2[br]# だから、[br]# a'(ct) = (1/c) * dot{a}, a''(ct) = (1/c^2) * ddot{a}[br]# 表示を分かりやすくするため、a1=dot{a}, a2=ddot{a} として代入表示[br]a1, a2 = symbols('a1 a2')[br]R_dot_notation = R_final.subs({diff(a, ct, 2): a2/c**2, diff(a, ct): a1/c})[br]print(f"スカラー曲率 R (時間微分ドット版): {simplify(R_dot_notation)}")[br] #=====================================[br][OUT][b]非ゼロのricciテンソル表示:[br][/b]R_ctct = -3.0*a2/(c**2*a(ct))[br]R_rr = -2.0*a1**2/(c**2*(k*r**2 - 1)) - 1.0*a2*a(ct)/(c**2*(k*r**2 - 1)) - 2.0*k/(k*r**2 - 1)[br]R_thetatheta = 2.0*a1**2*r**2/c**2 + 1.0*a2*r**2*a(ct)/c**2 + 2.0*k*r**2 - 1.0*(-tan(theta)**2 - 1)/tan(theta)**2 - 1.0 - 1.0/tan(theta)**2[br]R_phiphi = 2.0*a1**2*r**2*sin(theta)**2/c**2 + 1.0*a2*r**2*a(ct)*sin(theta)**2/c**2 + 1.0*k*r**2*sin(theta)**2 + 1.0*(k*r**2 - 1)*sin(theta)**2 + 0.5*sin(2*theta)/tan(theta) - 1.0*cos(2*theta)[br][b]スカラー曲率 R (時間微分ドット版): 6.0*(-a1**2 - a2*a(ct) - c**2*k)/(c**2*a(ct)**2)[br][/b][br]時間微分はダッシュでなくドットで表示しますが、geogebraの表記機能の関係からダッシュにします。 [br] [br][math]R_t_t=-3.0*a''/c^2a[/math] [br][math]R=-6.0\frac{c^2*k+a*a''+a'^2}{c^2a^2}[/math] [br][br]これらを見ると、[br]時間基底での曲率は「膨張の加速度a''」に連動していることがわかります。[br]スカラーの曲率は、「膨張の速度a'、その加速度a''」さらに「宇宙の形係数k」の混ざり合い[br]になってますね。
リッチテンソルRttとスカラー曲率Rが求められたから、[br]アインシュタイン方程式に代入できるね。[br][br][color=#0000ff]条件2・宇宙流体は静止しているので、エネルギーテンソルT[sup]ij[/sup]=diag(ρ,p,p,p)[br](ρはエネルギー密度、pは運動量。)[br][/color]条件3・アインシュタインの方程式は[br] R[sub]ij[/sub] ー1/2g[sub]ij[/sub] RーΛ[color=#0000ff]g[sub]ij[/sub] [/color]=[color=#0000ff]8πGT[/color][color=#0000ff][sub]ij[/sub]/c[sup]4[/sup][/color]。[br]条件3の時間座標から、R[sub]00[/sub]-1/2g[sub]00[/sub]R-Λ[color=#0000ff]g[sub]00[/sub] [/color]=[color=#0000ff]8πGT[/color][color=#0000ff][sub]00[/sub]/c[sup]4[/sup][/color][br]時間座標ctの成分だけで、代入してみよう。[br][br]・g[sub]00[/sub]=M[sub]00[/sub]=1[br]・条件2の時間座標T[sub]00[/sub]=ρ[br]条件1,5から得た、[br]・[math]R_{00}=-3.0*a''/c^2a[/math] [br]・[math]R=-6.0\frac{c^2*k+a*a''+a'^2}{c^2a^2}[/math] [br]を入れると、[br][br]アインシュタイン方程式は、次のようになるね。[br][math]-3.0*a''/c^2a-\frac{1}{2}\ast\left(-6.0\frac{c^2*k+a*a''+a'^2}{c^2a^2}\right)-Λ=\frac{8\pi G}{c^4}\rho[/math] [br]式を整理すると[br][math]\frac{c^2*k+a'^2}{c^2a^2}-Λ/3=\frac{8\pi G}{3c^4}\rho[/math] [br]両辺にc²をかけて、[br][math]\frac{c^2*k+a'^2}{a^2}-Λc^2/3=\frac{8\pi G}{3c^2}\rho[/math] [br]これが[b]フリードマン方程式[/b]だ。[br][br]膨張の加速度a''が消えました。[br]T[sub]00[/sub]=ρエネルギー密度ではなく、T[sup]00[/sup]=ρc²質量密度として、単位系を整理して、[br]K=c[sup]2[/sup]kとおこう。Kを宇宙曲率と呼ぶことがある。[br][math](\frac{a'}{a})^2+\frac{K}{a^2}=\frac{8\pi G}{3}\rho+\frac{Λ}{3}[/math] 第1方程式(宇宙のエネルギー)[br]これが[b]宇宙膨張のフリードマン方程式[/b]だ。[br][br]これを時間で微分すると、[br][math]\frac{a''}{a}=\frac{4\pi G}{3}(\rho'\frac{a}{a'}+2p)+\frac{Λ}{3}[/math] [br]計算の詳細は省略するが、Tijについても、同様に代入と整理をすることで、[br][math]\frac{a''}{a}=-\frac{4\pi G}{3}(\rho+3p)+\frac{Λ}{3}[/math]が得られる。第2方程式(宇宙空間の加速度)。[br]pは運動量、つまり、圧力pだ。圧力pのあるもの(光、ガス)が多いほど宇宙の膨張に負の影響、ブレーキがかかるということだね。[br][br]この2式の左辺が同じだから、比べることから次の式ができるね。[br][math]\rho'=-3\frac{a'}{a}(\rho+p)[/math]が得られる。[br]しかし、圧力pと質量エネルギーρの関係式がこの1個では、宇宙の状態構造は決まらない。[br][br][size=150][b]<宇宙状態方程式>[/b][br][/size]圧力pと質量エネルギーρの関係式、[br]状態方程式p=ωρによって、宇宙状態方程式ができるでしょう。[br]代表的な例をあげよう。[br][b]・p≪ρ(非相対論的物質中心モデル)[br][/b] 光、ガスによる圧力が少なく、塵や星の重力のみで宇宙膨張にブレーキがかかる宇宙。[br] ω=0つまり、p=0の場合の宇宙状態で、ρ∝a[sup]-3[/sup]。[br] [color=#0000ff][b]K=0、Λ=0は平坦で非相対論的物質で満ちた宇宙モデル(アインシュタインドジッター宇宙)[/b][/color]がある。[br][b]・p=ρ/3(相対論的物質モデル)[br][/b] ω=1/3の場合。[br] 光、ガスの圧力による体積増加の他に、波長の伸びでも宇宙膨張にブレーキがかかる宇宙。[br] m=0でもρ∝a[sup]-4[/sup]もありうる。[br][b]・p=-ρ(真空モデル)[br][/b] ω=-1の場合。[b]ρが一定で[/b]膨張しても密度が減らない。[br] 物質がなく宇宙定数Λだけが存在し、[color=#0000ff][b]指数関数的に膨張する[/b][/color][color=#0000ff][b]宇宙([/b][b]ジッター宇宙)[br][/b][/color][br]一般には、K、Λともに非ゼロで、非相対論的物質と相対論的物質の両方が存在して存在比をパラメータ化[br]したものがフリードマン・ルメートル宇宙。これは固定ではなく可変な宇宙像だ。[br]パラメータそのものを宇宙の成長とともに変えることで、宇宙の歴史を記述しようとするのが、今の標準宇宙論になっているようだ。