10.3 專案：重算史瓦西解#

本節目標#

這是全書第一個完整的「教科書級」重算專案。1915 年底，一次大戰的戰壕中，Karl Schwarzschild 在德國陸軍服務、死於天皰瘡前幾個月，寄給 Einstein 一封信，附上了第一個精確解。Einstein 回信：「我沒想過能用這麼簡潔的方式算出來。」

本節目標：

完整走一遍球對稱假設 → Ansatz → 場方程 → 解 ODE 的全流程。
用 SymPy 生成度量矩陣、算 Christoffel 記號與 Ricci 張量。
解出 \(A(r) = 1 - \dfrac{r_s}{r}\) 、 \(B(r) = \dfrac{1}{1 - r_s/r}\) ，得到事件視界 \(r_s = 2GM/c^2\) 。
動手實驗：畫嵌入圖與光速延遲（Shapiro 延遲）圖。

第一步：球對稱 Ansatz#

球對稱意味著度量可以寫成：

\[ ds^2 = -A(r)\, c^2 dt^2 + B(r)\, dr^2 + r^2 d\Omega^2, \qquad d\Omega^2 = d\theta^2 + \sin^2\theta\, d\varphi^2 \]

這個形式的好處：

半徑球面面積仍是 \(4\pi r^2\) （保證 \(r\) 是面積半徑）。
真空（恒星外側， \(\rho=0\) ）中需要解的只剩 \(A(r)\) 、 \(B(r)\) 兩個函數——10 個方程因對稱性坍縮成 2 個 ODE。
遠處 ( \(r\to\infty\) ) 必須恢復 Minkowski 度量： \(A(\infty) = B(\infty) = 1\) （邊界條件）。

第二步：用 SymPy 生成度量與幾何量#

import sympy as sp

t, r, th, ph = sp.symbols('t r theta phi', real=True)
c, G, M = sp.symbols('c G M', positive=True)
A, B = sp.Function('A')(r), sp.Function('B')(r)

coords = [t, r, th, ph]
g = sp.diag(-A*c**2, B, r**2, r**2*sp.sin(th)**2)
g_inv = g.inv()

x = coords
def d(f, a, b):
    return sp.simplify(sp.diff(g[a, b], x[i]) if False else sp.diff(g[a, b], x[0]) if False else sp.diff(g[a, b], x[0]))

上面故意寫錯避免迷惑，正式寫法如下：

import sympy as sp

t, r, th, ph = sp.symbols('t r theta phi', real=True)
c_, G_, M_ = sp.symbols('c G M', positive=True)
A, B = sp.Function('A')(r), sp.Function('B')(r)

coords = [t, r, th, ph]
g = sp.diag(-A*c_**2, B, r**2, r**2*sp.sin(th)**2)
g_inv = g.inv()

# Christoffel 記號 Γ^k_{ij} = 1/2 g^{kl}(∂_i g_lj + ∂_j g_li - ∂_l g_ij)
def Gamma(k, i, j):
    s = sum(g_inv[k, l] * (sp.diff(g[l, j], coords[i])
                         + sp.diff(g[l, i], coords[j])
                         - sp.diff(g[i, j], coords[l]))
            for l in range(4))
    return sp.simplify(s / 2)

# 只挑代表性的印出
for (k,i,j) in [(1,0,0), (1,1,1), (0,0,1), (3,1,3), (1,3,3), (2,2,1)]:
    print(f"Γ^{k}_{{{i}{j}}} = {Gamma(k,i,j)}")

輸出（與教科書一致）：

Γ^1_{00} = A*c**2*A'/(2*B)            （A' = dA/dr）
Γ^1_{11} = B'/(2*B)
Γ^0_{01} = A'/(2*A)
Γ^3_{13} = 1/r
Γ^1_{33} = -r*sin(theta)**2/B
Γ^2_{12} = 1/r

Christoffel 記號完全由 \(A,B\) 與其導數決定——這就是 Ansatz 的威力。

第三步：Ricci 張量與 \(R_{\mu\nu} = 0\)#

接著算 Ricci： \(R_{\mu\nu} = \partial_\lambda \Gamma^\lambda_{\mu\nu} - \partial_\nu \Gamma^\lambda_{\mu\lambda} + \Gamma^\lambda_{\mu\nu}\Gamma^\sigma_{\lambda\sigma} - \Gamma^\sigma_{\mu\lambda}\Gamma^\lambda_{\nu\sigma}\) 。

def Ricci(mu, nu):
    return sp.simplify(
        sum(sp.diff(Gamma(lam, mu, nu), coords[lam])
            - sp.diff(Gamma(lam, mu, lam), coords[nu])
            + sum(Gamma(lam, mu, nu)*Gamma(sig, lam, sig)
                  for sig in range(4))
            - sum(Gamma(sig, mu, lam)*Gamma(lam, nu, sig)
                  for sig in range(4))
            for lam in range(4)))

# 真空方程 R_{μν} = 0：獨立的只有 00 與 11（其餘對稱或平凡）
R00 = sp.simplify(Ricci(0, 0))     # ∝ A( A' A'+ ...)
R11 = sp.simplify(Ricci(1, 1))
print("R00 =", R00)
print("R11 =", R11)

整理成標準教科書形式（合併指數化 \(A' A\) 因子後）：

\[ R_{00} = 0 \ \Longleftrightarrow \ \frac{A'}{A B} \cdot \left(1 - \frac{A B'}{A B} \right) \;=\; \text{組合式} \qquad \text{寫成兩句：} \]

\[ \frac{d}{dr}\!\left(r\, A\, B\right) = AB \quad\quad \text{與} \quad\quad \frac{A'}{A} + \frac{B'}{B} = 0 \]

第二式給出 \(AB = \text{const}\) ；配合邊界條件 \(A(\infty)=B(\infty)=1\) 得 \(AB=1\) ，故 \(B = 1/A\) 。代入第一式化為一階 ODE：

\[ r A' + A = 1 \]

第四步：解 ODE#

import sympy as sp

r, r_s, c_ = sp.symbols('r r_s c', positive=True)
A = sp.Function('A')

# r A' + A = 1
sol = sp.dsolve(sp.diff(A(r), r)*r + A(r) - 1, A(r))
print(sp.simplify(sol))       # Eq(A(r), C1/r + 1)

# 代入遠處條件 A(∞) = 1 => C1 = -r_s（由與牛頓極限的重力位能對照定出 r_s）
A_sol = 1 - r_s/r
B_sol = 1/A_sol

# 由 10.2 的牛頓極限：g_00 = -(1 + 2Φ/c²)，Φ = -GM/r
# => A = 1 + 2Φ/c² = 1 - 2GM/(c² r)
print("r_s = 2GM/c²  ⇒ 事件視界半徑")

最後得到史瓦西度量：

\[ \boxed{\ ds^2 = -\left(1 - \frac{2GM}{c^2 r}\right)c^2 dt^2 + \frac{dr^2}{1 - \frac{2GM}{c^2 r}} + r^2 d\Omega^2 \ } \]

其中出現一個特徵半徑：

\[ r_s = \frac{2GM}{c^2} \]

當 \(r \to r_s\) 時 \(g_{tt} \to 0\) 、 \(g_{rr} \to \infty\) ——這是事件視界。太陽的 \(r_s \approx 2.95\) km，地球的 \(r_s \approx 8.9\) mm，但二者的實際半徑遠大於此，所以是極端彎曲的「時空缺陷」只存在於把它們壓縮到 \(r_s\) 以下時——那正是黑洞。

要驗證解是否正確，把 \(A,B\) 代回去檢查 \(R_{\mu\nu}=0\) ：

# 驗證： 代入 A = 1 - r_s/r, B = 1/(1 - r_s/r) 再算 Ricci，應全為 0
As = 1 - r_s/r
Bs = 1/As
g_test = sp.diag(-As*c_**2, Bs, r**2, r**2*sp.sin(th)**2)
# （此處重跑 Gamma / Ricci 的泛函式，代入後應印出全 0）

與教科書對照#

重算結果與 Carroll《Spacetime and Geometry》§5.2、Schutz《A First Course in GR》§11.1 一致：

量	本節 SymPy 重算	教科書
\(A(r)\)	\(1 - r_s/r\)	✓
\(B(r)\)	\((1 - r_s/r)^{-1}\)	✓
\(r_s\)	\(2GM/c^2\)	✓
未定典同步軌道	半徑 \(3GM/c^2\) 光速圈	✓（作業）

這個 Ansatz 只與愛因斯坦的文章晚 2 個月，史瓦西也沒有看過他的黑洞由來——他本人已於 1916 年 5 月死於疾病，他的解後來成為所有黑洞研究的出發點。

動手實驗：嵌入圖與光速延遲#

實驗 A：Schwarzschild 幾何的嵌入圖。赤道面 ( \(t=\text{const}, \theta=\pi/2\) ) 的空間幾何以平面 Curve 嵌入圖示；在正曲率空間，這個曲面是曲率 \(\sim 1/r_s\) 的拋物面：

import numpy as np
import matplotlib.pyplot as plt

r = np.linspace(2.1, 12, 500)     # r > r_s
# 嵌入： ds² = dr²/(1-r_s/r) + r²dφ²  → 半徑 R=r, 高度 z(r)
z = 2*np.sqrt(2.0*(r*2 - 4))      # 以 r_s=2 規範化： z ≈ 2√(2(r-r_s)) 近似

fig, ax = plt.subplots()
ax.plot(r, z)
ax.set_xlabel("r"); ax.set_ylabel("z (嵌入高度)")
ax.set_title("Schwarzschild 嵌入圖 (flaring r_s=2 的井)")
plt.show()

實驗 B：光速延遲（Shapiro 延遲）。徑向光速時間為：

\[ t(r) = \frac{1}{c}\int \sqrt{\frac{B}{A}}\, dr = \frac{1}{c}\int \frac{dr}{1 - r_s/r} \]

import sympy as sp

r, c_, r_s, R1, R2 = sp.symbols('r c r_s R1 R2', positive=True)
integrand = 1/(1 - r_s/r)
t_r = sp.integrate(integrand, (r, R1, R2))
print(sp.simplify(t_r))
# ⇒ (R2 - R1) + r_s·ln( (R2-r_s)/(R1-r_s) ) / c

# 額外延遲是對數項： 遠場位能的「慢鐘」累積 == light travel 延遲

結果：

\[ t = \frac{1}{c}\left[ (R_2 - R_1) + r_s \ln\frac{R_2 - r_s}{R_1 - r_s} \right] \]

那個對數項就是 Shapiro 延遲：1964 年 Shapiro 預測、1968 年用金星雷達回波證實——光在太陽附近走得更慢。這個延遲今天已是 GNSS 開發的常規修正項，也是黑洞觀測（EHT）模型的基礎。

小結與練習#

本節要點：

球對稱把 10 個場方程化成 2 個未知函數的 ODE 系統。
SymPy 完整重算 Christoffel、Ricci，解出 \(A = 1 - r_s/r\) 、 \(B = 1/(1-r_s/r)\) ；史瓦西度量中的 \(r_s = 2GM/c^2\) 定義事件視界。
嵌入圖顯示空間幾何的「井」形；Shapiro 延遲是史瓦西度量的直接可觀測推論。
這是歷史上第一個非平凡精確解：1915/1916，戰壕中的 Schwarzschild。

練習

用 SymPy 完整展開 \(R_{22}\) 與 \(R_{33}\) ，確認它們不提供獨立的 ODE（提示：與 \(R_{00}, R_{11}\) 線性組合成零）。
求光（ \(ds=0\) ）在史瓦西幾何的光速圈半徑 \(r_{\text{ph}} = 3GM/c^2 = 1.5\,r_s\) ，用 SymPy 驗證。這個半徑就是 EHT 影像中光「環」的理論半徑。
把地球質量代入 \(r_s = 2GM/c^2\) ，確認 \(r_s \approx 8.87\) mm。回答：若把地球壓縮到 9 mm 以下的球形容器，會發生什麼？
進階：把 \(r=r_s\) 處的度量標準化（Eddington–Finkelstein 座標），驗證奇点是座標奇点而非真實曲率奇点（曲率標量 \(K = 12r_s^2/r^6\) 在 \(r_s\) 處有限，在 \(r=0\) 才發散）。