Python（SymPy）を使って学ぶアインシュタイン方程式のライスナー=ノルドシュトロム解

今回は Python（SymPy）を活用する具体例として、４次元時空のReissner-Nordstrom（ライスナー=ノルドシュトロム）解が、電場が存在する場合のアインシュタイン方程式の解になっていることをPython（SymPy）を使って確認してみたいと思います。電場が存在しない場合には、前の記事

pianofisica.hatenablog.com

でみたように、宇宙定数がゼロの場合にはSchwarzschild（シュワルツシルト）解が、宇宙定数がゼロでない場合にはde Sitter-Schwarzschild（ド・ジッター=シュワルツシルト）解が、静的で球対称な場合のアインシュタイン方程式の真空解として求められました。今回はこれらを少し拡張したものになります。
　

ド・ジッター解
電磁場テンソル
電磁場のエネルギー・運動量テンソル
ライスナー=ノルドシュトロム解

ド・ジッター解

計算の詳細については、上でリンクを貼った過去の記事を参照してください。

宇宙定数 $\Lambda$ を含み、物質場の存在しない場合のEinstein方程式は

$\qquad\displaystyle{G_{ij}+\Lambda\,g_{ij}=0}$

です。その解は

$\qquad\displaystyle{g_{ij}=\left(\begin{array}{cccc} -\left(1-\frac{a}{r}-\frac{\Lambda}{3}r^2\right) & 0 & 0 & 0 \\ 0 & \frac{1}{1-\frac{a}{r}-\frac{\Lambda}{3}r^2} & 0 & 0 \\ 0 & 0 & r^2 & 0 \\ 0 & 0 & 0 & r^2\sin^2\theta \end{array} \right)}$

で与えられます。ただし、これらの成分は座標系

$\qquad\displaystyle{\xi^{i}=\left(\begin{array}{c} t \\ r \\ \theta \\ \phi \end{array} \right)}$

に関するものです。確認しておくと

%reset
import sympy as sp

sp.var('t, r, theta, phi, psi, a, Lambda')
xi = sp.Matrix([ [t],[r],[theta],[phi] ])
D = len(xi)

f = 1 - a/r - (Lambda/3)*r**2
g = sp.Matrix([
[-f, 0, 0, 0],
[0, 1/f, 0, 0],
[0, 0, r**2, 0],
[0, 0, 0, r**2*sp.sin(theta)**2] ])

ginv = g.inv()

def Gamma(k,i,j):
    gm = 0
    for n in range(D):
        gm = gm + sp.Rational(1,2)*ginv[k,n]*(g[n,i].diff(xi[j])\
                +g[n,j].diff(xi[i])-g[i,j].diff(xi[n]))
    return gm

def R(i,j,k,l):
    rm = Gamma(i,l,j).diff(xi[k]) - Gamma(i,k,j).diff(xi[l]) 
    for n in range(D):
        rm = rm + Gamma(n,l,j)*Gamma(i,k,n) - Gamma(n,k,j)*Gamma(i,l,n)
    return rm

def Ricc(j,k):
    ric = 0
    for n in range(D):
        ric = ric + R(n,j,n,k)
    return ric

RiccMat = sp.Matrix([
[Ricc(0,0), Ricc(0,1), Ricc(0,2), Ricc(0,3)],
[Ricc(1,0), Ricc(1,1), Ricc(1,2), Ricc(1,3)],
[Ricc(2,0), Ricc(2,1), Ricc(2,2), Ricc(2,3)],
[Ricc(3,0), Ricc(3,1), Ricc(3,2), Ricc(3,3)] ])

R = sp.simplify((ginv*RiccMat).trace())

GMat = sp.simplify(RiccMat - sp.Rational(1,2)*g*R)

EMat = sp.simplify(GMat + Lambda*g)
EMat

より、確かにアインシュタイン方程式の解になっていることがわかります。

電磁場テンソル

ここでは電磁場がゲージポテンシャルから得られることやベクトル解析の知識は既知としています。詳しく知りたい方は、例えば次のような電磁気学の教科書を参照してください。

ひとまず直交座標系 $x^a=(t,x,y,z)$ で考えることにします。電磁場 $(\vec{E},\vec{B})$ はゲージポテンシャル $A_a=(-\phi,\vec{A})$ （ただし $\vec{A}=(A_x,A_y,A_z)$ ）から次のようにして計算されます：

$\qquad\displaystyle{F_{ab}=\frac{\partial A_b}{\partial x^a}-\frac{\partial A_a}{\partial x^b}}$

$F_{ab}$ を行列のようにして表示すると

$\qquad\displaystyle{F_{ab}=\left(\begin{array}{cccc} 0 & \frac{\partial A_x}{\partial t}+\frac{\partial \phi}{\partial x} & \frac{\partial A_y}{\partial t}+\frac{\partial \phi}{\partial y} & \frac{\partial A_z}{\partial t}-\frac{\partial \phi}{\partial z} \\ -\frac{\partial \phi}{\partial x}-\frac{\partial A_x}{\partial t} & 0 & \frac{\partial A_y}{\partial x}-\frac{\partial A_x}{\partial y} & \frac{\partial A_z}{\partial x}-\frac{\partial A_x}{\partial z} \\ -\frac{\partial \phi}{\partial y}-\frac{\partial A_y}{\partial t} & \frac{\partial A_x}{\partial y}-\frac{\partial A_y}{\partial x} & 0 & \frac{\partial A_z}{\partial y}-\frac{\partial A_y}{\partial z} \\ -\frac{\partial \phi}{\partial z}-\frac{\partial A_z}{\partial t} & \frac{\partial A_x}{\partial z}-\frac{\partial A_z}{\partial x} & \frac{\partial A_y}{\partial z}-\frac{\partial A_z}{\partial y} &0 \end{array} \right)}$

であり、電場と磁場が

$\qquad\displaystyle{\vec{E}= -\vec{\nabla}\,\phi -\frac{\partial \vec{A}}{\partial t} \qquad\quad \vec{B}= \vec{\nabla}\times\vec{A}}$

で与えられることから

$\qquad\displaystyle{F_{ab}=\left(\begin{array}{cccc} 0 & -E_x & -E_y & -E_z \\ E_x & 0 & B_z & -B_y \\ E_y & -B_z & 0 & B_x \\ E_z & B_y & -B_x &0 \end{array} \right)}$

であることがわかります。実はこの量 $F_{ab}$ はテンソル場で、曲がった座標系 $\xi^i$ に移るには変換則

$\qquad\displaystyle{F_{ij}=\frac{\partial x^a}{\partial \xi^i}\frac{\partial x^b}{\partial \xi^j}F_{ab}}$

を適用すればよいことが知られています。極座標系 $\xi^i=(t,r,\theta,\phi)$ をとった場合、具体的には

$\qquad\displaystyle{t=t\\ x=r\sin\theta\cos\phi \\ y=r\sin\theta\sin\phi \\ z=r\cos\theta}$

という関係から、この座標系での電磁場テンソルを求めることができます。例えば第２式の両辺を $r$ で偏微分して

$\qquad\displaystyle{\frac{\partial x}{\partial r}=\sin\theta\cos\phi=\frac{x}{r}}$

です。同様にして

$\qquad\displaystyle{\frac{\partial y}{\partial r}=\frac{y}{r} \qquad\quad \frac{\partial z}{\partial r}=\frac{z}{r}}$

がわかります。さて、原点付近に局在する電荷 $Q$ が遠方に作る電磁場は

$\qquad\displaystyle{\vec{E}=-\vec{\nabla}\frac{Q}{r}\qquad\quad \vec{B}=0}$

で与えられます（クーロンの法則）。電場の各成分は

$\qquad\displaystyle{E_x=\frac{Q}{r^2}\frac{x}{r}\qquad E_y=\frac{Q}{r^2}\frac{y}{r}\qquad E_z=\frac{Q}{r^2}\frac{z}{r}}$

です。上で述べた変換則から、極座標系において

$\qquad\displaystyle{F_{tr}=\frac{\partial x^a}{\partial t}\frac{\partial x^b}{\partial r}F_{ab}=\frac{\partial x^b}{\partial r}F_{0b}=-\frac{\partial x}{\partial r}E_x-\frac{\partial y}{\partial r}E_y-\frac{\partial z}{\partial r}E_z=-\frac{Q}{r^2} }$

がわかります。少し紛らわしいですが、極座標系において電磁場テンソルはこのとき

$\qquad\displaystyle{F_{ij}=\left(\begin{array}{cccc} 0 & -\frac{Q}{r^2} & 0 & 0 \\ \frac{Q}{r^2} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 &0 \end{array} \right)}$

となります。
　

電磁場のエネルギー・運動量テンソル

上で導入した電磁場テンソル $F_{ij}$ を使うと、エネルギー・運動量テンソル $T_{ij}$ が

$\qquad\displaystyle{T_{ij}=\frac{1}{4} g_{ij} F_{kl} F^{kl}-F_{ik} F_{jl} g^{kl}}$

で与えられます。ただし

$\qquad\displaystyle{F^{ij}=g^{ik} g^{jl}F_{kl} }$

です。これがどのような量であるかを見るために、平坦な時空の直交座標系での表式を求めてみましょう。このとき電磁場テンソルは

$\qquad\displaystyle{F_{ab}=\left(\begin{array}{cccc} 0 & -E_x & -E_y & -E_z \\ E_x & 0 & B_z & -B_y \\ E_y & -B_z & 0 & B_x \\ E_z & B_y & -B_x &0 \end{array} \right)}$

また、時空の計量は

$\qquad\displaystyle{\eta_{ab}=\eta^{ab}=\left(\begin{array}{cccc} -1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 &1 \end{array} \right)}$

です。以上の情報から

%reset
import sympy as sp

sp.var('Ex, Ey, Ez, Bx, By, Bz')

F = sp.Matrix([
[0, -Ex, -Ey, -Ez],
[Ex, 0, Bz, -By],
[Ey, -Bz, 0, Bx],
[Ez, By, -Bx, 0] ])

eta = sp.Matrix([
[-1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1] ])

einv = eta.inv()

Finv = einv*F*einv

T = -sp.Rational(1,4)*((F*Finv).trace())*eta + F*einv*F

T

より、エネルギー・運動量テンソル $T_{ab}$ の各成分が以下のように求まります。

$\qquad\displaystyle{T_{ab}=\left(\begin{array}{cccc} -u & S_x & S_y & S_z \\ S_x & M_{xx} & M_{xy} & M_{xz}\\ S_y & M_{yx} & M_{yy} & M_{yz} \\ S_z & M_{zx} & M_{zy} & M_{zz} \end{array}\right)}$

ただし

$\qquad\displaystyle{u=\frac{1}{2}(|\vec{E}|^2+|\vec{B}|^2)\\ \vec{S}=\vec{E}\times\vec{B} \\ M_{ij}=E_iE_j+B_iB_j-\frac{1}{2}\delta_{ij}(|\vec{E}|^2+|\vec{B}|^2) \qquad(i,j\in\{x,y,z\}) }$

で、それぞれエネルギー密度、ポインティング・ベクトル、マクスウェルの応力テンソル
として知られている量です。
　

ライスナー=ノルドシュトロム解

ではいよいよ本題のライスナー=ノルドシュトロム解（以下RN解）に取り組みます。RN解は、電場に由来するエネルギー・運動量テンソル $T_{ij}$ が存在する場合のアインシュタイン方程式

$\qquad\displaystyle{G_{ij}+\Lambda\,g_{ij}=8\pi \mathcal{T}_{ij} \qquad\quad \left(\mathcal{T}_{ij}=-\frac{1}{4\pi}T_{ij}\right)}$

の解です。上で見たように、極座標系で電磁場テンソルが

$\qquad\displaystyle{F_{ij}=\left(\begin{array}{cccc} 0 & -\frac{Q}{r^2} & 0 & 0 \\ \frac{Q}{r^2} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 &0 \end{array} \right)}$

で与えられる場合について考えます。これまでの問題設定（Schwarzschild解、de Sitter-Schwarzschild解）と同様に、静的で球対称な場合であることに変わりはないので、計量の形として

$\qquad\displaystyle{g_{ij}=\left(\begin{array}{cccc} -f(r) & 0 & 0 & 0 \\ 0 & \frac{1}{f(r)} & 0 & 0 \\ 0 & 0 & r^2 & 0 \\ 0 & 0 & 0 & r^2\sin^2\theta \end{array} \right)}$

を仮定します。ただし、これらの成分は座標系

$\qquad\displaystyle{\xi^{i}=\left(\begin{array}{c} t \\ r \\ \theta \\ \phi \end{array} \right)}$

に関するものです。このとき、アインシュタイン方程式のエネルギー・運動量テンソルの項を左辺に移項して

$\qquad\displaystyle{E_{ij}=G_{ij}+\Lambda\,g_{ij}+2T_{ij}}$

を計算してみましょう。

%reset
import sympy as sp

sp.var('t, r, theta, phi, psi, Lambda, Q')
xi = sp.Matrix([ [t],[r],[theta],[phi] ])
D = len(xi)
f = sp.Function('f')(r)

g = sp.Matrix([
[-f, 0, 0, 0],
[0, 1/f, 0, 0],
[0, 0, r**2, 0],
[0, 0, 0, r**2*sp.sin(theta)**2] ])

ginv = g.inv()

def Gamma(k,i,j):
    gm = 0
    for n in range(D):
        gm = gm + sp.Rational(1,2)*ginv[k,n]*(g[n,i].diff(xi[j])\
                +g[n,j].diff(xi[i])-g[i,j].diff(xi[n]))
    return gm

def R(i,j,k,l):
    rm = Gamma(i,l,j).diff(xi[k]) - Gamma(i,k,j).diff(xi[l]) 
    for n in range(D):
        rm = rm + Gamma(n,l,j)*Gamma(i,k,n) - Gamma(n,k,j)*Gamma(i,l,n)
    return rm

def Ricc(j,k):
    ric = 0
    for n in range(D):
        ric = ric + R(n,j,n,k)
    return ric

RiccMat = sp.Matrix([
[Ricc(0,0), Ricc(0,1), Ricc(0,2), Ricc(0,3)],
[Ricc(1,0), Ricc(1,1), Ricc(1,2), Ricc(1,3)],
[Ricc(2,0), Ricc(2,1), Ricc(2,2), Ricc(2,3)],
[Ricc(3,0), Ricc(3,1), Ricc(3,2), Ricc(3,3)] ])

R = sp.simplify((ginv*RiccMat).trace())

GMat = sp.simplify(RiccMat - sp.Rational(1,2)*g*R)

# 電磁場テンソル
F = sp.Matrix([
[0, -Q/r**2, 0, 0],
[Q/r**2, 0, 0, 0],
[0, 0, 0, 0],
[0, 0, 0, 0] ])

Finv = ginv*F*ginv

# エネルギー・運動量テンソル
T = -sp.Rational(1,4)*((F*Finv).trace())*g + F*ginv*F

EMat = sp.simplify(GMat + Lambda*g + 2*T)
EMat

より、未知関数 $f(r)$ について、以下が満たされれば解になることがわかります。

$\qquad\displaystyle{0=\Lambda r^{4} + Q^{2} + r^{2} \left(r \frac{d}{d r} f{\left(r \right)} + f{\left(r \right)} - 1\right)\\ 0=\Lambda r^{4}-Q^{2} + \frac{r^{3}}{2} \left(r \frac{d^{2}}{d r^{2}} f{\left(r \right)} + 2 \frac{d}{d r} f{\left(r \right)}\right)}$

ここで第２式は第１式を仮定すると自然に満たされるので、第１式だけ考えればよいことになります。よって１階の常微分方程式

$\qquad\displaystyle{f'=\frac{1-f}{r}-\Lambda r-\frac{Q^2}{r^3}}$

を得ます。微分方程式の解の関数形として

$\qquad\displaystyle{f(r)=1-\frac{a}{r}-br^2-\frac{c}{r^2}}$

を考えてみましょう。

%reset
import sympy as sp

sp.var('r, Lambda, Q, a, b, c')

f = 1 - a/r - b*r**2 - c/(r**2)

sp.simplify( f.diff(r)-(1-f)/r+Lambda*r+Q**2/(r**3) )

より

$\qquad\displaystyle{b=\frac{\Lambda}{3}\qquad\quad c=-Q^2 }$

がわかります。よって $a$ を定数として

$\qquad\displaystyle{f(r)=1-\frac{a}{r}-\frac{\Lambda}{3}r^2+\frac{Q^2}{r^2}}$

によって定まる計量はアインシュタイン方程式の解であり、これがRN解です。
　　

キーワード：Python、SymPy、Riemann幾何、一般相対論、RN解