モンテカルロ放射伝達 MCRT シミュレーター — 光子追跡

吸収・散乱媒質の中を進む光子をモンテカルロ法でランダムウォーク追跡し、拡散反射率・透過率・浸透深さをリアルタイムに評価するツールです。生体組織光学・光線力学療法・薄膜光学・星間ダストなど、散乱優位の放射輸送に共通する MCML 流の枠組みをブラウザ上で体感できます。

パラメータ設定

吸収係数 μ_a

1/cm

媒質1cmあたりに吸収される確率の指標。組織で0.01〜1程度

散乱係数 μ_s

1/cm

単位長あたりの散乱回数。組織では10〜500/cm

異方性 g

Henyey-Greenstein位相関数。g=0等方、g→1前方散乱

媒質厚 d

光子数 log₁₀(N)

追跡する光子数。N=10^値。統計誤差は1/√N

屈折率 n

媒質の屈折率。組織で1.37〜1.45

波長 λ

入射光の波長。可視〜近赤外を想定

計算結果

—

平均自由行程 ℓ (cm)

—

輸送平均自由行程 ℓ* (cm)

—

光学的深さ τ'

—

拡散反射率 R (%)

—

拡散透過率 T (%)

—

浸透深さ δ_p (mm)

—

光子ランダムウォーク — MCRT アニメーション

スラブ媒質に上から入射した光子（白）が散乱（青）・吸収（赤）・反射（緑）・透過（橙）の各イベントを経て終端します。色は終端モードを示します。

反射率 R / 透過率 T vs 光学的深さ τ'

浸透深さ δ_p vs 吸収/散乱比 μ_a/μ_s

理論・主要公式

$$\mu_t = \mu_a + \mu_s,\qquad a = \frac{\mu_s}{\mu_t},\qquad \mu_s' = (1-g)\,\mu_s,\qquad \delta_p = \frac{1}{\sqrt{3\,\mu_a\,(\mu_a + \mu_s')}}$$

μ_a：吸収係数、μ_s：散乱係数、g：異方性（Henyey-Greenstein 位相関数）、μ_s'：減衰散乱係数、δ_p：拡散浸透深さ。τ' = (μ_a + μ_s')·d は減衰光学的深さ。

$$p_{HG}(\cos\theta) = \frac{1}{4\pi}\,\frac{1-g^2}{(1+g^2 - 2g\cos\theta)^{3/2}},\qquad s = -\frac{\ln(\xi)}{\mu_t}$$

Henyey-Greenstein 位相関数と自由行程 s のサンプリング式。ξ ∈ (0,1] は一様乱数。

モンテカルロ放射伝達 — 散乱・吸収媒質中の光子追跡

🙋

「モンテカルロ放射伝達」って、レーザーを当てたときに光がどう散らばるかを計算する手法ですか？聞いたことはあるんですが、なぜわざわざ乱数で解くんですか？

🎓

そう、いいところを突いてる。光が散乱と吸収を繰り返す媒質の輸送方程式（RTE）は7次元の積分微分方程式で、解析的にはほぼ解けないんだ。だから「個々の光子を1個ずつサイコロで追いかける」というのが MCRT の発想。1個の光子で起こることは単純で、(1)自由行程を指数分布からサンプル、(2)散乱角を位相関数からサンプル、(3)吸収か散乱かをアルベドで決める、これだけ。何百万個も流せば、結果としてRTEの解が得られる。

🙋

なるほど、ブルートフォースだけど発想は素直ですね。左の「異方性 g」を 0 から 0.9 にすると、輸送平均自由行程 ℓ* が急に大きくなりました。これって何が起きているんですか？

🎓

g は「散乱の前方バイアス」を表すパラメータで、g=0 が完全等方散乱、g→1 が真前方散乱だ。生体組織では g≈0.9 で、1回の散乱では平均しても 25° くらいしか方向が変わらない。だから方向の記憶が消えるまでに何回も散乱する必要があって、その合計距離が ℓ* = 1/(μ_a + (1-g)μ_s) になる。g=0 と g=0.9 では ℓ* が約 10 倍違う。「光は深く入るけど、すぐ広がる」のはこの効果のせいなんだ。

🙋

浸透深さ δ_p が 17.4 mm って結構深いですね。皮膚に当てた赤色光がここまで届くんですか？

🎓

デフォルトは μ_a=0.1, μ_s=10 の弱吸収・強散乱を想定した値で、現実の皮膚（μ_a≈0.3〜2/cm）よりやや「光が通りやすい」設定だ。実皮膚だと δ_p は数 mm 程度になる。それでも 670 nm の赤色光は皮下までかなり届くから、低出力レーザー治療（LLLT）や光線力学療法（PDT）に使われる。波長を変えると μ_a・μ_s が大きく変わるから、「治療したい深さ」に合わせて波長を選ぶのが現場のコツだよ。

🙋

光子数 N を上げると統計誤差が下がりますよね。でもなぜ MCRT は遅いと言われるんですか？

🎓

統計誤差は 1/√N でしか下がらない、ここがネックなんだ。誤差を 1/10 にするには光子を 100 倍流す必要がある。標準的な MCML は 1 光子あたり数十マイクロ秒で、N=10⁶ で数十秒〜数分。3D 不均質モデルだとさらに重くなる。だから (1) GPU 並列、(2) 分散低減（Russian roulette や photon splitting）、(3) 重要度サンプリング、といった高速化技術が研究されている。最近は CUDA-MCML や MCX のような GPU 実装で、1秒に数億光子流せるようになったよ。

🙋

生体以外でも MCRT は使われるんですか？

🎓

むしろ MCRT は分野横断のメタ手法だよ。星間ダストの放射輻射（DUSTY/HYPERION）、CFD 燃焼解析の放射熱伝達、太陽光発電セルの光取り込み、薄膜光学・OLED 設計、PET/SPECT の校正、海洋光学（クロロフィル分布の遠隔測定）など、「散乱と吸収が同時に起きるあらゆる輸送問題」に適用されている。中性子輸送の MCNP も同じ枠組みで、こちらは原子炉設計の標準ツール。物理過程が変わるだけで、サンプリングの構造は全く同じなんだ。

よくある質問

MCRT は放射輸送方程式（RTE）を解析的に解く代わりに、個々の光子を確率モデルで一つずつ追跡するシミュレーション手法です。光子が次に散乱・吸収されるまでの距離は平均自由行程 1/(μ_a + μ_s) の指数分布、散乱角は Henyey-Greenstein 位相関数（非等方性 g）、吸収か散乱かはアルベド a = μ_s/(μ_a + μ_s) でサンプリングします。多数の光子を追跡してアンサンブル平均をとると、拡散反射率・透過率・吸収分布が得られます。Wang・Jacques・Zheng の MCML が生体組織光学の事実上の標準コードです。

平均自由行程 ℓ = 1/(μ_a + μ_s) は次の散乱・吸収までに光子が直進する平均距離です。一方、輸送平均自由行程 ℓ* = 1/(μ_a + μ_s') は、初期方向の記憶が失われるまでに必要な距離で、減衰した散乱係数 μ_s' = (1-g)μ_s を使います。前方散乱が強い（g→1）と、1回の散乱ではほとんど方向が変わらないため ℓ* ≫ ℓ となり、媒質は事実上「もっと厚く」見えます。生体組織では g≈0.9 のため、ℓ* は ℓ の約10倍です。拡散近似はこの ℓ* に基づきます。

拡散浸透深さ δ_p = 1/√(3 μ_a (μ_a + μ_s')) は、散乱媒質の中で拡散光フルエンスが 1/e に減衰する深さです。半無限の散乱媒質では、深さ z における内部光強度はおおよそ exp(-z/δ_p) に従います。レーザー治療や光線力学療法（PDT）では、薬剤が活性化される深さを決める重要な指標で、波長を選ぶことで δ_p を制御します。例えば 670 nm の赤色光で皮膚の δ_p は数mmオーダーです。

モンテカルロ法では N 個の光子を独立に追跡し、その平均を推定値とします。中心極限定理により、推定値の標準誤差は 1/√N で減少します。つまり N=10⁴ なら誤差 1%、N=10⁶ なら 0.1%。誤差を 10 分の 1 にするには光子数を 100 倍にする必要があり、これが MCRT の最大の弱点です。実務では分散低減法（Russian roulette、photon splitting、importance sampling）や GPU 並列化で計算時間を短縮します。

実世界での応用

生体組織光学・光線力学療法（PDT）：皮膚癌や食道癌などに対する PDT では、光感受性薬剤（PpIX、ベルテポルフィン等）を集積させた組織にレーザーを照射し、活性酸素で病変部を破壊します。治療効果は照射光が「どこまで届くか」で決まるため、MCML 系の MCRT で δ_p を事前計算し、波長（緑色 532 nm vs 赤色 630 nm vs 近赤外 800 nm）と出力を決定します。皮膚科では Dermatology Photonics、眼科 AMD 治療でも標準的なシミュレーション手法です。

核医学イメージング（PET/SPECT/光学トモグラフィ）：PET/SPECT の検出器応答や、近赤外光トモグラフィ（DOT、fNIRS）の画像再構成は、MCRT の出力（感度行列）を順問題として使います。脳血流イメージング fNIRS では、頭皮〜頭蓋骨〜脳の多層 MCRT で「どの深さの信号を拾っているか」を可視化し、ノイズ源（頭皮血流）と本信号（脳活動）を分離します。

CFD 燃焼放射熱伝達と天体物理：ガスタービン燃焼器・ボイラ・ロケットエンジンの炉内放射熱は、煤粒子と気体分子の散乱・吸収を含むため、CFD と結合した MCRT で計算します。ANSYS Fluent の DTRM/Discrete Ordinates と並び、Photon Monte Carlo が高精度解として使われます。天体物理ではダスト円盤や原始星周辺の放射輻射輸送（DUSTY/HYPERION/RADMC-3D）の主力手法で、観測スペクトルとモデルを直接比較できます。

太陽光発電・薄膜光学・OLED 設計：テクスチャ加工された Si 太陽電池、ペロブスカイト薄膜、OLED 多層構造での光取り込み・取り出し効率（吸収率・量子収率）は、波動光学 FDTD と幾何光学 MCRT を使い分けて評価します。リフレクタやテクスチャ層が大きい場合は MCRT が有効で、Sunpower や Tandem セルの設計最適化に使われています。

よくある誤解と注意点

第一の落とし穴は「光子数 N を十分大きくすれば誤差は消える」と思い込むことです。MCRT の統計誤差は 1/√N でしか下がりません。本ツールで N=10⁵ なら誤差約 0.32%ですが、これは「グローバル平均」の誤差で、深さ方向の吸収分布（ローカル量）の誤差はもっと大きくなります。特に深部や弱吸収領域では光子が到達しにくく、ローカル誤差が 10% を超えることがよくあります。深部分布まで精度を出すには分散低減法（Russian roulette、photon splitting、bias scattering）が必須です。

第二に、「拡散近似（diffusion approximation）」と「MCRT」は別物であることに注意してください。本ツールの拡散反射率・透過率の解析式は、δ_p や半無限スラブの近似式で、散乱が吸収より十分大きい（μ_s' ≫ μ_a）かつ媒質厚 ≫ ℓ* の条件で初めて成立します。表面から ℓ* 以内の浅い領域、強吸収領域（皮膚のメラニン層、血管）、媒質界面では拡散近似は大きく外れます。研究レベルの設計では必ずフル MCML を回し、拡散近似は概算用に留めます。

最後に、「Henyey-Greenstein 位相関数は完全ではない」こと。HG 関数は実際の Mie 散乱を 1 パラメータ g で近似した便利な式ですが、後方散乱や 90° 付近の振る舞いは実測と数十%ずれます。皮膚のように後方散乱が重要な系では、Modified HG（Reynolds-McCormick）や Mie 計算による位相関数テーブルを使う方が正確です。生体光学では g 単独だけでなく、g₁・g₂ の 2 項モデルや、波長依存性 g(λ) を別途測定して入れるのが学術標準です。

モンテカルロ放射伝達 MCRT シミュレーター — 光子追跡

モンテカルロ放射伝達 — 散乱・吸収媒質中の光子追跡

よくある質問

実世界での応用

よくある誤解と注意点

使い方ガイド

具体的な計算例

実務での注意点