リーマンからアインシュタインへ
このワークシートはMath by Codeの一部です。
アリも人間も平坦な道をまっすぐすすんでいるときは、世界はまっすぐでできている
と感じているでしょう。でも、その距離が伸びていくと地面の曲がりを補正する必要が出てきます。
今回は、曲面の曲がりを見失うという局所的な視点を抜け出して、曲面の曲がり具合を読み取り、
大局的な物理量を知るための数学の道具（曲面での微分）を探ってみよう。
2次元曲面を平行移動するベクトル
１．クリストッフェル記号に親しもう
複雑な仕組みを学ぶときは、一般や原理から入るよりも道具や具体から入る方が
イメージがわきやすくなることがある。
そこで、ベクトル、テンソルの変化を調べるための便利な道具、偏微分を使った関数から入ってみよう。
ベクトル解析では、grad,div,rotなど、ベクトルの成分を変数varsで偏微分する関数があったね。
コードでかくと、
import numpy as np
from sympy import symbols, diff,Matrix
f,x,y,z,r,t ,u,v,vx,vy,vz= symbols('f x y z r t u v vx vy vz')
def grad(w,vars): #スカラー場ｗの勾配はベクトル
return Matrix([diff(w, i) for i in vars])
def div(v,vars): #ベクトル場ｖの発散はスカラー
return sum([diff(v[i],vars[i]) for i in range(len(vars))])
def laplacian(w,vars):#ラプラシアンはスカラー場からスカラーへ
return div(grad(w,vars),vars)
def rot(F, vars):#回転はcurlとすることが多い。日本ではrot
Fx,Fy,Fz = F
x,y,z = vars
return Matrix([diff(Fz,y) - diff(Fy,z),diff(Fx,z) - diff(Fz,x),diff(Fy,x) - diff(Fx,y)])
vars= [x,y,z]
などがあった。
テンソル解析でも、テンソルGの成分を変数varsで偏微分することは大切だ。
その土台になる関数は
第１種クリストッフェル記号Γijk、略してクリ１、コード上はChr1(i,j,k,G)
第２種クリストッフェル記号Γijk、略してクリ２、コード上はChr2(i,j,k,G)
の２つがあるね。
＜クリ１＞
クリ１の対象となるテンソルは2軸で、型は(2,2),(3,3),....などの正方行列になる。
そのサイズ2,3,....はdim=len(vars)などとして、可変にしておく。
たとえば、dim=3の場合、vars=[x,y,z]でも、vars=[x1,x2,x3]でも、vars=[r,θ,φ]でも同じことだ。
クリ1自身は3軸あるので、Chr1(i,j,k,G)の3軸i,j,kともにサイズがdimとなるね。
ということは、クリ１の指定は全部で、dimの3乗通りあることなる。
テキストでは
Γijk＝1/2(∂/∂xi(G_jk)+∂/∂xj(G_ki)-∂/∂xk(G_ij))
のように書いてあることが多い。
例を作るとわかりやすくなる。
たとえばvars=[x1,x2,x3]とすると、i,j,kは1スタートの番号1，2，3のどれかという意味でしかない。
i,jが同じなったり、i,j,kともに同じ数字にあっても全然問題ない。
Chr1(1,2,3,G)=1/2(diff(G[2,3],x1)+diff(G[3,1],x2)-diff(G[1,2],x3))
Chr1(1,2,2,G)=1/2(diff(G[2,2],x1)+diff(G[2,1],x2)-diff(G[1,2],x2))
Chr1(2,2,2,G)=1/2(diff(G[2,2],x2)+diff(G[2,2],x2)-diff(G[2,2],x2))
G_mnは計量テンソルなので、基底ベクトルの内積だから、対称行列になるね。
すると、クリ１の定義がi,j,kのサイクリックな式になっているけれど、G_mnのmnは入れ替えても同じだ。
さらに、G_mnが直交座標系の計量なら、m≠nのときは０になる。
これを使うと、
Chr1(1,2,3,G)=1/2(diff(G[2,3],x1)+diff(G[3,1],x2)-diff(G[1,2],x3))=1/2(0+0+0)=0
Chr1(1,2,2,G)=1/2(diff(G[2,2],x1)+diff(G[2,1],x2)-diff(G[1,2],x2))=1/2(diff(G[2,2],x1)
Chr1(2,2,2,G)=1/2(diff(G[2,2],x2)+diff(G[2,2],x2)-diff(G[2,2],x2))=1/2(diff(G[2,2],x2)
となるね。
これら、3の3乗＝27個のデータのかたまりクリ１は3次元データなので、プログラミング上はテンソルと言えるけれど、数学的にはテンソル則（積、商、変換）が適用できないので、テンソルとは言いません。
座標系の基底の関係を表す「計量テンソル」を網羅的に偏微分することで、「座標系のメモリのゆがみ、勾配」をデータ化できてます。
番号の位置に着目すると、先頭の２つの添え字は交換可能だね。Γijk=Γjik
課題：クリ１を求めるコードはどうやればできますか。
i,j,kがvarsの順番をさす、1スタートの数だけれど、配列のインデックスは０スタートなので、
順番がｎのときvars[n-1]となることに注意しよう。
また、次元を可変にするために、dim=len(vars)などを使うといいね。
[IN]Python
#=====================================================
import sympy as sp
from sympy import symbols, sin, cos, Matrix, simplify, diag,diff
#第一種クリストッフェル記号
def Chr1(i,j,k,M):
i=i-1;j=j-1;k=k-1
di = diff(M[j,k],xs[i])
dj = diff(M[k,i],xs[j])
dk = diff(M[i,j],xs[k])
res = 1/2*(di+dj-dk)
#Mはメトリックテンソルg_mnの行列
return res
# 変数
x, y, z, u, v = symbols('x y z u v')
xs = Matrix([x,y,z])
M = Matrix([[1,0,0],[0,x**2,0],[0,0,x**2*(sin(y))**2]])
print("球座標の計量テンソルg_mn:")
display(M)
dim = len(xs)
K_ijkM=[f"Γ{i}{j}{k} = {Chr1(i,j,k,M)}" for i in range(1,dim+1) for j in range(1,dim+1) for k in range(1,dim+1)]
for item in K_ijkM:
print(item)
[OUT]
#===================================================
Γ111 = 0
Γ112 = 0
Γ113 = 0
Γ121 = 0
Γ122 = 1.0*x
Γ123 = 0
Γ131 = 0
Γ132 = 0
Γ133 = 1.0*x*sin(y)**2
Γ211 = 0
Γ212 = 1.0*x
Γ213 = 0
Γ221 = -1.0*x
Γ222 = 0
Γ223 = 0
Γ231 = 0
Γ232 = 0
Γ233 = 1.0*x**2*sin(y)*cos(y)
Γ311 = 0
Γ312 = 0
Γ313 = 1.0*x*sin(y)**2
Γ321 = 0
Γ322 = 0
Γ323 = 1.0*x**2*sin(y)*cos(y)
Γ331 = -1.0*x*sin(y)**2
Γ332 = -1.0*x**2*sin(y)*cos(y)
Γ333 = 0
＜クリ２＞
クリ１と計量テンソルG_mnの逆行列G^mnを使った積との和、微分の積和の構造によって、
「変化の集計、ズレの大きさ」という視点で座標系のゆがみをつかまえるためにクリ２があります。
クリ２もdimの3乗個のデータでできてます。その点から別の点へいくときの基底と基底のつながり・ゆがみを求めているので、「テンソルではなく接続（係数）」と呼ばれています。
テキストでは
Γijk＝Σr gir (Γjkr)
dim=3ならば、
Γijk＝gi1 Γjk1+gi2 Γjk2+gi3 Γjk3ということです。
クリ２の下つきの２添え字は、クリ１の先頭の２添え字と同じなので、入れ替え可能だとわかるね。
つまり、Γijk=Γikj
課題：クリ２を求めるコードはどうやればできますか。
dim=len(vars)などを使うといいね。
クリ１とセットにして使えば、Pythonのカーネルをクリアしたあとに、
部分的にセルが使えるね。
[IN]Python
import sympy as sp
from sympy import symbols, sin, cos, Matrix, simplify, diag,diff
#第2種Chr2(i,j,k,M)=Γ^i_jkは第1種Chr(j,k,i,M)に対して3番目の変数iを使ったg^ir*Chr(j,k,r,M)のrを動かした和になる。
#第一種クリストッフェル記号
def Chr1(i,j,k,M):
i=i-1;j=j-1;k=k-1
di = diff(M[j,k],xs[i])
dj = diff(M[k,i],xs[j])
dk = diff(M[i,j],xs[k])
res = 1/2*(di+dj-dk)
#Mはメトリックテンソルg_mnの行列
return res
#第2種クリストッフェル記号
def Chr2(i,j,k,M):
g_ir = M.inv()
res=0
#行列には０スタートで渡し、記号には１スタートで渡す
for r in range(dim):
res += g_ir[i-1,r]*Chr1(j,k,r+1,M)
return res
# 変数
x, y, z, u, v = symbols('x y z u v')
xs = Matrix([x,y])
M = Matrix([[1,0],[0,x**2]])
print("極座標の計量テンソルg_mn:")
display(M)
dim = len(xs)
K2_ijkM=[f"Γ^{i}_{j}{k} = {Chr2(i,j,k,M)}" for i in range(1,dim+1) for j in range(1,dim+1) for k in range(1,dim+1)]
for item in K2_ijkM:
print(item)
#====================================================
[OUT]
極座標の計量テンソルg_mn:
Γ^1_11 = 0
Γ^1_12 = 0
Γ^1_21 = 0
Γ^1_22 = -1.0*x
Γ^2_11 = 0
Γ^2_12 = 1.0/x
Γ^2_21 = 1.0/x
Γ^2_22 = 0
２．クリストッフェルから共変微分へ
＜ベクトルの平行移動＞
ベクトルがある世界に制限がなければ、平行移動してもベクトルは変わらないね。
でも、ベクトルの始点を時空の世界に制限すると、平行移動の意味が変わる。
たとえば、地球の表面をはなれないようにベクトルAをxからまで(ｘ＋dx)まで進む。
表面だけの世界では２つのベクトルは「平行」と見られるので、平行成分をA||とかく。
しかし、始点を合わせるとA(x)とずれが終点にできる。この垂直な補正成分をA⊥とかく。
つまり、細かくかくと、A(x)=表面にそったベクトル＋ベクトルの終点のズレベクトル
=A||(x->x+dx)＋A⊥(x->x+dx)となるね。
これは、基底ベクトルemをｘ＋dxからｘに戻したときのようすを図にかくと、
em(x + dx) = e||m(x + dx → x) （平行成分）+ e⊥m(x + dx → x)（垂直成分）
表面が曲がっているので、ベクトルはもとのベクトルからずれる。
それを調整するためのものが必要になる。
基底ベクトル em を xk 方向に動かしたとき、それは元の基底 enたちの線形結合でどう表されるか？
e||m(x + dx → x) = em(x) + dxk ∂/∂xk(em)
= em(x) + dxk Γnkm en(x)
とできる。
Γnkmを（アフィン）接続係数という。
これはさっきやった第２種クリストッフェル記号、クリ２だね。
クリ２は3軸で、3次元ならば３つの 基底en(x)それぞれにかけ算する係数になっている。
3次元ならば、3つのどの基底を、３つのどの方向に動かし、それが３つのどの基底成分に影響するのか？
Γnkmはそのすべての27個の係数により、あらゆる『ねじれ』のパターンを網羅している、基底の線形結合の係数になっている。だから、接続係数と呼ばれているんだね。
ｎは影響を受ける基底、m は動かす基底、ｋが動かす方向を表す。
変数xkに関する偏微分の演算子を
∂k≡: ∂/∂xkとかくことがあります。
Amがスカラーなら、xkに関する偏微分は∂kAm
Amがベクトルだと、ベクトルの基底自体が変化する。
だから、さっきベクトルの平行移動での基底の微分のときのように補正が必要になる。
ベクトルの普通の偏微分の他に、補正をした微分を共変微分という。
＜反変ベクトル＞
ベクトルA=Amemをxk方向に微分すると、積の微分法則から、
∇kA＝(∇kAm)em+Am(∇kem ) =(Am,k+ΓmknAn)em =Am,kem
成分だけかくと、Am,k ＝ Am,k+ΓmnkAn
反変ベクトルの共変微分という。
つまり、座標系（ｘm)の反変ベクトルA＝(An)のｘkに関する共変微分はテンソルとなる。
▽k Am＝∂kAm ＋ Γmnk An
反変ベクトルの共変微分をとする書き方もありますね。
Am；k≡: ▽k Am
この場合は、反変ベクトルの普通の偏微分を
Am,k≡: ∂k（Am）ともかいて区別したりします。
しかし、テキストによっては、微分も共変微分も両方Am,kで表すものもある。文脈で判断する。
では、実際にどのように計算するのでしょうか。第２項の共通軸はnですから3次元ならば、n=1,2,3です。
▽k Am＝∂kAm ＋Γm1k A1＋Γm2k A2＋Γm3k A3となり、㎞を指定した微分になります。
クリ２でm番の基底が、ｋ番目の基底をベクトルの成分を1，2，3の方向に動かした影響を加えています。
その調整によって、座標系のゆがみが取り除かれて、座標変換できる物理量、テンソルになったのです。
＜共変ベクトル＞
Aν と Bν を共にベクトルとすると Aν Bν はスカラーだから∇µ(Aν Bν) = ∂µ(Aν Bν)。
左辺＝(∂µAν + ΓνλµAλ)Bν + Aν ∇µBν
右辺＝ (∂µAν )Bν + Aν ∂µBν
式変形により、∇µBλ = ∂µBλ − ΓνλµBν
つまり、座標系（ｘm)の共変ベクトルA＝(Am)のｘkに関する共変微分はテンソルとなる。
▽k Am＝∂kAm ー Γnmk An
では、実際にどのように計算するのでしょうか。第２項の共通軸はnですから3次元ならば、n=1,2,3です。
▽k Am＝∂kAm ー(Γ1mk A1＋Γ2mk A2＋Γ3mk A3)となり、㎞を指定した微分になります。
クリ２でベクトルの基底１，２，３が、k番の基底をmの方向に動かした影響を差し引いています。
課題：共変微分をコードで実現するにはどうしたらよいでしょうか。
ベクトルの普通の微分に対して、クリ２の合計をたしひきするようにすればよいですね。
上記のクリ１、クリ２をかいたセルに下におくか、
同じプログラムコードの中におきましょう。
[IN]Python
#============================================
def CovariantDiff_Up(m, k, A_up, M):
"""反変ベクトル A^m の座標 k による共変微分"""
# 1. 普通の偏微分
res = diff(A_up[m-1], xs[k-1])
# 2. クリストッフェル項（接続による補正）を足す
for n in range(1, dim + 1):
res += Chr2(m, n, k, M) * A_up[n-1]
return simplify(res)
def CovariantDiff_Down(m, k, A_down, M):
"""共変ベクトル A_m の座標 k による共変微分"""
# 1. 普通の偏微分
res = diff(A_down[m-1], xs[k-1])
# 2. クリストッフェル項（接続による補正）を引く
for n in range(1, dim + 1):
res -= Chr2(n, m, k, M) * A_down[n-1]
return simplify(res)
３．曲率のテンソル
＜絶対微分から測地線へ＞
共変微分は局所的な点の近くでの微分でした。
それを曲線にそって動かしたときの微分を絶対微分といいます。
絶対微分にも、ふつうの微分に加えて、クリ２を使った補正項がつきますね。
反変ベクトルTiのrによる共変微分Ti;rと曲線C：xi=xi(t)の接ベクトル(dxr/dt)との内積
（Ti;r・dxr/dt)はTiと同じ型のテンソルで、曲線Cにそう絶対微分といい、
Ti=Ti(x(t))とするとき、（δTi/δt)≡：（dTi/dt + Γirs Tr dxs/dt) とかく。
たとえば、直交座標系の加速度ベクトルは、位置ベクトルｘ＝{xi(t)}の2階微分d2xi/dt2。
ｔにおける距離はa=sqrt(δij aiaj)。
一般座標の加速度は、 a=
おもしろいことに、tを弧長sに変えた式が曲率を表す。
一般座標の曲率は、 b=
その時空で曲がってない線が測地線(その時空の直線）だから、
測地線の方程式は曲率＝０となるはず。
つまり、 （i=1,2,...,n）
＜曲率テンソル＞
さて、曲率と言えば外せないものがあります。
それは、リーマンテンソル（リーマンの曲率テンソル）です。
ここでも、クリ２が大活躍しますよ。
リーマン幾何学では、大局的な平行の概念は存在しない。
リーマンの曲率テンソルは、
Rnm,kl = (∂kΓnlm + Γnki Γilm) − (∂lΓnkm + Γnli Γikm)
4軸なので、4元では４の4乗の256個のデータで時空の曲がりを表現しているのです。
その重いデータをコンパクトにしたものがあります。
まず、曲率テンソルを次のように縮約したテンソル、いわゆるリッチテンソルがあります。
Rml = Rkm,kl
まだ2軸なので、４の2乗で16個のデータがあります。
さらに、これを縮約したものがスカラー曲率です。1個のデータです。ちょっと軽すぎか？
R = gml Rml = gml Rkm,kl
さらに、アインシュタインはこの材料から美味しいテンソルを作りました。
Gij = Rij −1/2gij R
これがアインシュタインテンソルです。
種類が多いですが、流れとしてはこうです。
リーマンで全方向しらべたものをリッチテンソルで１方向につぶし、スカラーで数値にして単純化した。
最終的にアインシュタインがエネルギー保存則（ΔGij＝０）が成り立つように、曲率を集大成した。
ということです。
アインシュタインテンソルという幾何学的な量が不変であることが、
エネルギーの保存則を表しているというは、おもしろいね。
計算は省略しますが、ビアンキの恒等式(Rijkl;m+Rijlm;k+Rijmk;l=0)を2回縮約して、計量テンソルの共変微分が０という事実を使うと、アインシュタインテンソルの発散が０になることがわかる。