公開したブラウン運動シミュレータでは、無数の粒子が1次元上をランダムに揺れ動きながら、少しずつ左右に広がっていく様子を観察できます。
1個1個の粒子の動きを追うと、次にどちらへ動くかはまったく予測できません。ところが粒子の「分布」に目を移すと話が一変します。ヒストグラムの山は時間とともになめらかに広がり、その形は「拡散方程式」という決まった数式にぴったり重なっていきます。個々の動きはランダムなのに、全体の振る舞いは決定論的。この記事ではその仕組みを、高校で習うニュートン力学と確率の基本だけを使って、一歩ずつ組み立ててみます。
1. なぜ粒子はランダムに動くのか — 力学的な出発点
水の中に浮かべた花粉の粒など、流体に浮んだ小さな粒を想像してください。この粒は、周囲を飛び回る水分子から、毎秒何兆回という桁違いの回数で衝突を受けています。以前の気体分子運動論の記事で扱ったのと同じ、F=maで動く分子たちの世界です。
1回1回の衝突はニュートンの運動方程式に従う完全に決定的な現象です。しかし、
- 1秒間に起きる衝突の数がとてつもなく多い
- それぞれの衝突がどの方向から来るかはばらばらで規則性がない
- 1回の衝突が与える運動量の変化は、大きな粒子にとっては極めて小さい
という三つの事情が重なる結果、
1つひとつは決定論的でも、合計はランダムに見える
という状況が生まれます。無数の衝突を1つずつ追いかけるのは事実上不可能なので、短い時間Δtの間に粒子が受ける影響をまとめて、「平均は0だが大きさはランダムな1回のキック」として近似する—これがブラウン運動を「ランダムウォーク」としてモデル化することの力学的な根拠です。
2. 1ステップの変位を確率変数として式にする
シミュレータは、各粒子の位置を次の更新式で進めています。
\[ X_{n+1} = X_n + \sqrt{2D\Delta t}\, Z_n, \qquad Z_n \sim N(0, 1) \]
\(Z_n\) は平均0・分散1の標準正規乱数です。つまり1ステップあたりの変位 \(\Delta X = X_{n+1}-X_n\) は、
- 平均0（左右どちらに動く確率も等しく、ドリフトしない）
- 分散 \(2D\Delta t\)（\(\sqrt{2D\Delta t}\) を2乗すればそのまま出てきます）
という確率分布に従います。ここで \(D\)（拡散係数）は、流体が粒子をどれだけ激しく揺さぶるかを表すパラメータです。シミュレータで \(D\) を大きくすると、同じ \(\Delta t\) でも1ステップの典型的な移動量が大きくなり、粒子群がより速く広がっていきます。
なお、このシミュレータでは \(D\) を自由に設定できるパラメータとして扱っています。実際の物理では、\(D\) は温度と流体の粘性・粒子の大きさから \(D = kT / (6\pi\eta r)\)（アインシュタインの関係式）として決まりますが、ここでは「拡散のしやすさ」そのものを直接動かして、その効果を確かめる形になっています。
3. ステップを重ねると分散は時間に比例する
1ステップごとの変位 \(\Delta X_1, \Delta X_2, \dots, \Delta X_N\) は、それぞれ独立で平均０・分散 \(2D\Delta t\) の確率変数です。確率論には、
独立な確率変数の和の分散は、各々の分散の和に等しい
という基本法則があります。これをNステップ後の位置 \(X_N = X_0 + \sum_{n=1}^{N}\Delta X_n\) に適用すると、
\[ \mathrm{Var}(X_N) = N \cdot 2D\Delta t = 2D \cdot (N\Delta t) = 2Dt \]
となります（\(t = N\Delta t\) は経過時間）。これが、シミュレータの統計パネルに表示される平均二乗変位（MSD）の理論値、
\[ \mathrm{MSD}(t) = 2Dt \]
そのものです。ポイントは、平均位置は0のまま（ドリフトしない）なのに、広がり方を表す分散は時間に比例して増え続けることです。各ステップはNという回数の足し算だけですが、その単純な加法性が、拡散という現象のスピードを決めています。
さらに、中心極限定理によれば、多数の独立な小さな確率変数の和は、元々の1ステップの分布の形によらず、Nが大きければ正規分布に近づきます。シミュレータでは各ステップのキックをもともと正規乱数（Box-Muller法）で作っていますが、仮に各キックが正規分布でなくても、積み重ねた先の分布は同じです—これが次の節で出てくるガウス分布の伏線となります。
4. 個々はランダムでも、分布は決定論的な方程式に従う — 拡散方程式
ここがこの話の中心です。一番単純なモデル、左右に均等の確率で一歩進む格子上のランダムウォークで考えてみます。位置\(x\)にいる確率密度\(p(x,t)\)は、一つ前の時刻で左隣\((x-\Delta x)\)または右隣\((x+\Delta x)\)にいた粒子が確率\(\tfrac{1}{2}\)ずつで移ってきた結果なので、
\[ p(x, t+\Delta t) = \tfrac{1}{2}\,p(x-\Delta x, t) + \tfrac{1}{2}\,p(x+\Delta x, t) \]
と書けます。右辺を\(x\)周りでテイラー展開して両者を平均すると、一次の項（\(\partial p/\partial x\)）は打ち消し合い、二次の項だけが残ります。
\[ \tfrac{1}{2}\,p(x-\Delta x,t) + \tfrac{1}{2}\,p(x+\Delta x,t) \approx p(x,t) + \tfrac{1}{2}\Delta x^2\,\frac{\partial^2 p}{\partial x^2} \]
一方、左辺は\(t\)方向にテイラー展開すれば \(p(x,t) + \Delta t\,\partial p/\partial t\) です。両辺を等しくおくと、
\[ \frac{\partial p}{\partial t} = \frac{\Delta x^2}{2\Delta t}\,\frac{\partial^2 p}{\partial x^2} \]
ここで \(D = \Delta x^2 / (2\Delta t)\) と置き換えれば（格子間隔を細かくしながらこの比を一定に保つ極限を取る）、シミュレータが実際に解いている
\[\dfrac{\partial p}{\partial t} = D\,\dfrac{\partial^2 p}{\partial x^2}\]
という拡散方程式が得られます。個々の粒子の動きは残さずランダムなまま、その確率分布の時間発展だけが、熱伝導方程式とも呼ばれるこの滑らかな偏微分方程式に属する—これが、カオスなミクロから法則正しいマクロが生まれる仕組みです。
5. 拡散方程式の解 — 2つの道のりが一致する
原点を初期位置とする拡散方程式の解は、
\[ p(x, t) = \frac{1}{\sqrt{4\pi D t}}\,\exp\!\left(-\frac{x^2}{4Dt}\right) \]
です。これは、平均０・分散\(\sigma^2\)の正規分布の密度関数 \(\frac{1}{\sqrt{2\pi\sigma^2}}\exp(-x^2/2\sigma^2)\) に、第3節で導いた \(\sigma^2 = 2Dt\) を代入したものとちょうど一致します。つまり、
- 「多数の独立なランダムステップを積み重ねる」（確率の見方）
- 「格子ランダムウォークの連続極限を取る」（微分方程式の見方）
という二つの全く違う道のり方が、同じガウス分布で合流するわけです。シミュレータの「分布」タブでは、実際の粒子のヒストグラムと、この理論式から描いた曲線を重ねて表示しています。両者がよく重なるのを確かめてみてください。
まとめ
| 力学・確率の要素 | 対応する拡散の量 |
| 無数の分子衝突による瞬間的なランダムな力 | 1ステップの変位 \(\Delta X \sim N(0, 2D\Delta t)\) |
| 独立な確率変数の分散の加法性（確率論の基本法則） | 分散が時間に比例：平均二乗変位 MSD(t) = 2Dt |
| 格子ランダムウォークの連続極限（テイラー展開） | 拡散方程式 \(\partial p/\partial t = D\,\partial^2 p/\partial x^2\) |
| 中心極限定理（多数の小さなランダムな影響の積み重ね） | 確率密度がガウス分布 \(p(x,t)=\frac{1}{\sqrt{4\pi Dt}}\exp(-x^2/4Dt)\) に収束 |
ブラウン運動の面白さは、「個々のランダムさ」と「全体の滑らかな法則」が矛盾しないところにあります。ニュートン力学と確率の基本法則だけから、熱伝導方程式と同じ形の方程式が導かれるというのは、少し不思議ですが確かに成り立つ事実です。
この記事で導いた式を手に、ぜひシミュレータで拡散係数や粒子数を変えながら、ヒストグラムの形やMSDグラフの伸びを理論曲線と見比べてみてください。
最後まで読んでいただきありがとうございます。ご質問やコメントはお問い合わせからよろしくおねがいします。