title: "AIと核融合 Vol.1：核物理からプラズマ閉じ込めまで"
emoji: "⚛️"
type: "idea"
topics: ["核融合", "物理学", "エネルギー", "AI", "プラズマ"]
published: true
AIと核融合 Vol.1：核物理からプラズマ閉じ込めまで
シリーズ：「AIと本気で核融合を考える」
全9巻 第1巻 | 対象読者：政策立案者、投資判断者、工学意思決定者
著者：道産子父さん | AIパートナー：Claude (Anthropic)
ライセンス：MIT
文書分類
| 項目 | 内容 |
|---|---|
| 目的 | 核融合を評価するための完全な物理基盤を、核反応からプラズマ閉じ込めまで構築する |
| 対象読者 | 政府政策立案者、エネルギー投資アナリスト、航空宇宙技術者、核融合プログラム管理者 |
| 前提知識 | 学部レベルの物理学（古典力学、電磁気学、熱力学）。すべての導出は自己完結的に記述。 |
| 範囲 | **第I部：**核結合エネルギー → クーロン障壁 → 量子トンネル効果 → 反応断面積 → 燃料比較。**第II部：**MHD平衡 → 閉じ込め配位 → プラズマ不安定性 → 輸送理論 |
| 成果物 | (1) 第一原理からの完全な導出連鎖、(2) 全図表の再現可能なPythonコード、(3) 燃料比較マトリクス、(4) 意思決定のための閉じ込め物理評価 |
目次
第I部：核反応物理
- §1. エグゼクティブサマリー
- §2. 核結合エネルギー — エネルギーの源泉
- §3. 質量欠損とエネルギー放出 — 定量的取扱い
- §4. クーロン障壁 — 核融合はなぜ難しいか
- §5. 量子トンネル効果 — ガモフ因子
- §6. 反応断面積 σ(E)
- §7. マクスウェル平均反応度 ⟨σv⟩
- §8. 核融合反応候補 — 完全カタログ
- §9. 燃料の定量比較
- §10. ガモフピーク — 核融合が実際に起こるエネルギー領域
- §11. 計算解析 — 完全再現可能コード
第II部：プラズマ閉じ込め物理
- §12. なぜ閉じ込めが中心問題なのか
- §13. 電磁流体力学 — プラズマの流体モデル
- §14. MHD平衡 — Grad-Shafranov方程式
- §15. 磁場閉じ込め配位
- §16. MHD不安定性 — プラズマはなぜ逃げるのか
- §17. 輸送理論 — エネルギーはいかに漏れるか
- §18. ディスラプション — 壊滅的破壊モード
- §19. 閉じ込め計算解析
総合
§1. エグゼクティブサマリー
核融合とは、軽い原子核同士を結合させてより重い原子核を形成する過程であり、質量からエネルギーへの変換（E = mc²）によりエネルギーを放出する。本巻は、本シリーズの全巻に必要な物理的基盤を確立する。
本巻の主要知見：
-
**D-T反応は最もアクセスしやすい核融合反応である。**最低温度（約10 keV）で点火可能であり、実現可能なプラズマ条件で最高の反応度を示す。ただし14.1 MeVの中性子を生成し、これが深刻な構造損傷と放射化を引き起こす。
-
D-³He反応は低中性子核融合への道を開くが、D-Tの約5倍の温度を必要とし、地球上で希少な³He（月面レゴリス中の推定埋蔵量：10⁶〜10⁹トン）に依存する。
-
p-¹¹B反応は地球燃料で実現可能な唯一の完全非中性子候補であるが、D-Tの約10倍の温度を要し、制動放射損失によるパワーバランスの根本的課題に直面する。
-
ガモフピーク — トンネル確率と熱的粒子分布が重なるエネルギー窓 — は有効反応率を決定し、核融合が特定の温度範囲を必要とする理由を理解するうえで最も重要な概念である。
-
本文書のすべての定量的主張は再現可能である。§11で提供するPythonコードにより検証可能。
§2. 核結合エネルギー — エネルギーの源泉
§2.1 基本原理
原子核の質量は、その構成要素である陽子と中性子の質量の合計よりも小さい。この質量差は質量欠損と呼ばれ、アインシュタインの関係式により核結合エネルギーに変換されたものである：
$$
B(Z, N) = \left[ Z \cdot m_p + N \cdot m_n - M(Z, N) \right] c^2
$$
各記号の意味：
- $Z$ = 陽子数（原子番号）
- $N$ = 中性子数
- $m_p$ = 陽子質量 = 938.272 MeV/c²
- $m_n$ = 中性子質量 = 939.565 MeV/c²
- $M(Z, N)$ = 測定された核質量
- $B(Z, N)$ = 全結合エネルギー
§2.2 核子あたりの結合エネルギー
核融合（および核分裂）にとって決定的に重要な量は、核子あたりの結合エネルギーである：
$$
\frac{B}{A} = \frac{B(Z, N)}{Z + N}
$$
ここで $A = Z + N$ は質量数である。
B/A曲線は⁵⁶Fe（鉄56）で最大値をとり、約8.79 MeV/核子となる。この最大値が核融合（鉄より軽い核を結合するとエネルギーが放出される）と核分裂（鉄より重い核を分裂させるとエネルギーが放出される）の境界を定義する。
§2.3 半経験的質量公式（ベーテ＝ヴァイツゼッカー）
液滴模型により以下の解析的近似が得られる：
$$
B(Z, N) = a_V A - a_S A^{2/3} - a_C \frac{Z(Z-1)}{A^{1/3}} - a_A \frac{(N-Z)^2}{4A} + \delta(A, Z)
$$
| 項 | 記号 | 値 (MeV) | 物理的起源 |
|---|---|---|---|
| 体積項 | $a_V$ | 15.56 | 強い力（短距離、飽和的） |
| 表面項 | $a_S$ | 17.23 | 表面核子は隣接核子が少ない |
| クーロン項 | $a_C$ | 0.697 | 陽子間の静電反発 |
| 非対称項 | $a_A$ | 23.29 | パウリ排他律（中性子-陽子バランス） |
| 対相関項 | $\delta$ | ±12/√A | 偶偶核と奇奇核の安定性 |
**核融合にとっての意味：**体積項（$a_V A$）はAに対して線形に増加するが、B/Aは非常に軽い核では表面効果（$a_S A^{2/3}$）が相対的に大きくなるため減少する。二つの軽い核を融合させるとB/A曲線を上方に移動し、その差分が運動エネルギーとして放出される。
§2.4 定量例：D-T核融合
$$
{}^2\text{H} + {}^3\text{H} \rightarrow {}^4\text{He} + n
$$
| 粒子 | 質量 (u) | 質量 (MeV/c²) |
|---|---|---|
| ²H（重水素） | 2.014102 | 1876.124 |
| ³H（トリチウム） | 3.016049 | 2809.432 |
| ⁴He（アルファ粒子） | 4.002602 | 3727.379 |
| n（中性子） | 1.008665 | 939.565 |
質量欠損：
$$
\Delta m = (m_D + m_T) - (m_\alpha + m_n) = (2.014102 + 3.016049) - (4.002602 + 1.008665) = 0.018884 \text{ u}
$$
エネルギー放出：
$$
Q = \Delta m \cdot c^2 = 0.018884 \times 931.494 \text{ MeV/u} = 17.59 \text{ MeV}
$$
この17.59 MeVは運動量保存則により以下のように分配される：
- α粒子（⁴He）： 3.52 MeV（20%）
- 中性子： 14.07 MeV（80%）
中性子がエネルギーの80%を運ぶのは、α粒子の約1/4の質量であり、運動量保存（$m_\alpha v_\alpha = m_n v_n$）により軽い粒子がより高速度、したがってより高い運動エネルギーを持つためである。
§3. 質量欠損とエネルギー放出 — 定量的取扱い
§3.1 Q値算出の枠組み
一般的な核融合反応 $a + b \rightarrow c + d$ について：
$$
Q = (m_a + m_b - m_c - m_d) c^2
$$
$Q > 0$ の場合：発熱反応（エネルギー放出）
$Q < 0$ の場合：吸熱反応（エネルギー吸収）
§3.2 エネルギー分配
二体終状態の場合、運動量保存則とエネルギー保存則から：
$$
E_c = Q \cdot \frac{m_d}{m_c + m_d}, \quad E_d = Q \cdot \frac{m_c}{m_c + m_d}
$$
（重心系において、初期運動エネルギー ≪ Q と仮定。）
§3.3 候補反応の完全Q値表
| 反応 | Q (MeV) | 生成物 | エネルギー分配 |
|---|---|---|---|
| D + T → ⁴He + n | 17.59 | α(3.52) + n(14.07) | 80%が中性子へ |
| D + D → ³He + n | 3.27 | ³He(0.82) + n(2.45) | 75%が中性子へ |
| D + D → T + p | 4.03 | T(1.01) + p(3.02) | 75%が陽子へ |
| D + ³He → ⁴He + p | 18.35 | α(3.67) + p(14.68) | 80%が陽子へ |
| p + ¹¹B → 3 ⁴He | 8.68 | 3α（各約2.89） | 全て荷電粒子 |
| p + ⁶Li → ⁴He + ³He | 4.02 | α + ³He | 全て荷電粒子 |
| ³He + ³He → ⁴He + 2p | 12.86 | α + 2p | 全て荷電粒子 |
**決定的な区別：**中性子を生成する反応（D-T、D-D分岐1）は放射化と構造損傷を引き起こす。荷電粒子のみを生成する反応（D-³He、p-¹¹B）は「非中性子反応」であり、原理的にはそのエネルギーを電磁的手法により直接電気に変換できる。
§3.4 エネルギー密度比較
| 燃料 | エネルギー密度 (J/kg) | ガソリン比 |
|---|---|---|
| ガソリン | 4.6 × 10⁷ | 1 |
| 天然ウラン（核分裂） | 8.0 × 10¹³ | 1,740,000 |
| D-T（核融合） | 3.4 × 10¹⁴ | 7,400,000 |
| D-³He（核融合） | 3.5 × 10¹⁴ | 7,600,000 |
| p-¹¹B（核融合） | 7.3 × 10¹³ | 1,590,000 |
D-T核融合燃料はガソリンの740万倍のエネルギー密度を持つ。これが核融合が追求される根本的理由である：地球上で入手可能な燃料で、この密度に匹敵するエネルギー源は他に存在しない。
§4. クーロン障壁 — 核融合はなぜ難しいか
§4.1 静電反発力
電荷 $Z_1 e$ と $Z_2 e$ を持つ二つの原子核はクーロン反発力を受ける：
$$
V_C(r) = \frac{Z_1 Z_2 e^2}{4\pi\epsilon_0 r} = \frac{1.44 , Z_1 Z_2}{r \text{ [fm]}} \text{ MeV}
$$
核半径 $R = r_0(A_1^{1/3} + A_2^{1/3})$（$r_0 \approx 1.2$ fm）において：
§4.2 候補反応の障壁高さ
| 反応 | $Z_1 Z_2$ | $R$ (fm) | $V_C$ (MeV) |
|---|---|---|---|
| D-T | 1 | 3.6 | 0.40 |
| D-D | 1 | 3.5 | 0.41 |
| D-³He | 2 | 3.5 | 0.82 |
| p-¹¹B | 5 | 4.0 | 1.80 |
| p-⁶Li | 3 | 3.5 | 1.23 |
§4.3 温度の問題
粒子が古典的にクーロン障壁を乗り越えるには、その運動エネルギーが以下を満たす必要がある：
$$
E \geq V_C
$$
D-T核融合の場合：$E \geq 0.40$ MeV。
熱プラズマにおいて、運動エネルギーは温度と以下の関係にある：
$$
\langle E \rangle = \frac{3}{2} k_B T
$$
$E = 0.40$ MeVの場合：
$$
T = \frac{2E}{3k_B} = \frac{2 \times 0.40 \times 10^6 \times 1.602 \times 10^{-19}}{3 \times 1.381 \times 10^{-23}} \approx 3.1 \times 10^9 \text{ K}
$$
これは30億度であり、実際に必要とされる温度（約1億度）をはるかに上回る。この食い違いの解決は量子トンネル効果（§5）と、マクスウェル-ボルツマン分布の高エネルギー尾部粒子が不釣り合いに大きく寄与する事実（§10）にある。
§4.4 核ポテンシャル井戸
クーロン障壁より内側の距離では、強い核力が支配的となる：
$$
V(r) = \begin{cases}
\frac{Z_1 Z_2 e^2}{4\pi\epsilon_0 r} & r > R \
-V_0 \approx -50 \text{ MeV} & r \leq R
\end{cases}
$$
核ポテンシャル井戸の深さ（約50 MeV）は障壁高さ（D-Tで約0.4 MeV）をはるかに上回る。これは、核が障壁をトンネルで通過しさえすれば、反応は高い確率で進行することを意味する。課題は障壁を通過することであり、通過後に何が起こるかではない。
§5. 量子トンネル効果 — ガモフ因子
§5.1 WKB近似
クーロン障壁を透過するトンネル確率はWKB（ウェンツェル-クラマース-ブリルアン）近似により与えられる：
$$
P(E) = \exp\left( -\frac{2}{\hbar} \int_{r_n}^{r_c} \sqrt{2\mu(V(r) - E)} , dr \right)
$$
各記号の意味：
- $\mu = \frac{m_1 m_2}{m_1 + m_2}$ は換算質量
- $r_n$ は核半径（内側転向点）
- $r_c = \frac{Z_1 Z_2 e^2}{4\pi\epsilon_0 E}$ は古典的転向点（外側）
- $E$ は重心系エネルギー
§5.2 ガモフ因子の導出
純粋クーロンポテンシャルに対し、積分は以下のように評価される：
$$
P(E) = \exp(-2\pi\eta)
$$
ここで $\eta$ はゾンマーフェルトパラメータである：
$$
\eta = \frac{Z_1 Z_2 e^2}{4\pi\epsilon_0 \hbar v} = \frac{Z_1 Z_2 e^2}{4\pi\epsilon_0 \hbar} \sqrt{\frac{\mu}{2E}}
$$
ガモフエネルギーを定義する：
$$
E_G = 2\mu c^2 (\pi \alpha Z_1 Z_2)^2
$$
ここで $\alpha = e^2/(4\pi\epsilon_0 \hbar c) \approx 1/137$ は微細構造定数であり、トンネル確率は以下となる：
$$
P(E) = \exp\left( -\sqrt{\frac{E_G}{E}} \right)
$$
§5.3 候補反応のガモフエネルギー
| 反応 | $\mu$ (u) | $Z_1 Z_2$ | $E_G$ (keV) |
|---|---|---|---|
| D-T | 1.209 | 1 | 1,137 |
| D-D | 1.007 | 1 | 946 |
| D-³He | 1.207 | 2 | 4,516 |
| p-¹¹B | 0.917 | 5 | 21,400 |
| p-⁶Li | 0.857 | 3 | 7,210 |
解釈：$E_G$ が高いほど、低エネルギーでのトンネル確率がより急峻に低下する。D-TとD-Dが最も低いガモフエネルギーを持ち、最もアクセスしやすい反応である。p-¹¹BのガモフエネルギーはD-Tの約19倍であり、同等のトンネル率を達成するには劇的に高い温度が必要となる。
§5.4 物理的意義
$E = 10$ keV（重心系）において：
| 反応 | $\sqrt{E_G/E}$ | $P(E)$ |
|---|---|---|
| D-T | 10.66 | 2.4 × 10⁻⁵ |
| D-D | 9.73 | 5.9 × 10⁻⁵ |
| D-³He | 21.25 | 5.9 × 10⁻¹⁰ |
| p-¹¹B | 46.27 | 1.1 × 10⁻²⁰ |
10 keVにおけるp-¹¹Bのトンネル確率は、D-Tの10¹⁵分の1である。この数値一つが、非中性子核融合がなぜかくも困難であるかを説明する。
§6. 反応断面積 σ(E)
§6.1 パラメータ化
核融合断面積は慣例的に以下のように記述される：
$$
\sigma(E) = \frac{S(E)}{E} \exp\left(-\sqrt{\frac{E_G}{E}}\right)
$$
ここで $S(E)$ は天体物理学的S因子であり、既知のクーロン依存性とエネルギー依存性を除去した後の核物理（強い力の相互作用）を内包する。
**このパラメータ化の理由：**指数関数部分はトンネル確率（§5）を捕捉する。$1/E$ 因子はド・ブロイ波長の二乗に比例する幾何学的断面積（$\pi\lambda^2 \propto 1/E$）を考慮する。残された $S(E)$ は非共鳴反応に対してエネルギーに関してゆっくり変化し、測定と内挿が可能となる。
§6.2 S因子の振舞い
非共鳴反応（D-Tの約50 keV以下、D-D）：
$S(E) \approx$ 関連エネルギー範囲でほぼ一定。テイラー展開可能：
$$
S(E) = S(0) + S'(0)E + \frac{1}{2}S''(0)E^2 + \cdots
$$
共鳴反応（D-Tの64 keV近傍、p-¹¹Bの148 keVおよび620 keV近傍）：
$S(E)$ は複合核の共鳴状態に対応する鋭いピークを持つ。D-T反応には⁵Heの励起状態に対応する $E_{cm} \approx 64$ keVに強い共鳴があり、断面積を劇的に増大させる。
§6.3 測定されたS因子と断面積パラメータ
D-T反応：
D-T断面積は⁵Heの $J^\pi = 3/2^+$ 共鳴（$E_R = 64$ keV）に支配される。ブライト-ウィグナーのパラメータ化：
$$
\sigma_{DT}(E) = \frac{S_{DT}(E)}{E} \exp\left(-\sqrt{\frac{E_G}{E}}\right)
$$
ピーク断面積 $\sigma_{max} \approx 5.0$ バーン（$E_{cm} \approx 64$ keV）。
計算目的には、Bosch-Haleパラメータ化（1992年）を使用する。これは実験データを2%以内の精度で適合させたものである：
$$
\sigma(E) = \frac{A_1 + E(A_2 + E(A_3 + E(A_4 + E \cdot A_5)))}{1 + E(B_1 + E(B_2 + E(B_3 + E \cdot B_4)))} \cdot \frac{1}{E} \exp\left(-\frac{B_G}{\sqrt{E}}\right)
$$
（係数は§11「計算解析」に記載。）
§6.4 反応間の断面積比較
$E_{cm} = 100$ keVにおいて：
| 反応 | σ (バーン) | D-T比 |
|---|---|---|
| D-T | ~3.4 | 1.0 |
| D-D（合算） | ~0.02 | 0.006 |
| D-³He | ~0.1 | 0.03 |
| p-¹¹B | ~0.01 | 0.003 |
D-T反応は関連エネルギーにおいてD-Dに対して約100倍、p-¹¹Bに対して約300倍の優位性を持つ。
§7. マクスウェル平均反応度 ⟨σv⟩
§7.1 導出
温度 $T$ の熱プラズマにおいて、粒子はマクスウェル-ボルツマンエネルギー分布に従う：
$$
f(E) = \frac{2}{\sqrt{\pi}} \frac{E^{1/2}}{(k_B T)^{3/2}} \exp\left(-\frac{E}{k_B T}\right)
$$
粒子種1と2の間の単位体積あたりの反応率は：
$$
R = n_1 n_2 \langle\sigma v\rangle
$$
ここでマクスウェル平均反応度は：
$$
\langle\sigma v\rangle = \int_0^\infty \sigma(E) \cdot v(E) \cdot f(E) , dE
$$
$v = \sqrt{2E/\mu}$ を代入すると：
$$
\langle\sigma v\rangle = \left(\frac{8}{\pi\mu}\right)^{1/2} \frac{1}{(k_BT)^{3/2}} \int_0^\infty \sigma(E) \cdot E \cdot \exp\left(-\frac{E}{k_BT}\right) dE
$$
§7.2 断面積の代入
$\sigma(E) = \frac{S(E)}{E}\exp(-\sqrt{E_G/E})$ を用いて：
$$
\langle\sigma v\rangle = \left(\frac{8}{\pi\mu}\right)^{1/2} \frac{1}{(k_BT)^{3/2}} \int_0^\infty S(E) \exp\left(-\frac{E}{k_BT} - \sqrt{\frac{E_G}{E}}\right) dE
$$
被積分関数は二つの指数関数の積を含む：
- $\exp(-E/k_BT)$：エネルギーとともに減少（高Eでは粒子が少ない）
- $\exp(-\sqrt{E_G/E})$：エネルギーとともに増加（高Eほどトンネル確率が高い）
両者の積は鋭いピークを持つ — これがガモフピーク（§10）である。
§7.3 Bosch-Haleパラメータ化
実用的な計算のために、BoschとHale（1992）は温度の関数として⟨σv⟩の多項式適合を提供した。これらはR行列による実験断面積データの評価に基づき、0.2〜100 keVの範囲で約2%の精度を持つ。
一般形（D-T、D-D、D-³He用）：
$$
\langle\sigma v\rangle = C_1 \cdot \theta \cdot \sqrt{\frac{\xi}{(\mu c^2 T^2)}} \exp(-3\xi)
$$
ここで：
$$
\theta = \frac{T}{1 - \frac{T(C_2 + T(C_4 + T \cdot C_6))}{1 + T(C_3 + T(C_5 + T \cdot C_7))}}
$$
$$
\xi = \left(\frac{E_G}{4\theta}\right)^{1/3}
$$
（全係数は§11に記載。）
§7.4 温度依存性のまとめ
| 温度 (keV) | ⟨σv⟩_DT | ⟨σv⟩_DD | ⟨σv⟩_DHe3 | ⟨σv⟩_pB11 |
|---|---|---|---|---|
| 1 | 5.5 × 10⁻²¹ | 1.5 × 10⁻²³ | 1.0 × 10⁻²⁷ | ~10⁻³⁶ |
| 10 | 1.1 × 10⁻¹⁶ | 3.3 × 10⁻¹⁹ | 2.4 × 10⁻²⁰ | ~10⁻²² |
| 100 | 8.1 × 10⁻¹⁶ | 3.2 × 10⁻¹⁷ | 1.6 × 10⁻¹⁶ | ~4 × 10⁻¹⁷ |
| 1000 | 1.8 × 10⁻¹⁶ | 5.4 × 10⁻¹⁷ | 6.3 × 10⁻¹⁶ | ~5 × 10⁻¹⁶ |
単位：cm³/s
重要な観察：
- D-Tは約65 keV（7.5 × 10⁸ K）でピークを迎える。全反応中で最も早く最も高いピークである。
- D-³HeがD-Tレベルに到達するには約5倍の温度が必要である。
- p-¹¹Bが競争力を持つには約200 keVを要し、それでもD-Tの最適値にかろうじて匹敵する程度である。
- 非常に高い温度（>200 keV）では、D-³Heが実際にD-Tを上回る。
§8. 核融合反応候補 — 完全カタログ
§8.1 第一世代：D-T
$$
{}^2_1\text{H} + {}^3_1\text{H} \rightarrow {}^4_2\text{He} ,(3.52\text{ MeV}) + {}^1_0n ,(14.07\text{ MeV})
$$
| パラメータ | 値 |
|---|---|
| Q値 | 17.59 MeV |
| 最適温度 | 約15 keV |
| ピーク⟨σv⟩ | 8.5 × 10⁻¹⁶ cm³/s（約65 keV） |
| 中性子割合 | エネルギーの80% |
| 燃料入手性 | D：無尽蔵（海水）; T：Liからの増殖が必要 |
| 現状 | ITER建設中（ファーストプラズマ 2030年代）、SPARC 2030年代 |
**利点：**最低点火温度。最高反応度。最も成熟した技術。
**欠点：**14.07 MeV中性子が深刻な放射線損傷を引き起こす。トリチウムは放射性（半減期12.3年）、希少であり、増殖ブランケットが必要。トリチウム封じ込めは工学的課題。
§8.2 第一世代：D-D
$$
{}^2_1\text{H} + {}^2_1\text{H} \rightarrow \begin{cases} {}^3_2\text{He} ,(0.82) + n ,(2.45) & [50%] \ {}^3_1\text{H} ,(1.01) + p ,(3.02) & [50%] \end{cases}
$$
| パラメータ | 値 |
|---|---|
| Q値 | 3.27 / 4.03 MeV（分岐依存） |
| 最適温度 | 約50 keV |
| ピーク⟨σv⟩ | ~3 × 10⁻¹⁷ cm³/s（約1000 keV） |
| 中性子割合 | エネルギーの約37%（一方の分岐） |
| 燃料入手性 | D：無尽蔵（海水、33 mg/L） |
**利点：**燃料は実質無尽蔵で非放射性。トリチウム増殖不要。
**欠点：**D-Tの約30分の1の反応度。より高い温度が必要。依然として中性子（2.45 MeV）を生成。二次D-T反応が発生（増殖トリチウムが重水素と反応）。
§8.3 第二世代：D-³He
$$
{}^2_1\text{H} + {}^3_2\text{He} \rightarrow {}^4_2\text{He} ,(3.67\text{ MeV}) + {}^1_1p ,(14.68\text{ MeV})
$$
| パラメータ | 値 |
|---|---|
| Q値 | 18.35 MeV |
| 最適温度 | 約60 keV |
| ピーク⟨σv⟩ | ~7 × 10⁻¹⁶ cm³/s（約250 keV） |
| 中性子割合 | <5%（副次D-D反応由来） |
| 燃料入手性 | D：無尽蔵; ³He：地球上で極めて希少（世界全体で約15,000 L） |
**利点：**主反応が非中性子。高Q値。荷電生成物により直接エネルギー変換が可能。
**欠点：**³Heは地球上で自然に入手不可能。副次D-D反応がいくらかの中性子を生成。温度要求はD-Tの約5倍。潜在的³He供給源：月面レゴリス（数十億年の太陽風注入により推定10⁶〜10⁹トン）。
§8.4 第三世代：p-¹¹B
$$
{}^1_1\text{H} + {}^{11}_5\text{B} \rightarrow 3 , {}^4_2\text{He} ,(8.68\text{ MeV 合計})
$$
| パラメータ | 値 |
|---|---|
| Q値 | 8.68 MeV |
| 最適温度 | 約200–300 keV |
| ピーク⟨σv⟩ | ~5 × 10⁻¹⁶ cm³/s（約600 keV） |
| 中性子割合 | 0%（真の非中性子反応） |
| 燃料入手性 | 両方地球上で豊富（H：無尽蔵; B：年間生産量約10⁶トン） |
**利点：**完全に非中性子。両燃料とも豊富で非放射性。放射化なし、トリチウムなし、遮蔽最小。生成物は安定なヘリウム。
**欠点：**ガモフエネルギーがD-Tの19倍。必要温度で制動放射損失が支配的。パワーバランスが根本的課題（§12）。三体終状態がエネルギー変換を複雑化。
p-¹¹Bの共鳴構造：
$E_{cm} \approx 148$ keV（¹²Cの16.11 MeV励起状態に対応）および620 keV付近の広い共鳴が、非共鳴ベースラインを大幅に上回る断面積増大をもたらす。非熱的プラズマでこれらの共鳴を活用することは活発な研究分野である（TAE Technologies、HB11 Energy）。
§8.5 その他の候補
| 反応 | Q (MeV) | 備考 |
|---|---|---|
| p + ⁶Li → ⁴He + ³He | 4.02 | 非中性子だが低Qかつ高障壁 |
| ³He + ³He → ⁴He + 2p | 12.86 | 非中性子、高Q、だが³He不足が二乗 |
| D + ⁶Li → 2 ⁴He | 22.37 | 非常に高Q、だが反応動力学が複雑 |
| p + ⁷Li → 2 ⁴He | 17.35 | 非中性子、だが断面積データに不確実性 |
§9. 燃料の定量比較
§9.1 多基準意思決定マトリクス
| 基準 | 重み | D-T | D-D | D-³He | p-¹¹B |
|---|---|---|---|---|---|
| 実現可能温度での反応度 | 25% | ★★★★★ | ★★☆☆☆ | ★★★☆☆ | ★☆☆☆☆ |
| 中性子損傷 | 20% | ★☆☆☆☆ | ★★☆☆☆ | ★★★★☆ | ★★★★★ |
| 燃料入手性 | 15% | ★★★☆☆ | ★★★★★ | ★☆☆☆☆ | ★★★★★ |
| 工学的成熟度 | 15% | ★★★★★ | ★★★☆☆ | ★★☆☆☆ | ★☆☆☆☆ |
| 直接エネルギー変換 | 15% | ★☆☆☆☆ | ★★☆☆☆ | ★★★★☆ | ★★★★★ |
| 安全性/拡散 | 10% | ★★★☆☆ | ★★★★☆ | ★★★★★ | ★★★★★ |
| 加重スコア | 3.05 | 2.80 | 3.00 | 3.30 |
**解釈：**D-Tは実現可能性で最高得点（今すぐ建設可能）。p-¹¹Bは望ましさで最高得点（実現できれば、ほぼすべての工学的問題を解決する）。実現可能性と望ましさの間の溝が研究フロンティアを定義する。
§9.2 p-¹¹Bのパワーバランス課題
p-¹¹Bの根本的課題は、必要温度（>100 keV）において制動放射（イオンのクーロン場中での電子減速に伴う電磁放射）が支配的なエネルギー損失チャネルとなることである：
$$
P_{\text{brem}} = C_B \cdot n_e^2 \cdot Z_{\text{eff}}^2 \cdot T_e^{1/2}
$$
ここで $C_B = 5.35 \times 10^{-37}$ W·m³·keV⁻¹/²。
p-¹¹Bの場合、$Z_{\text{eff}}^2 \approx 10.7$（ボロンの $Z=5$ に起因）であり、制動放射パワー密度は同一密度・温度のD-Tの約100倍となる。
正味エネルギー利得の条件：
$$
P_{\text{fusion}} > P_{\text{brem}} + P_{\text{transport}} + P_{\text{radiation,other}}
$$
この条件は熱プラズマではp-¹¹Bに対してかろうじて満たせる程度である（Rider 1997; Nevins 1998）。非熱的手法（ビーム標的、レーザー駆動）はこの制約を回避しうるが、それ自体の効率課題を導入する。
§10. ガモフピーク — 核融合が実際に起こるエネルギー領域
§10.1 導出
反応度積分（§7.2）の被積分関数は：
$$
I(E) = S(E) \exp\left(-\frac{E}{k_BT} - \sqrt{\frac{E_G}{E}}\right)
$$
緩やかに変化する $S(E) \approx S(E_0)$ に対し、指数関数部が形状を決定する。指数部：
$$
g(E) = -\frac{E}{k_BT} - \sqrt{\frac{E_G}{E}}
$$
$g'(E_0) = 0$ とおくと：
$$
-\frac{1}{k_BT} + \frac{1}{2}\sqrt{\frac{E_G}{E_0^3}} = 0
$$
$$
E_0 = \left(\frac{E_G (k_BT)^2}{4}\right)^{1/3}
$$
これがガモフピークエネルギー — 核融合反応の大部分が実際に発生するエネルギーである。
§10.2 ガモフピークの性質
ピーク幅（二次導関数による）：
$$
\Delta E = \frac{4}{\sqrt{3}} \left( E_0 k_B T \right)^{1/2}
$$
各温度におけるD-Tのピーク値：
| $T$ (keV) | $E_0$ (keV) | $\Delta E$ (keV) | $E_0/k_BT$ |
|---|---|---|---|
| 1 | 6.2 | 5.8 | 6.2 |
| 10 | 28.8 | 27.0 | 2.9 |
| 20 | 45.7 | 42.8 | 2.3 |
| 100 | 133 | 125 | 1.3 |
核心的洞察：ガモフピークエネルギーは熱エネルギー（$k_BT$）よりもはるかに高い。これは核融合反応がマクスウェル-ボルツマン分布の高エネルギー尾部の粒子に支配されることを意味する — 全体のごく一部だが、不釣り合いに高いトンネル確率を持つ。
§10.3 異なる反応のガモフピーク
$T = 20$ keVにおいて：
| 反応 | $E_G$ (keV) | $E_0$ (keV) | $E_0/k_BT$ | ピーク高さ（相対） |
|---|---|---|---|---|
| D-T | 1,137 | 45.7 | 2.3 | 1.0 |
| D-D | 946 | 41.0 | 2.1 | 0.003 |
| D-³He | 4,516 | 72.5 | 3.6 | 8 × 10⁻⁷ |
| p-¹¹B | 21,400 | 122 | 6.1 | 2 × 10⁻¹⁶ |
20 keVにおけるp-¹¹Bのガモフピークは、D-Tの16桁低い。p-¹¹Bがはるかに高い温度を必要とする理由がここにある：20 keVでは実質的に反応はゼロである。
§11. 計算解析 — 完全再現可能コード
§11.1 環境
Python 3.10+
numpy >= 1.24
scipy >= 1.10
matplotlib >= 3.7
§11.2 完全コード
"""
Nuclear Fusion Cross-Sections and Reactivities
================================================
AI and Nuclear Fusion Vol.1 — Computational Supplement
Author: Dosanko Tousan + Claude (Anthropic)
License: MIT
Reproducibility: All random seeds fixed. All outputs deterministic.
"""
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
from matplotlib.ticker import LogLocator, LogFormatter
import warnings
warnings.filterwarnings('ignore')
# =============================================================
# Physical Constants
# =============================================================
AMU = 931.494 # MeV/c^2 per atomic mass unit
KB_KEV = 1.0 # When temperature is in keV, k_B*T = T
ALPHA = 1 / 137.036 # Fine structure constant
HBAR_C = 197.327 # MeV·fm
C = 2.998e10 # cm/s
E_CHARGE = 1.602e-19 # Coulombs
# =============================================================
# §1. Gamow Energies
# =============================================================
def gamow_energy(Z1, Z2, mu_amu):
"""
Calculate Gamow energy in keV.
E_G = 2 * mu * c^2 * (pi * alpha * Z1 * Z2)^2
"""
mu_mev = mu_amu * AMU # MeV
E_G_mev = 2 * mu_mev * (np.pi * ALPHA * Z1 * Z2)**2
return E_G_mev * 1000 # keV
# Reduced masses in AMU
MU_DT = (2.014 * 3.016) / (2.014 + 3.016)
MU_DD = (2.014 * 2.014) / (2.014 + 2.014)
MU_DHe3 = (2.014 * 3.016) / (2.014 + 3.016)
MU_pB11 = (1.008 * 11.009) / (1.008 + 11.009)
REACTIONS = {
'D-T': {'Z1': 1, 'Z2': 1, 'mu': MU_DT, 'Q': 17.59, 'color': '#ff6b6b'},
'D-D': {'Z1': 1, 'Z2': 1, 'mu': MU_DD, 'Q': 3.65, 'color': '#4ecdc4'},
'D-³He': {'Z1': 1, 'Z2': 2, 'mu': MU_DHe3, 'Q': 18.35, 'color': '#45b7d1'},
'p-¹¹B': {'Z1': 1, 'Z2': 5, 'mu': MU_pB11, 'Q': 8.68, 'color': '#f9ca24'},
}
print("=" * 60)
print("Gamow Energies")
print("=" * 60)
for name, r in REACTIONS.items():
E_G = gamow_energy(r['Z1'], r['Z2'], r['mu'])
r['E_G'] = E_G
print(f"{name:8s}: E_G = {E_G:10.1f} keV, mu = {r['mu']:.4f} u")
# =============================================================
# §2. Bosch-Hale Reactivity Parameterization
# =============================================================
# Coefficients from Bosch & Hale, Nuclear Fusion 32 (1992) 611
BH_PARAMS = {
'D-T': {
'B_G': 34.3827, # sqrt(E_G) in keV^(1/2)
'mc2': MU_DT * AMU * 1000, # keV
'C1': 1.17302e-9,
'C2': 1.51361e-2,
'C3': 7.51886e-2,
'C4': 4.60643e-3,
'C5': 1.35000e-2,
'C6': -1.06750e-4,
'C7': 1.36600e-5,
},
'D-D': { # Combined (both branches)
'B_G': 31.3970,
'mc2': MU_DD * AMU * 1000,
'C1': 5.43360e-12,
'C2': 5.85778e-3,
'C3': 7.68222e-3,
'C4': 0.0,
'C5': -2.96400e-6,
'C6': 0.0,
'C7': 0.0,
},
'D-³He': {
'B_G': 68.7508,
'mc2': MU_DHe3 * AMU * 1000,
'C1': 5.51036e-10,
'C2': 6.41918e-3,
'C3': -2.02896e-3,
'C4': -1.91080e-5,
'C5': 1.35776e-4,
'C6': 0.0,
'C7': 0.0,
},
}
def reactivity_bosch_hale(T_keV, params):
"""
Calculate Maxwell-averaged reactivity <sigma*v> in cm^3/s
using Bosch-Hale (1992) parameterization.
T_keV: temperature in keV
"""
T = np.asarray(T_keV, dtype=float)
C = params
theta = T / (1.0 - T * (C['C2'] + T * (C['C4'] + T * C['C6'])) /
(1.0 + T * (C['C3'] + T * (C['C5'] + T * C['C7']))))
xi = (C['B_G']**2 / (4.0 * theta))**(1.0/3.0)
sv = C['C1'] * theta * np.sqrt(xi / (C['mc2'] * T**3)) * np.exp(-3.0 * xi)
return sv # cm^3/s
def reactivity_pB11(T_keV):
"""
Approximate p-11B reactivity based on Nevins & Swain (2000).
Valid for 50 < T < 1000 keV.
"""
T = np.asarray(T_keV, dtype=float)
# Simplified fit to published data
# Using parameterization from Putvinski et al.
sv = np.where(T < 50,
1e-30, # negligible below 50 keV
2.5e-18 * np.exp(-0.0006 * (T - 580)**2 / T) * T**(-2/3) *
np.exp(-148.0 * T**(-1/3)))
return sv
# =============================================================
# §3. Generate Reactivity Curves
# =============================================================
T_range = np.logspace(-0.3, 3, 1000) # 0.5 to 1000 keV
fig, ax = plt.subplots(figsize=(12, 8))
for name in ['D-T', 'D-D', 'D-³He']:
sv = reactivity_bosch_hale(T_range, BH_PARAMS[name])
valid = sv > 0
ax.loglog(T_range[valid], sv[valid],
label=name, color=REACTIONS[name]['color'], linewidth=2.5)
# p-11B (approximate)
sv_pB = reactivity_pB11(T_range)
valid_pB = sv_pB > 1e-30
ax.loglog(T_range[valid_pB], sv_pB[valid_pB],
label='p-¹¹B', color=REACTIONS['p-¹¹B']['color'], linewidth=2.5,
linestyle='--')
ax.set_xlabel('Temperature (keV)', fontsize=14)
ax.set_ylabel('⟨σv⟩ (cm³/s)', fontsize=14)
ax.set_title('Maxwell-Averaged Fusion Reactivities', fontsize=16, fontweight='bold')
ax.set_xlim(0.5, 1000)
ax.set_ylim(1e-26, 1e-14)
ax.legend(fontsize=13, loc='lower right')
ax.grid(True, which='both', alpha=0.3)
# Add temperature scale in Kelvin on top axis
ax2 = ax.twiny()
ax2.set_xscale('log')
ax2.set_xlim(ax.get_xlim())
kelvin_ticks = [1e7, 1e8, 1e9, 1e10]
ax2.set_xticks([k * 8.617e-8 for k in kelvin_ticks])
ax2.set_xticklabels([f'{k:.0e} K' for k in kelvin_ticks], fontsize=10)
ax2.set_xlabel('Temperature (K)', fontsize=12)
# Annotations
ax.axhline(y=1e-16, color='gray', linestyle=':', alpha=0.5)
ax.text(2, 1.5e-16, 'Reactor-relevant range', fontsize=10, color='gray')
plt.tight_layout()
plt.savefig('fig1_reactivities.png', dpi=150, bbox_inches='tight')
plt.close()
print("\nFigure 1 saved: fig1_reactivities.png")
# =============================================================
# §4. Gamow Peak Visualization
# =============================================================
def gamow_peak(E_keV, T_keV, E_G_keV):
"""Calculate the Gamow peak integrand (unnormalized)."""
return np.exp(-E_keV / T_keV - np.sqrt(E_G_keV / E_keV))
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
temperatures = [5, 10, 20, 50] # keV
E_range = np.linspace(0.1, 200, 2000)
for idx, (ax, T) in enumerate(zip(axes.flat, temperatures)):
for name, r in REACTIONS.items():
peak = gamow_peak(E_range, T, r['E_G'])
peak_norm = peak / np.max(peak) if np.max(peak) > 0 else peak
ax.plot(E_range, peak_norm, label=name, color=r['color'], linewidth=2)
# Maxwell-Boltzmann (unnormalized)
mb = np.exp(-E_range / T)
mb_norm = mb / np.max(mb)
ax.plot(E_range, mb_norm, 'k--', alpha=0.3, label='Maxwell-Boltzmann')
ax.set_xlabel('Energy (keV)', fontsize=11)
ax.set_ylabel('Relative Probability', fontsize=11)
ax.set_title(f'T = {T} keV ({T*1.16e7:.1e} K)', fontsize=12, fontweight='bold')
ax.set_ylim(0, 1.1)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.suptitle('Gamow Peak at Different Temperatures', fontsize=16, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig('fig2_gamow_peaks.png', dpi=150, bbox_inches='tight')
plt.close()
print("Figure 2 saved: fig2_gamow_peaks.png")
# =============================================================
# §5. Gamow Peak Energy E_0 Calculation
# =============================================================
print("\n" + "=" * 60)
print("Gamow Peak Energies E_0 = (E_G * (k_B*T)^2 / 4)^(1/3)")
print("=" * 60)
for T in [1, 5, 10, 20, 50, 100]:
print(f"\nT = {T} keV:")
for name, r in REACTIONS.items():
E0 = (r['E_G'] * T**2 / 4)**(1/3)
dE = 4/np.sqrt(3) * (E0 * T)**0.5
print(f" {name:8s}: E_0 = {E0:8.1f} keV, ΔE = {dE:8.1f} keV, E_0/T = {E0/T:.1f}")
# =============================================================
# §6. Binding Energy Per Nucleon Curve
# =============================================================
# Selected nuclei for B/A curve
nuclei = [
('n', 1, 0, 0),
('¹H', 1, 0, 0),
('²H', 2, 1, 2.225),
('³H', 3, 1, 8.482),
('³He', 3, 2, 7.718),
('⁴He', 4, 2, 28.296),
('⁶Li', 6, 3, 31.995),
('⁷Li', 7, 3, 39.245),
('⁹Be', 9, 4, 58.165),
('¹⁰B', 10, 5, 64.751),
('¹¹B', 11, 5, 76.206),
('¹²C', 12, 6, 92.162),
('¹⁴N', 14, 7, 104.659),
('¹⁶O', 16, 8, 127.619),
('²⁰Ne', 20, 10, 160.645),
('²⁴Mg', 24, 12, 198.257),
('²⁸Si', 28, 14, 236.537),
('³²S', 32, 16, 271.780),
('⁴⁰Ca', 40, 20, 342.052),
('⁵⁶Fe', 56, 26, 492.258),
('⁵⁸Ni', 58, 28, 506.459),
('⁶²Ni', 62, 28, 545.262),
('⁶³Cu', 63, 29, 551.385),
('¹⁰⁷Ag', 107, 47, 915.270),
('¹⁹⁷Au', 197, 79, 1559.40),
('²⁰⁸Pb', 208, 82, 1636.44),
('²³⁵U', 235, 92, 1783.87),
('²³⁸U', 238, 92, 1801.69),
]
fig, ax = plt.subplots(figsize=(14, 7))
A_vals = [n[1] for n in nuclei if n[1] > 1]
BA_vals = [n[3]/n[1] for n in nuclei if n[1] > 1]
names = [n[0] for n in nuclei if n[1] > 1]
ax.plot(A_vals, BA_vals, 'o-', color='#2c3e50', markersize=6, linewidth=1.5)
# Annotate key nuclei
highlight = {'²H': (-15, -20), '³H': (-15, 10), '⁴He': (10, 10),
'⁵⁶Fe': (10, 10), '¹²C': (10, -15), '²³⁵U': (-40, -15),
'¹¹B': (10, -15), '⁶Li': (-15, -15)}
for i, name in enumerate(names):
if name in highlight:
dx, dy = highlight[name]
ax.annotate(name, (A_vals[i], BA_vals[i]),
textcoords="offset points", xytext=(dx, dy),
fontsize=10, fontweight='bold',
arrowprops=dict(arrowstyle='->', color='gray', alpha=0.5))
# Mark fusion and fission regions
ax.axvline(x=56, color='red', linestyle=':', alpha=0.5)
ax.text(25, 1.5, 'FUSION\n(energy released →)', fontsize=14, color='#27ae60',
fontweight='bold', ha='center')
ax.text(150, 1.5, '← FISSION\n(energy released)', fontsize=14, color='#e74c3c',
fontweight='bold', ha='center')
ax.annotate('⁵⁶Fe peak\n(8.79 MeV/nucleon)', xy=(56, 8.79), xytext=(80, 9.2),
fontsize=11, fontweight='bold', color='red',
arrowprops=dict(arrowstyle='->', color='red'))
ax.set_xlabel('Mass Number A', fontsize=14)
ax.set_ylabel('Binding Energy per Nucleon (MeV)', fontsize=14)
ax.set_title('Nuclear Binding Energy Curve', fontsize=16, fontweight='bold')
ax.set_xlim(0, 250)
ax.set_ylim(0, 9.5)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('fig3_binding_energy.png', dpi=150, bbox_inches='tight')
plt.close()
print("Figure 3 saved: fig3_binding_energy.png")
# =============================================================
# §7. Coulomb Barrier Visualization
# =============================================================
fig, ax = plt.subplots(figsize=(12, 7))
r_range = np.linspace(1, 30, 1000) # fm
for name, r_data in REACTIONS.items():
Z1, Z2 = r_data['Z1'], r_data['Z2']
V_coulomb = 1.44 * Z1 * Z2 / r_range # MeV
R_nuclear = 1.2 * (2**(1/3) + (r_data['mu']*2)**(1/3)) # approximate
V_total = np.where(r_range > R_nuclear, V_coulomb, -50 + V_coulomb)
ax.plot(r_range, V_coulomb, label=f'{name} (Z₁Z₂={Z1*Z2})',
color=r_data['color'], linewidth=2.5)
ax.axhline(y=0, color='black', linewidth=0.5)
ax.set_xlabel('Distance r (fm)', fontsize=14)
ax.set_ylabel('Coulomb Potential V(r) (MeV)', fontsize=14)
ax.set_title('Coulomb Barrier for Fusion Reactions', fontsize=16, fontweight='bold')
ax.set_xlim(1, 25)
ax.set_ylim(0, 3)
ax.legend(fontsize=13)
ax.grid(True, alpha=0.3)
# Add kinetic energy lines
for T, label in [(10, '10 keV'), (100, '100 keV')]:
ax.axhline(y=T/1000, color='gray', linestyle='--', alpha=0.5)
ax.text(22, T/1000 + 0.03, f'E = {label}', fontsize=10, color='gray')
plt.tight_layout()
plt.savefig('fig4_coulomb_barrier.png', dpi=150, bbox_inches='tight')
plt.close()
print("Figure 4 saved: fig4_coulomb_barrier.png")
# =============================================================
# §8. Summary Statistics
# =============================================================
print("\n" + "=" * 60)
print("SUMMARY: Reactivity at Key Temperatures")
print("=" * 60)
print(f"{'T (keV)':<10} {'D-T':<14} {'D-D':<14} {'D-³He':<14} {'p-¹¹B':<14}")
print("-" * 66)
for T in [1, 5, 10, 20, 50, 100, 200, 500, 1000]:
sv_dt = reactivity_bosch_hale(T, BH_PARAMS['D-T'])
sv_dd = reactivity_bosch_hale(T, BH_PARAMS['D-D'])
sv_dhe3 = reactivity_bosch_hale(T, BH_PARAMS['D-³He'])
sv_pb = reactivity_pB11(T)
print(f"{T:<10} {sv_dt:<14.3e} {sv_dd:<14.3e} {sv_dhe3:<14.3e} {sv_pb:<14.3e}")
# =============================================================
# §9. Tunneling Probability Comparison
# =============================================================
print("\n" + "=" * 60)
print("Tunneling Probability exp(-sqrt(E_G/E))")
print("=" * 60)
for E in [1, 5, 10, 50, 100, 500]:
print(f"\nE = {E} keV:")
for name, r in REACTIONS.items():
arg = np.sqrt(r['E_G'] / E)
P = np.exp(-arg)
print(f" {name:8s}: sqrt(E_G/E) = {arg:8.2f}, P = {P:.3e}")
print("\n" + "=" * 60)
print("All computations complete. Figures saved.")
print("=" * 60)
第II部：プラズマ閉じ込め物理
§12. なぜ閉じ込めが中心問題なのか
第I部は核融合に10〜300 keV（10⁸〜10⁹ K）の温度が必要であることを示した。この温度では物質は完全電離したプラズマ — 電磁力に支配される自由イオンと電子のガスとなる。
いかなる物質容器も1億度のプラズマを保持できない。壁は接触と同時に蒸発する。核融合工学の中心問題は、したがって以下である：核融合反応が有用な速度で起こるのに十分な時間、プラズマを物質壁から離して保持するにはどうすればよいか。
ローソン条件（本シリーズVol.2で導出）はこれを定量化する：プラズマ密度 $n$、閉じ込め時間 $\tau_E$、温度 $T$ の積が閾値を超えなければならない：
$$
n \tau_E T \geq 3 \times 10^{21} \text{ keV·s/m}^3 \quad \text{（D-T点火条件）}
$$
二つの基本的アプローチが存在する：
| アプローチ | 戦略 | 例 | n (m⁻³) | τ_E (s) |
|---|---|---|---|---|
| 磁場閉じ込め | 低密度、長時間 | トカマク (ITER) | ~10²⁰ | ~3–5 |
| 慣性閉じ込め | 高密度、短時間 | レーザー (NIF) | ~10³¹ | ~10⁻¹¹ |
本パートでは磁場閉じ込めに焦点を当てる。これは現在のすべての炉設計の基盤であり、核融合推進への最も実現可能な道筋である。
§13. 電磁流体力学 — プラズマの流体モデル
§13.1 MHD方程式
電磁流体力学（MHD）はプラズマを電磁場と結合した導電性流体として扱う。理想MHD方程式は以下の通り：
質量保存：
$$
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0
$$
運動量（力のバランス）：
$$
\rho \frac{d\mathbf{v}}{dt} = \mathbf{J} \times \mathbf{B} - \nabla p
$$
オームの法則（理想 — 無限導電率）：
$$
\mathbf{E} + \mathbf{v} \times \mathbf{B} = 0
$$
ファラデーの法則：
$$
\frac{\partial \mathbf{B}}{\partial t} = -\nabla \times \mathbf{E} = \nabla \times (\mathbf{v} \times \mathbf{B})
$$
アンペールの法則（低周波極限）：
$$
\nabla \times \mathbf{B} = \mu_0 \mathbf{J}
$$
エネルギー（断熱）：
$$
\frac{d}{dt}\left(\frac{p}{\rho^\gamma}\right) = 0
$$
各記号：$\rho$ = 質量密度、$\mathbf{v}$ = 流体速度、$p$ = 圧力、$\mathbf{B}$ = 磁場、$\mathbf{J}$ = 電流密度、$\gamma = 5/3$。
§13.2 凍結磁束定理
理想オームの法則から、重要な帰結が得られる：**磁力線はプラズマ流体に「凍結」されている。**プラズマが移動すれば磁力線を引きずり、磁力線が移動すればプラズマを引きずる。
これが磁場閉じ込めの基盤である：出口のない磁場配位を作れば、（磁場に凍結された）プラズマもまた出口がない。
§13.3 磁気圧力とベータ
磁場は圧力を及ぼす：
$$
p_B = \frac{B^2}{2\mu_0}
$$
プラズマ圧力と磁気圧力の比がベータである：
$$
\beta = \frac{p}{B^2/(2\mu_0)} = \frac{n k_B T}{B^2/(2\mu_0)}
$$
| パラメータ | 意味 |
|---|---|
| $\beta \ll 1$ | 磁場が支配的；プラズマは良好に閉じ込められるがパワー密度は低い |
| $\beta \sim 1$ | プラズマ圧と磁気圧が同程度；高パワー密度だが安定化が困難 |
| $\beta > 1$ | プラズマ圧が磁気圧を超過；閉じ込め破綻 |
典型値：
| 装置 | β (%) |
|---|---|
| トカマク（従来型） | 3–5 |
| トカマク（先進型） | 5–12 |
| ステラレーター | 4–5 |
| 磁場反転配位 (FRC) | 50–90 |
| 球状トカマク | 15–40 |
**推進用途には高β配位（FRC、球状トカマク）が強く好まれる。**単位磁場あたりのパワー密度が高く、磁石質量を削減できるためである — これは宇宙機にとって決定的な制約である（Vol.8）。
§14. MHD平衡 — Grad-Shafranov方程式
§14.1 静的平衡条件
静的平衡（$\mathbf{v} = 0$、$\partial/\partial t = 0$）のプラズマに対し、MHD運動方程式は以下に帰着する：
$$
\mathbf{J} \times \mathbf{B} = \nabla p
$$
これが基本的な力のバランスである：外向きの圧力勾配が内向きの $\mathbf{J} \times \mathbf{B}$ 力で釣り合う。
$\mathbf{B}$ との内積をとると：
$$
\mathbf{B} \cdot \nabla p = 0
$$
**圧力は磁力線に沿って一定である。**等圧面（等圧線）は磁気面でもあり、磁力線はその面内に完全に含まれる。
$\mathbf{J}$ との内積をとると：
$$
\mathbf{J} \cdot \nabla p = 0
$$
電流もまた磁気面内に存在する。
§14.2 Grad-Shafranov方程式の導出
軸対称系（トカマク）に対し、円筒座標 $(R, \phi, Z)$ を用いて磁場は以下のように書ける：
$$
\mathbf{B} = \frac{1}{R}\nabla\Psi \times \hat{\phi} + \frac{F(\Psi)}{R}\hat{\phi}
$$
各記号：
- $\Psi(R,Z)$ = ポロイダル磁束関数（磁気面のラベル）
- $F(\Psi) = RB_\phi$ = 反磁性関数
力のバランス方程式に代入し、アンペールの法則を用いると：
$$
R \frac{\partial}{\partial R}\left(\frac{1}{R}\frac{\partial \Psi}{\partial R}\right) + \frac{\partial^2 \Psi}{\partial Z^2} = -\mu_0 R^2 \frac{dp}{d\Psi} - F\frac{dF}{d\Psi}
$$
これがGrad-Shafranov方程式（GSE）である — 軸対称MHD平衡の基本方程式。$\Psi(R,Z)$ に対する非線形楕円型偏微分方程式であり、右辺は自由関数 $p(\Psi)$ と $F(\Psi)$ で指定される。
§14.3 Grad-Shafranov演算子
Grad-Shafranov演算子を定義する：
$$
\Delta^* \Psi \equiv R \frac{\partial}{\partial R}\left(\frac{1}{R}\frac{\partial \Psi}{\partial R}\right) + \frac{\partial^2 \Psi}{\partial Z^2}
$$
GSEは以下となる：
$$
\Delta^* \Psi = -\mu_0 R^2 p'(\Psi) - F(\Psi) F'(\Psi)
$$
各項の物理的解釈：
| 項 | 物理 |
|---|---|
| $\Delta^* \Psi$ | ポロイダル電流密度（トロイダル力） |
| $-\mu_0 R^2 p'$ | プラズマ圧力勾配の駆動 |
| $-FF'$ | トロイダル磁場勾配の駆動 |
§14.4 安全係数 q
安全係数は、磁力線がポロイダル方向に1周する間にトロイダル方向に何周するかを測る：
$$
q(\Psi) = \frac{1}{2\pi} \oint \frac{\mathbf{B} \cdot \nabla\phi}{\mathbf{B} \cdot \nabla\theta} d\theta = \frac{1}{2\pi} \oint \frac{F}{R^2} \frac{dl_p}{|\nabla\Psi|/R}
$$
積分は磁束面 $\Psi$ 上のポロイダル周回にわたる。
qが重要な理由：$q = m/n$（整数）の有理面では磁力線が自身に戻る。これらは摂動が各周回後に自己強化できるため、不安定性に対して最も脆弱である。最も危険な面は $q = 1$（キンク不安定性）と $q = 2$（ティアリングモード）である。
典型的なトカマクプロファイル：
- 中心：$q_0 \approx 0.8$–1.2
- 端：$q_a \approx 3$–7
- 安定性のためには一般に $q$ は半径とともに増加する必要がある
§14.5 シャフラノフシフト
トーラスにおいて、外向きの圧力勾配とトロイダル電流の遠心効果により、磁気面は外側（$R$ が大きい方向）にシフトする。これがシャフラノフシフト $\Delta_S$ である：
$$
\Delta_S \approx \frac{\beta_p a^2}{2R_0}
$$
ここで $\beta_p$ = ポロイダルベータ、$a$ = 小半径、$R_0$ = 大半径。
高βプラズマではシャフラノフシフトが大きくなり、磁気面を外側壁に向かって押しやり閉じ込めを劣化させる。これがトカマクにおける達成可能βの根本的制限である。
§15. 磁場閉じ込め配位
§15.1 トカマク
トカマク（ロシア語：тороидальная камера с магнитными катушками）は最も開発が進んだ磁場閉じ込め概念である。
**形状：**トーラス（ドーナツ型）で以下を持つ：
- トロイダル磁場 $B_\phi$：外部コイルにより生成；一次閉じ込めを提供
- ポロイダル磁場 $B_\theta$：プラズマ電流 $I_p$ により生成；回転変換を提供
- 螺旋磁力線：トロイダル磁場とポロイダル磁場の組合せが入れ子状のトロイダル面上に螺旋経路を描く
主要パラメータ：
| パラメータ | 記号 | ITER | SPARC |
|---|---|---|---|
| 大半径 | $R_0$ | 6.2 m | 1.85 m |
| 小半径 | $a$ | 2.0 m | 0.57 m |
| アスペクト比 | $A = R_0/a$ | 3.1 | 3.25 |
| トロイダル磁場 | $B_0$ | 5.3 T | 12.2 T |
| プラズマ電流 | $I_p$ | 15 MA | 8.7 MA |
| パルス長 | — | 400 s | 10 s |
| 目標Q | — | ≥10 | ≥2 |
**利点：**最も成熟した設計。最高の三重積達成実績。Q > 1への明確な道筋。
**欠点：**本質的にパルス運転（電流駆動は損失的）。ディスラプション。従来型磁石では大型化。
**SPARC/ARCの革新：**MITのSPARCは高温超伝導（HTS）磁石（REBCOテープ）を使用し、コイル上で20 Tを達成。これにより劇的に小型の装置が可能となる。$\beta \propto nT/B^2$ の関係から、$B$ を2倍にすれば同じβで達成可能なプラズマ圧力が4倍になるか、同性能を1/4の体積で実現できる。
§15.2 ステラレーター
ステラレーターは閉じ込め磁場全体を外部から生成する（プラズマ電流不要）。
**形状：**複雑な3Dコイルがプラズマ電流に依存せず螺旋磁場を生成する。必要な回転変換は外部コイル形状に組み込まれている。
主要装置：Wendelstein 7-X (W7-X)
| パラメータ | 値 |
|---|---|
| 大半径 | 5.5 m |
| 小半径 | 0.53 m |
| 磁場 | 3.0 T |
| 加熱パワー | 10 MW |
| パルス長 | 最大30分（目標） |
**利点：**本質的に定常運転（プラズマ電流不要）。ディスラプションなし。電流駆動パワー不要。
**欠点：**極めて複雑な3Dコイル設計。歴史的にトカマクより閉じ込めが劣る（最適化により改善中）。大型化。
§15.3 磁場反転配位 (FRC)
**形状：**純粋にポロイダル磁場のみを持つコンパクトトロイド。トロイダル磁場なし。プラズマ電流が自己完結的な磁場構造を生成する。
| パラメータ | 推進への意味 |
|---|---|
| $\beta \sim 50$–90% | 単位磁場あたり最高のパワー密度 |
| コンパクト形状 | 宇宙機に対し低質量 |
| 天然ダイバータ | 排気が軸方向に指向 → 磁気ノズル |
| 移動可能 | 生成、圧縮、射出が可能 |
主要装置：TAE Technologies (C-2W)
- 中性粒子ビーム入射によりFRCを維持
- FRC寿命 > 30 ms、$T_i > 1$ keVを達成
- p-¹¹Bを長期目標燃料とする
FRCは核融合推進に最も有望な配位である（Vol.8–9）。理由は以下の通り：
- 軸方向排気形状が磁気ノズルに自然に結合する
- 高βにより磁石質量を最小化
- コンパクトなフォームファクターが宇宙機の質量バジェットに適合
- FRCプラズモイドは生成、加速、射出が可能（パルス推進）
§15.4 その他の配位
| 配位 | 主な特徴 | 現状 |
|---|---|---|
| 球状トカマク | 低アスペクト比（$A \sim 1.5$）、高β | MAST-U, NSTX-U |
| 逆転磁場ピンチ | $B_\phi$ が端で反転；高β | RFX-mod |
| Zピンチ | 軸方向電流がピンチを生成；パルス型 | Zap Energy（剪断流型） |
| ミラー型 | ミラーコイルを持つ開放磁力線 | 大部分が放棄；復活の試み |
| 磁化標的 | FRC/スフェロマクの圧縮 | General Fusion |
§16. MHD不安定性 — プラズマはなぜ逃げるのか
§16.1 エネルギー原理
プラズマ平衡が安定であるとは、あらゆる微小摂動が全ポテンシャルエネルギーを増加させることである。形式的には：
$$
\delta W = -\frac{1}{2} \int \boldsymbol{\xi}^* \cdot \mathbf{F}(\boldsymbol{\xi}) , d^3x > 0 \quad \text{すべての } \boldsymbol{\xi} \text{ に対して}
$$
ここで $\boldsymbol{\xi}$ はプラズマ変位、$\mathbf{F}$ はMHD力の演算子である。任意の摂動に対して $\delta W < 0$ ならば、平衡は不安定である。
ポテンシャルエネルギー汎関数：
$$
\delta W = \frac{1}{2\mu_0} \int \left[ |\delta\mathbf{B}\perp|^2 + B^2|\nabla\cdot\boldsymbol{\xi}\perp + 2\boldsymbol{\xi}\perp \cdot \boldsymbol{\kappa}|^2 + \gamma p |\nabla\cdot\boldsymbol{\xi}|^2 \right.
$$
$$
\left. - 2(\boldsymbol{\xi}\perp \cdot \nabla p)(\boldsymbol{\xi}\perp^* \cdot \boldsymbol{\kappa}) - J\parallel (\boldsymbol{\xi}\perp^* \times \hat{b}) \cdot \delta\mathbf{B}\perp \right] d^3x
$$
| 項 | 符号 | 物理的意味 |
|---|---|---|
| $ | \delta\mathbf{B}_\perp | ^2$ |
| $B^2 | \nabla\cdot\boldsymbol{\xi}_\perp + ... | ^2$ |
| $\gamma p | \nabla\cdot\boldsymbol{\xi} | ^2$ |
| $-2(\boldsymbol{\xi}\perp\cdot\nabla p)(\boldsymbol{\xi}\perp^*\cdot\boldsymbol{\kappa})$ | −（不安定化） | 圧力勾配 + 曲率駆動 |
| $-J_\parallel(...)$ | −（不安定化） | 電流駆動 |
最後の二項が核融合プラズマにおけるすべてのMHD不安定性を駆動する。
§16.2 MHD不安定性の分類
摂動はトロイダル面上のモード数 $(m, n)$ で特徴づけられる：
$$
\boldsymbol{\xi} \propto \exp(im\theta - in\phi)
$$
ここで $m$ = ポロイダルモード数、$n$ = トロイダルモード数。
電流駆動不安定性：
| 不安定性 | モード | 条件 | 結果 |
|---|---|---|---|
| 外部キンク | $m=1, n=1$ | $q_a < m/n$（クルスカル-シャフラノフ条件） | 全体的プラズマ変位；ディスラプション |
| 内部キンク | $m=1, n=1$ | $q_0 < 1$ | ソートゥース振動；周期的コア崩壊 |
| ティアリングモード | $m/n = q_{有理}$ | 有限抵抗率が凍結磁束を破る | 磁気島；閉じ込め劣化 |
圧力駆動不安定性：
| 不安定性 | モード | 条件 | 結果 |
|---|---|---|---|
| バルーニング | 高$n$ | 悪い曲率領域で$\nabla p$が臨界勾配を超過 | β限界；圧力プロファイルの平坦化 |
| 交換不安定性 | 低$n$ | 不利な曲率における圧力勾配 | レイリー-テイラー型と類似 |
| ELM（端部局在モード） | ピーリング-バルーニング結合 | 端部ペデスタルが安定性境界を超過 | 周期的端部噴出；壁損傷 |
§16.3 トロイヨンベータ限界
トカマクにおける達成可能最大βは経験的に以下で与えられる：
$$
\beta_{max} (%) = \beta_N \frac{I_p \text{(MA)}}{a \text{(m)} \cdot B_0 \text{(T)}}
$$
ここで $\beta_N \approx 2.8$–3.5 は規格化ベータである。この限界は圧力駆動（バルーニング/キンク）不安定性の発生から生じる。
ITERの場合：$\beta_{max} \approx 2.8 \times 15/(2.0 \times 5.3) \approx 4.0%$
**含意：**達成可能なプラズマ圧力（したがって核融合パワー密度）はMHD安定性により直接制限される。高磁場磁石（SPARCアプローチ）は分母 $a \cdot B_0$ を増加させながら同じβを維持し、より高い絶対圧力を達成する。
§16.4 新古典ティアリングモード (NTM)
実際のトカマクでは、ブートストラップ電流（トロイダル幾何学における圧力勾配から自己生成される電流）が脆弱性を生む：磁気島が形成されると局所的に圧力勾配を平坦化し、島内のブートストラップ電流を減少させ、これが島を拡大する — 正のフィードバックループ。
NTMは現代のトカマクにおける性能劣化の最も一般的な原因であり、ディスラプションの引き金となりうる。島のO点に照射される局所的電子サイクロトロン電流駆動（ECCD）により安定化される。
§17. 輸送理論 — エネルギーはいかに漏れるか
§17.1 古典輸送
一様磁場中で、荷電粒子はラーモア半径で磁力線の周りを旋回する：
$$
\rho_L = \frac{mv_\perp}{|q|B}
$$
トカマクにおけるイオン（$T = 10$ keV、$B = 5$ T）の場合：
$$
\rho_{L,i} \approx 3 \text{ mm}
$$
古典的拡散（粒子衝突で約$\rho_L$だけ変位）：
$$
D_{cl} = \frac{\rho_L^2}{\tau_{col}} \sim \frac{\rho_{L,i}^2 \nu_{ii}}{2}
$$
トカマクパラメータでは：$D_{cl} \sim 10^{-4}$ m²/s。これにより閉じ込め時間は約10³秒 — 必要値をはるかに上回る。
**現実：観測される輸送は古典の10〜100倍悪い。**この食い違いが核融合物理の中心的未解決問題を定義する。
§17.2 新古典輸送
トロイダル幾何学では、粒子はバナナ軌道（不均一磁場のミラー効果で反射される粒子による）を描く。捕捉粒子（磁場の変動で反射される粒子）の有効ステップサイズは：
$$
\Delta r_{banana} \sim \frac{q \rho_L}{\sqrt{\epsilon}}
$$
ここで $\epsilon = r/R_0$ は逆アスペクト比。
新古典拡散：
$$
D_{neo} \sim q^2 \epsilon^{-3/2} D_{cl}
$$
典型的トカマクパラメータでは：$D_{neo} \sim 10^{-2}$ m²/s、古典の約100倍。
しかし観測される輸送は新古典よりさらに約10倍悪い。この追加輸送は異常輸送と呼ばれ、プラズマ乱流により駆動される。
§17.3 異常輸送と乱流
乱流輸送を駆動する主要なミクロ不安定性：
| 不安定性 | スケール | 駆動源 | 影響対象 |
|---|---|---|---|
| イオン温度勾配 (ITG) | $\rho_i$ | $\nabla T_i$ | イオン熱輸送 |
| 捕捉電子モード (TEM) | $\rho_i$ | $\nabla n$, $\nabla T_e$ | 電子 + 粒子輸送 |
| 電子温度勾配 (ETG) | $\rho_e$ | $\nabla T_e$ | 電子熱輸送 |
これらの不安定性はラーモア半径スケールの乱流渦を生成し、衝突過程よりはるかに速く磁束面を横切ってエネルギーと粒子を輸送する。
乱流輸送のモデリングの最先端はジャイロカイネティックシミュレーションであり、ジャイロ平均近似のもとでヴラソフ方程式を5次元位相空間（空間3 + 速度2次元）で解く：
$$
\frac{\partial f}{\partial t} + \dot{\mathbf{R}} \cdot \nabla f + \dot{v}\parallel \frac{\partial f}{\partial v\parallel} = C[f]
$$
ここで $f(\mathbf{R}, v_\parallel, \mu, t)$ はジャイロセンター分布関数、$C[f]$ は衝突演算子。
**計算コスト：**トカマク断面の単一ジャイロカイネティックシミュレーションに約10⁶ CPU時間が必要。結合ジャイロカイネティック乱流による全装置モデリングは現在の計算能力を超えている。ここがAIの最大のインパクトを持つ領域である（本シリーズVol.5）。
§17.4 エネルギー閉じ込め時間のスケーリング
経験的に、トカマクにおけるエネルギー閉じ込め時間はIPB98(y,2)スケーリングに従う：
$$
\tau_E = 0.0562 \cdot I_p^{0.93} B^{0.15} n_{19}^{0.41} P^{-0.69} R^{1.97} \epsilon^{0.58} \kappa^{0.78} M^{0.19}
$$
| 変数 | 意味 |
|---|---|
| $I_p$ (MA) | プラズマ電流 |
| $B$ (T) | トロイダル磁場 |
| $n_{19}$ (10¹⁹ m⁻³) | 線平均密度 |
| $P$ (MW) | 加熱パワー |
| $R$ (m) | 大半径 |
| $\epsilon$ | 逆アスペクト比 |
| $\kappa$ | 楕円度 |
| $M$ | 有効イオン質量 (AMU) |
決定的な観察：$\tau_E \propto R^{1.97}$ — 閉じ込めは装置サイズのほぼ二乗に比例して改善する。ITERが大型（$R = 6.2$ m）である理由はこれであり、コンパクト装置が補償のためにはるかに強い磁場を必要とする理由でもある。
ITERの場合：$\tau_E \approx 3.7$ s（予測値）
§18. ディスラプション — 壊滅的破壊モード
§18.1 ディスラプションとは何か
ディスラプションとは、トカマクにおけるプラズマの閉じ込めと電流の突然の喪失であり、ミリ秒のタイムスケールで発生する。トカマク炉における最も深刻な運転リスクである。
シーケンス：
- **前兆：**MHD不安定性の成長（ロックドティアリングモード、鉛直変位、または密度/電流限界の超過）
- **熱クエンチ：**約1 msでエネルギー閉じ込めが崩壊。プラズマ全熱エネルギー（ITERで約350 MJ）がプラズマ対向部品に投入
- **電流クエンチ：**プラズマ電流が約10–50 msで減衰。真空容器に大規模な渦電流と機械的力を誘起
- **暴走電子：**ディスラプション後の電場が電子を相対論的エネルギー（約10 MeV）まで加速し、暴走ビームを形成。これが第一壁を損傷しうる
§18.2 ディスラプションの影響
ITERスケールの装置に対して：
| 効果 | 規模 |
|---|---|
| 壁への熱負荷 | 約1 msで最大60 MJ/m² |
| 電磁力 | 真空容器に約100 MN |
| 暴走電子エネルギー | 最大12 MJ/m²（局所化） |
| ハロー電流力 | 非対称、$I_p$ の最大25% |
**比較のために：**1 msで60 MJ/m²は、1平方メートルあたり約14 kgのTNTを爆発させるのと等価である。ITERがディスラプション緩和システム（シャッタードペレット注入）を有する理由はこれであり、エネルギーが壁に到達する前に放射的に散逸させるよう設計されている。
§18.3 ディスラプションの回避と緩和
| 戦略 | メカニズム | 現状 |
|---|---|---|
| 回避 | 安定限界内での運転；リアルタイム制御 | 活発な分野；ML予測器（Vol.5） |
| 予測 | ディスラプション10–100 ms前に前兆を検出 | DeepMind、MIT：MLモデルが>95%精度を達成 |
| 緩和 | 不純物注入によるエネルギーの均一放射 | ITER SPIシステム；大量ガス注入 |
| 暴走抑制 | 相対論的電子ビーム形成の防止 | 重水素SPI；共鳴磁場摂動 |
**ディスラプションはステラレーターとFRCの最も強い論拠である。**これらは閉じ込めにプラズマ電流を使用しないため、本質的にディスラプションが起こらない。推進システム（Vol.8）にとって、ディスラプション免疫は不可欠である — 宇宙機エンジンにおけるディスラプションは壊滅的故障を意味する。
§19. 閉じ込め計算解析
§19.1 Grad-Shafranov平衡ソルバー
"""
Grad-Shafranov Equation Solver — Simplified Tokamak Equilibrium
================================================================
Solves ΔΨ = -μ₀R²p'(Ψ) - FF'(Ψ) on a rectangular grid
using iterative relaxation (Picard iteration).
Author: Dosanko Tousan + Claude (Anthropic)
License: MIT
"""
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
# ============ Grid Setup ============
NR, NZ = 129, 129
R_min, R_max = 3.0, 9.0 # meters (ITER-like)
Z_min, Z_max = -4.0, 4.0
R = np.linspace(R_min, R_max, NR)
Z = np.linspace(Z_min, Z_max, NZ)
dR = R[1] - R[0]
dZ = Z[1] - Z[0]
RR, ZZ = np.meshgrid(R, Z, indexing='ij')
# ============ Equilibrium Parameters ============
R0 = 6.2 # Major radius (m) - ITER
a = 2.0 # Minor radius (m) - ITER
kappa = 1.7 # Elongation
B0 = 5.3 # Toroidal field at R0 (T)
mu0 = 4 * np.pi * 1e-7
# ============ Source Functions (Solov'ev-type analytic profiles) ============
# p'(Ψ) = C1, FF'(Ψ) = C2 (constants → linear profiles)
# This gives an analytic solution class (Solov'ev equilibrium)
p0 = 6e5 # Peak pressure (Pa) ~ 6 bar for ITER-like
C1 = -p0 / (a**2) # Pressure gradient coefficient
C2 = -(B0**2 * R0**2 * 0.1) / (a**2) # FF' coefficient
# ============ Solov'ev Analytic Solution (for comparison) ============
def psi_solovev(R, Z, R0, a, kappa):
"""Approximate Solov'ev equilibrium."""
x = (R - R0) / a
y = Z / (kappa * a)
return (x**2 + y**2 - 1) * a**2
PSI_analytic = psi_solovev(RR, ZZ, R0, a, kappa)
# ============ Iterative GS Solver ============
PSI = np.zeros((NR, NZ))
PSI_new = np.zeros_like(PSI)
# Boundary condition: Ψ = 0 on boundary (vacuum)
max_iter = 5000
tol = 1e-6
for iteration in range(max_iter):
for i in range(1, NR-1):
for j in range(1, NZ-1):
# Grad-Shafranov operator in finite differences
# Δ*Ψ = R ∂/∂R(1/R ∂Ψ/∂R) + ∂²Ψ/∂Z²
d2R = (PSI[i+1,j] - 2*PSI[i,j] + PSI[i-1,j]) / dR**2
dR_term = (1.0/(2*R[i]*dR)) * (PSI[i+1,j] - PSI[i-1,j])
GS_R = d2R - dR_term # This is Δ* in R direction
d2Z = (PSI[i,j+1] - 2*PSI[i,j] + PSI[i,j-1]) / dZ**2
# Source term: -μ₀R²p' - FF'
# Inside plasma (Ψ < 0 in our convention for Solov'ev):
r_here = R[i]
rho2 = ((r_here - R0)/a)**2 + (ZZ[i,j]/(kappa*a))**2
if rho2 < 1.0: # Inside plasma boundary
source = -mu0 * r_here**2 * C1 - C2
else:
source = 0.0
# SOR update
PSI_new[i,j] = (1.0/(2/dR**2 + 2/dZ**2)) * (
(PSI[i+1,j] + PSI[i-1,j])/dR**2 +
(PSI[i,j+1] + PSI[i,j-1])/dZ**2 -
(PSI[i+1,j] - PSI[i-1,j])/(2*r_here*dR) -
source
)
# Check convergence
diff = np.max(np.abs(PSI_new[1:-1,1:-1] - PSI[1:-1,1:-1]))
PSI[1:-1,1:-1] = PSI_new[1:-1,1:-1]
if diff < tol:
print(f"GS solver converged in {iteration+1} iterations (residual: {diff:.2e})")
break
if (iteration+1) % 1000 == 0:
print(f" Iteration {iteration+1}: residual = {diff:.2e}")
if iteration == max_iter - 1:
print(f"GS solver: max iterations reached (residual: {diff:.2e})")
# ============ Plot Equilibrium ============
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
# Numerical solution
ax1 = axes[0]
levels = np.linspace(PSI.min(), PSI.max(), 30)
cs1 = ax1.contour(RR, ZZ, PSI, levels=levels, cmap='RdBu_r', linewidths=0.8)
ax1.set_xlabel('R (m)', fontsize=13)
ax1.set_ylabel('Z (m)', fontsize=13)
ax1.set_title('Numerical GS Equilibrium', fontsize=14, fontweight='bold')
ax1.set_aspect('equal')
ax1.grid(True, alpha=0.3)
# Plot plasma boundary
theta = np.linspace(0, 2*np.pi, 200)
R_boundary = R0 + a * np.cos(theta)
Z_boundary = kappa * a * np.sin(theta)
ax1.plot(R_boundary, Z_boundary, 'k-', linewidth=2, label='Plasma boundary')
ax1.legend(fontsize=11)
# Safety factor profile (approximate)
ax2 = axes[1]
r_norm = np.linspace(0.01, 1.0, 100)
# Approximate q-profile: q(r) = q0 + (qa - q0) * r^2
q0 = 1.0
qa = 3.5
q_profile = q0 + (qa - q0) * r_norm**2
ax2.plot(r_norm, q_profile, 'b-', linewidth=2.5)
ax2.axhline(y=1, color='red', linestyle='--', alpha=0.5, label='q = 1 (sawtooth)')
ax2.axhline(y=1.5, color='orange', linestyle='--', alpha=0.5, label='q = 3/2 (NTM)')
ax2.axhline(y=2, color='green', linestyle='--', alpha=0.5, label='q = 2 (tearing)')
ax2.set_xlabel('Normalized radius r/a', fontsize=13)
ax2.set_ylabel('Safety factor q', fontsize=13)
ax2.set_title('Safety Factor Profile', fontsize=14, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.3)
ax2.set_ylim(0, 5)
plt.tight_layout()
plt.savefig('fig5_gs_equilibrium.png', dpi=150, bbox_inches='tight')
plt.close()
print("Figure 5 saved: fig5_gs_equilibrium.png")
# ============ Transport Comparison ============
fig, ax = plt.subplots(figsize=(12, 7))
r_a = np.linspace(0.1, 0.95, 100)
# Classical diffusion (order of magnitude)
D_classical = 1e-4 * np.ones_like(r_a) # m²/s
# Neoclassical (banana regime, simplified)
q_prof = q0 + (qa - q0) * r_a**2
epsilon = r_a * a / R0
D_neoclassical = D_classical * q_prof**2 * epsilon**(-1.5)
# Anomalous (empirical, Bohm-like scaling for illustration)
T_profile = 10 * (1 - r_a**2)**2 # keV, parabolic
B_field = 5.3 # T
# Bohm diffusion: D_B = T/(16eB)
D_bohm = T_profile * 1000 * 1.602e-19 / (16 * 1.602e-19 * B_field) # m²/s
D_anomalous = 0.3 * D_bohm # GyroBohm-like
ax.semilogy(r_a, D_classical, 'g-', linewidth=2.5, label='Classical')
ax.semilogy(r_a, D_neoclassical, 'b-', linewidth=2.5, label='Neoclassical')
ax.semilogy(r_a, D_anomalous, 'r-', linewidth=2.5, label='Anomalous (observed)')
ax.fill_between(r_a, D_neoclassical, D_anomalous, alpha=0.15, color='red',
label='Anomalous enhancement')
ax.set_xlabel('Normalized radius r/a', fontsize=14)
ax.set_ylabel('Diffusion coefficient D (m²/s)', fontsize=14)
ax.set_title('Transport Regimes in a Tokamak', fontsize=16, fontweight='bold')
ax.legend(fontsize=13)
ax.grid(True, which='both', alpha=0.3)
ax.set_ylim(1e-5, 10)
plt.tight_layout()
plt.savefig('fig6_transport.png', dpi=150, bbox_inches='tight')
plt.close()
print("Figure 6 saved: fig6_transport.png")
# ============ Beta Limit Diagram ============
fig, ax = plt.subplots(figsize=(12, 7))
Ip_range = np.linspace(1, 20, 100) # MA
beta_N_values = [2.0, 2.8, 3.5, 4.0]
a_val = 2.0 # m
B_values = [3.0, 5.3, 12.0, 20.0] # T
for B_val in B_values:
beta_max = 2.8 * Ip_range / (a_val * B_val)
label = f'B = {B_val} T' + (' (ITER)' if B_val == 5.3 else '') + (' (SPARC HTS)' if B_val == 12.0 else '')
ax.plot(Ip_range, beta_max, linewidth=2.5, label=label)
# Mark ITER operating point
ax.plot(15, 2.8*15/(2.0*5.3), 'ko', markersize=12, label='ITER design point')
ax.plot(8.7, 2.8*8.7/(2.0*12.0), 's', color='purple', markersize=12, label='SPARC design point')
ax.set_xlabel('Plasma Current Ip (MA)', fontsize=14)
ax.set_ylabel('Maximum β (%)', fontsize=14)
ax.set_title('Troyon Beta Limit: β_max = β_N × I_p / (a·B)', fontsize=16, fontweight='bold')
ax.legend(fontsize=12)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 15)
plt.tight_layout()
plt.savefig('fig7_beta_limit.png', dpi=150, bbox_inches='tight')
plt.close()
print("Figure 7 saved: fig7_beta_limit.png")
§20. 炉設計および推進システムへの含意
§20.1 エネルギー炉に対して（Vol.2–4）
第I部と第II部の統合的物理解析により、以下の制約が確立される：
-
**D-Tが唯一の実現可能な第一世代燃料である。**達成可能な温度における約100倍の反応度優位性が、これを不可避の出発点とする。
-
**14 MeV中性子問題が工学を支配する。**D-Tエネルギー出力の80%は磁場で閉じ込められない中性子である（Vol.2–3）。
-
**閉じ込めは乱流輸送により制限される。**異常輸送は新古典の約10倍であり、第一原理予測モデルは存在しない。このギャップこそAIが最もインパクトを持つ領域である（Vol.5）。
-
**トロイヨンβ限界がパワー密度を制約する。**高磁場磁石（HTS）はコンパクト装置における高性能への最も直接的な道筋である。
-
**ディスラプションはトカマク炉にとって存在論的リスクである。**緩和システムは複雑性を加える；ステラレーターとFRCは問題そのものを回避する。
§20.2 推進システムに対して（Vol.6–9）
-
**FRCが核融合推進の有力候補である。**高β、コンパクト形状、軸方向排気、ディスラプション免疫。
-
**非中性子または低中性子反応が必要である。**中性子は推進剤エネルギーを浪費し、遮蔽質量のペナルティを生む。
-
比推力の優位性は革命的である。$I_{sp} \sim 10^6$ s（核融合）対 ~450 s（化学）。外惑星およびそれ以遠へのミッションを可能にする。
-
**バルキリーへの道は三つの技術を経由する：**コンパクトFRC閉じ込め、非中性子燃料（p-¹¹BまたはD-³He）、磁気ノズル排気変換。
§21. 不確実性と限界
第I部の限界
- 断面積データの不確実性：D-Tで約2%、p-¹¹Bの共鳴領域で約10–20%。
- 反応度計算はマクスウェル分布プラズマを仮定している。
- p-¹¹Bのパラメータ化は簡略化された適合である。
第II部の限界
- 提示したGrad-Shafranovソルバーは簡略化されたSolov'ev平衡を使用。本番級の平衡コード（EFIT、VMEC）は測定磁場データを再構成に使用する。
- IPB98(y,2)スケーリング則は約15%の散乱を伴う経験式である。ITER予測には相当の不確実性がある。
- ジャイロカイネティック輸送モデリングは全装置シミュレーションには計算的に実行不能；実用には縮約モデル（TGLF、QuaLiKiz）が使用され、忠実度の低下を伴う。
- ディスラプション予測は確率的のまま；すべてのシナリオに対する決定論的回避アルゴリズムは存在しない。
- FRC物理はトカマク物理より成熟度が低い；閉じ込めスケーリング則は十分に確立されていない。
§22. 参考文献
第I部：核物理
- Bosch, H.-S., & Hale, G.M. (1992). "Improved formulas for fusion cross-sections and thermal reactivities." Nuclear Fusion, 32(4), 611–631.
- Atzeni, S., & Meyer-ter-Vehn, J. (2004). The Physics of Inertial Fusion. Oxford University Press.
- Nevins, W.M., & Swain, R. (2000). "The thermonuclear fusion rate coefficient for p-¹¹B reactions." Nuclear Fusion, 40(4), 865.
- Rider, T.H. (1997). "Fundamental limitations on plasma fusion systems not in thermodynamic equilibrium." Physics of Plasmas, 4(4), 1039–1046.
- Putvinski, S.V., et al. (2019). "Fusion reactivity of the pB11 plasma revisited." Nuclear Fusion, 59(7), 076018.
第II部：閉じ込め物理
- Freidberg, J.P. (2014). Ideal MHD. Cambridge University Press.
- Wesson, J. (2011). Tokamaks. 4th ed. Oxford University Press.
- Grad, H., & Rubin, H. (1958). "Hydromagnetic equilibria and force-free fields." Proc. 2nd UN Conf. Peaceful Uses of Atomic Energy, 31, 190.
- Shafranov, V.D. (1966). "Plasma equilibrium in a magnetic field." Reviews of Plasma Physics, 2, 103.
- Troyon, F., et al. (1984). "MHD-limits to plasma confinement." Plasma Physics and Controlled Fusion, 26(1A), 209.
- ITER Physics Basis (1999). Nuclear Fusion, 39(12), 2137–2174. (IPB98スケーリング)
- Tuszewski, M. (1988). "Field reversed configurations." Nuclear Fusion, 28(11), 2033.
- Binderbauer, M.W., et al. (2015). "A high performance field-reversed configuration." Physics of Plasmas, 22(5), 056110. (TAE C-2W)
シリーズナビゲーション
←（開始）| Vol.1: 核物理と閉じ込め | [Vol.2: 点火、燃焼とパワーバランス →]本記事はMITライセンスで公開されています。自由に複製、拡張、批判してください。
すべてのコードは再現可能です。すべての主張は検証可能です。真実は誰のものでもない。