マクスウェル・ボルツマン分布

古典統計力学で頻繁に出てくるマクスウェル・ボルツマン分布を求めます。
粒子数 N 、体積 V 、エネルギー E を与えて状態数を求めるので、ここでの話はミクロカノニカルです。
体積 V の箱の中に N 個の自由粒子が入っているとし、粒子間の相互作用は十分小さいとし全体のエネルギー Eは各粒子のエネルギーの和になっているとします。この場合の状態数は、各粒子が合計で E になるようにエネルギーを持っている組み合わせの数なので、それを求めます。まず、6 次元位相空間を考えて、微小な体積 a で分割します。この分割した領域に番号をつけて、微小体積 a1 にはエネルギー ϵ1 、微小体積 a2 にはエネルギー ϵ2 、…
が対応するとします。これによって、エネルギー ϵi を持った粒子を、a1 には n1 個、a2 には n2 個、…というように分配できます (N = n1 + n2 + · · · , E = ϵ1 n1 + ϵ2 n2 + · · · )。この分配の仕方の数が可能な状態の数になります。このとき、粒子が微小体積に入る確率は全て同じだとします (エルゴード仮説)。というわけで、状態数を求めるには N 個の粒子を各微小体積に分配する場合の数がどうなっているのかを考えればよく、これは単純な順列の問題なので (各粒子は区別できるとして)w=

N!
n1 !n2 ! · · ·となります。注意ですが、正確にはこれはまだ状態数ではないです。w は微小体積 ai への粒子の分配の仕方なので、ai の中にいくつ状態があるのか定義しないといけないです。例えば a1 の中に m1 個の状態があるとしたら、m1 個の状態に n1 個の粒子を分配するので、その仕方は mn1 1 個あります。なので、

w′ =

N!
mn1 mn2 · · ·n1 !n2 ! · · · 1 2となります。状態数 W にするにはさらに、ni の分配の仕方の全パターンを足す必要があるのでW =∑w′ (ni )niとなります。ni の和は ni の分配の仕方の全パターンの和を取れという意味です。
しかし、今の設定では W を求めるのは難しいので、w の段階で考え、w が一番大きくなる ni を求めます。熱平衡状態ではこの場合 (可能な組み合わせの数が一番多い場合) が最も実現しやすくなっています。なんでかは、状態数 W によるエントロピーの式 (k はボルツマン定数)

S = k log W′が近似的に、w′ が最大のときの wmaxによって′S = k log W ≃ k log wmaxとできるからです。実際に差を計算してみると、log nmax程度の差が出るだけなので (nmaxは w′ を最大にする
ii′ni )、十分近似として成立しています。で、w でなく w で行っても同じ結果になるので、w でやっていくということです。
計算をしやすくするために w の対数をとってlog w = log N ! − log n1 ! − log n2 ! − · · ·∑= log N ! −log ni !
i=1

1

N が十分大きいとしてスターリングの公式
log A! ≃ A(log A − 1)
を使うとlog w = N (log N − 1) −∑ni (log ni − 1)i=1

= N log N − N −∑ni log ni +i=1

= N log N −

∑

∑

nii=1

ni log nii=1

W が最大となる ni が知りたいので、極値を求めます。これを ni で微分したのが 0 になればいいので、各 ni で微分すれば−∑(log ni + 1) = 0(1)i=1

これをもう一度 ni で微分したものは負なので、極大値になっています。ただし、今は∑ni = N(2)i=1と、系の全エネルギー E がエネルギー ϵi の和になっている必要があるために∑ϵi ni = E(3)i=1という制限がかかります。この 2 つの条件を満たすように極値を選ぶ必要があります。このように拘束条件があるときの極値はラグランジュの未定乗数法 (解析力学の「拘束条件」参照) を使えばいいです。
拘束条件を ϕ(xi ) としたとき、関数 f (xi ) の極値は

∂ϕ∂f+λ=0∂xi∂xiで求められます。今の場合では、∂f /∂xi は (1) の左辺で、拘束条件は (2) と (3) なので、未定乗数を α, β として∑∑ ∂ni∑ ∂ni∑(log ni + 1) + α+βϵi=(log ni + 1 + α + βϵi )∂ni
∂nii=1i=1i=1i=1
これが 0 になればいいのでlog ni + 1 + α + βϵi = 0ni = exp[−1 − α − βϵi ]−1 − α は定数なので、この部分を分離することで ni は2

ni = C exp[−βϵi ](4)となります。明確に書きませんが、ここでの ni は w を最大にする場合です。C は (2) からC∑

e−βϵi = Ni=1に従う必要があります。C = N/f とすればN ∑ −βϵie=N
f i=1
f=

∑e−βϵii=1

このように見やすくなります。また、E はE=C∑ϵi e−βϵi =i=1N ∑ −βϵiϵi ef i=1

未定乗数として導入した β がなんなのか決めないと意味のある量を計算できないので、決定させます。そのためにエネルギーの平均を計算してみます。
分割した体積 a とは違う大きさとして、d3 xd3 p という領域を考えます (d3 x = dxdydz, d3 p = dpx dpy dpz )。ai内にいる粒子数 ni は (4) で与えられていて、d3 xd3 p にいるエネルギー ϵi を持った粒子数は ni d3 xd3 p/a です (a は6 次元位相空間の体積)。全粒子数 N は ni の和ですが、d3 xd3 p にいる密度 ni d3 xd3 p/a から、ni が i でなく x, pに依存していると考えれば、和は位相空間の全領域にわたる積分だと見なせて∫∫
3

d xd3 pn=Naと書けます。そうすると、d3 xd3 p に粒子がいる確率は

1 n 3 3nd3 xd3 pd xd p = ∫ 3 ∫ 3Nad x d pnとなります。これは粒子がエネルギー ϵ(i 番目のエネルギー ϵi を x, p に依存するとした ϵ) を持つ確率と言う事ができます (ni はエネルギー ϵi を持った粒子の個数)。なので、これにエネルギー ϵ をつけて位相空間で積分することで平均エネルギー < ϵ > が∫ 3 ∫ 3d x d p ϵe−βϵ< ϵ >= ∫ 3 ∫ 3 −βϵd x d peと求まります。今は自由粒子を考えているので、エネルギー ϵ を自由粒子でのϵ=|p|21 2=(p + p2y + p2z )2m2m xとすることで、x 積分はどこにも引っかからずにただの空間体積 V になって、打ち消しあい3∫ 3 −βϵd pϵe< ϵ >= ∫ 3 −βϵd pe
よって、後は∫d3 p|p|2 −β|p|2 /2me,2m∫

d3 pe−β|p| /2m2

これらの積分を実行すればいいです。ただ、式変形をすると∫ 3 −βϵ∫d pϵe∂< ϵ >= ∫ 3 −βϵ = −log[ d3 pe−βϵ ]∂βd peとできるので、ϵ が外にいない方を考えるだけでいいです。この積分は∫d3 pe−β|p| /2m =2∫ ∞dpx dpy dpz e−βpx /2m e−βpy /2m e−βpz /2m22dpx dpy dpz e−βpx /2m e−βpy /2m e−βpz /2m = (2mπ 3/2)β−∞2なので、個別にガウス積分を行うことで∫ ∞−∞222となります。よって平均エネルギーは< ϵ >= −∂2mπ 3/2log()∂ββ=

3 β 2mπ2 2mπ β 2=32β

ここで気体運動論での理想気体 (自由粒子の集まり) のベルヌーイの関係（P は圧力）PV =2N <ϵ>3を持ち出すとPV =Nβそして、理想気体の状態方程式 (n0 は気体のモル数、R は気体定数、T は温度)

P V = n0 RT
から

4N
βn0 RT =Nn0 RT

β=n0 を 1 モルとすれば、粒子の総数 N はアボガドロ数 NA になるので

β=NA 1R T

ここで k = R/NA とすれば

β=1kTとなります。k はボルツマン定数の定義そのものです。下の補足で別の導き方を載せています。というわけで、

ni =N −βϵie(β = 1/kT )
f
∑
e−βϵi
f=i=1

これをマクスウェル・ボルツマン分布といいます。
ここまでは空間座標と運動量で考えてきましたが、運動量を速度に書き換えてみます。運動量と速度は (太字は三次元ベクトルです)

p = mv
なので、d3 p = m3 d3 v となります。このことから、結局比例関係になっていることが分かります。つまり、d3 xd3 pにいる粒子数と同じように d3 xd3 v にいる粒子数も e−βϵ に比例した形になります。このことから、d3 xd3 v にいる粒子数を

f (v)d3 xd3 v = Ae−βϵ d3 xd3 vと書くことにし、これを全領域で積分すれば全粒子数 N になるので∫
f (v)d3 xd3 v = Nとします。速度に書き換えるときに出てくる比例係数 A を e−βϵ に含めさせたものを f (v) と定義しているだけです。
自由粒子だとして A を求めます。自由粒子では∫∫
33f (v)d xd v =Ae−mβ|v| /2 d3 xd3 v = N2x 積分はここでも体積 V になるだけです。v 積分はまたガウス積分を使うことで5∫ ∞−∞

Ae−mβ(vx +vy +vz )/2 dvx dvy dvz = (2222π 3/2)mβ

よって

V A(2π 3/2)
=NmβA=

N 2π −3/2()
V mβとなるので

f (v)d3 xd3 v =N βm 3/2βm|v|2 3 3() exp[−]d xd v
V 2π2(5)これをマクスウェルの速度分布 (velocity disribution) と呼びます。マクスウェル・ボルツマンの分布から求めましたが、歴史的にはマクスウェルの速度分布の方が先に導出されています。速度とついているように空間座標による V と d3 x を外している場合が多いです (速度空間だけで行えば最初から空間座標は出てこない)。また、確率に書き換えるなら全粒子数 N で割ってn1 βm 3/2βm|v|2 3 3= () exp[−]d xd vNV 2π2とすればいいです。
マクスウェルの速度分布で勘違いしやすいことを注意しておきます。f (v) は 3 次元の速度 v の関数ですが、右辺を見てみると明らかに成分ごとに分けることができて、例えば x 成分の f (vx ) はf (vx ) =N βm 1/2βmvx2() exp[−]V 2π2のように書けます。他の成分も同様になっているので、状況を簡単に考えるために x 成分のみを考えるとします。
そうすると、f (vx ) は exp の形から vx2 = 0 のときを最大値とし、|vx | の増加にともなって減少していることがすぐに分かります。このことを f (v)d3 xd3 v にそのまま適用させてみると、f (v)d3 xd3 v は d3 xd3 v にいる粒子数なので、速度が 0 の粒子が一番多くいるように思えてしまいます。しかし、vx = 0 はあくまで速度の x 成分のことをいっているのにすぎなく、他の成分 vy , vz が速度 0 になっているとは言っていません。
実際に、全成分があるときの粒子数を最大にする速度を求めてみます。このとき速度の方向には特別な向きがないとします (等方性)。言い換えれば、どの方向を向いているかは無関係だとするということです。そうすると、各成分ごとに分けて考える必要がなくなります。なので、|v| ∼ |v| + d|v| の間の粒子数について考えます。
|v| ∼ |v| + d|v| の間の体積は、3 次元の微小体積|v|2 sin θdϕdθd|v|の角度積分を行ったもの (半径 |v| の球の厚さ d|v| の球殻の体積) なので4π|v|2 d|v|
この領域に入っている粒子数を求めればいいです。粒子数は (5) の d3 v をこれに置き換え、今は空間座標の位置については特定しないので体積 V と d3 x をはずした (もしくは x 積分を実行した)

n = 4πN (βm|v|2βm 3/2) exp[−]|v|2 d|v|
2π26(6)

これが |v| ∼ |v| + d|v| 間にいる粒子数です。これを見てみると、速度が 0 となる |v| = 0 では n = 0 になっています。なので、実際に速度が 0 のときが粒子数が多いということを言っていないのが分かります。粒子数 n が最大値となる速度は極値の計算

dnβm 3/2βm|v|2βm 3/2βm|v|2= 8πN () |v| exp[−]d|v| − 4πN () βm|v|3 exp[−]d|v|d|v|
2π22π2= 4πN (βm 3/2βm|v|2) (2|v| − βm|v|3 ) exp[−]d|v|
2π2=0から√|v| =2βmとなります。
ついでにマクスウェルの速度分布での平均速度 < |v| > も求めておきます。速度の絶対値の平均であることに注意してください。速度ベクトル v の平均はマクスウェルの速度分布の分布の仕方が等方的 (イメージ的には 3 次元球の内部で放射状に速度ベクトルが散っている) であるために < v >= 0 (< vx >=< vy >=< vz >= 0) です。
平均は (6) に |v| をかけて積分して N で割ればいいので< |v| >= 4π(βm 3/2)2π

∫ ∞exp[−0βm|v|2]|v|3 d|v|2を計算すればいいです。βm|v|2 /2 = x として∫ ∞exp[−02βm|v|2]|v|3 d|v| =2(βm)2∫ ∞x exp[−x]dx0(dx = βm|v|d|v|)
この積分はガンマ関数 Γ(n) によって∫ ∞xn e−ax dx =0

Γ(n + 1)
an+1となることから2(βm)2∫ ∞x exp[−x]dx =022
Γ(2) =2(βm)
(βm)2

よって

βm 3/2 22< |v| >= 4π()=√22π(βm)πとなります。
また、平均エネルギー < E >=< m|v|2 /2 >= 3kT /2 から求まる

7√2βm

√vt = < |v|2 > =√3kTmのことを熱速度 (thermal valocity) と言います。定義によっては 3 を省いている場合もあります。
・補足
平均エネルギーの比較で β = 1/kT としましたが他の方法もあるので、それを示します。熱力学で熱量 Q とエントロピー S の関係は、系に熱量 δQ を加えた時のエントロピーの変化として

δQ = T dS(7)と与えられます。これは系に熱量を加えることで各状態の粒子数が変化したと見ることができます。そして、log wは大雑把には状態数に対応するので、log w の ni による変化を求めれば対応する関係が出てくると予想できます。
log w の変化は上で求めたように

δ log w = −∑(log ni + 1)δnii=1ni として
ni = C exp[−βϵi ]を使えば

δ log w = −∑(log C − βϵi + 1)δnii=1

そして、全粒子数は変化しないために∑δni = 0i=1であることと、ϵi δni は系のエネルギー E の変化分なので∑ϵi δni = δEi=1と書けることを使えば

δ log w = βδEとなります。
系の変化は熱量によるものとしかしていないので、系のエネルギーの変化は熱量の変化 δQ と同じことから

δ log w = βδE = βδQ
そして、w を状態数だとみなしてしまえば δS = kδ log w なのでδS = kδ log w = kβδQ8よって (7) から

β=1kTとなります。
もっと手っ取り早く温度に対応することをみることができます。1 つの系で考えてきましたが、同じ系をもう 1つ用意して、その間でエネルギーの交換が起きているとします。2 つの系を A, B とすれば、くっつけた系ではw=N A!

NB!

AB B
nA
1 !n2 ! · · · n1 !n2 ! · · ·B
これを最大にする nAi , ni はnA
i =N A −βϵA
N B −βϵBi , nB =ee iifAfB

このようになります。1 つだけの場合に同じものがくっついてくるだけなので、こうなることは分かると思います。ただし、拘束条件は∑i=1nAi = NA ,∑nB
i = NB ,i=1

∑
i=1A
ϵA
i ni +∑B
ϵB
i ni = Ei=1

このようになっていて、E による拘束条件での未定乗数 β が、A, B 両方に対して同じになっています。そして、
2 つの系で同じ β が使われているということと、熱平衡では温度が等しいということを合わせれば、β は温度に対応することが分かります。

9