リウヴィル方程式

ミクロカノニカル、カノニカル、グランドカノニカルは時間依存性を持たない分布関数によるものでしたが、ここでは分布関数が時間依存している場合を考え、それに対する基本的な方程式であるリウヴィル方程式を古典的、量子論的な両方の場合で求めます。
演算子での場合でもハットをつけたりしていないので古典論の場合と記号が混ざっていますが、分けて話をしているので多分混乱しないと思います。
最後にリウヴィル方程式を使った線形応答を見ます。

1 つの粒子が正準方程式 (ハミルトン方程式)

∂H dp
∂H
dq=,=−dt∂pdt∂qに従っているとし、N 個の粒子によって形成される系がdqi∂H dpi∂H==−,dt∂pidt∂qi(i = 1, 2, . . . , N )によって記述される場合を考えます。p と q が時間に依存していることは明確に書かないので注意してください(ただし、力学変数である p, q の時間の偏微分は 0 なので陽に依存していない。解析力学の「ハミルトン形式」参照)。力学の問題としてはこの 6N 個の方程式の解を求めることになりますが、これを変更します。N 個の粒子に対する運動方程式を直に解くことから振る舞いを知ろうとするのではなく、確率分布 (分布関数) によって平均的な振る舞いを求めることを考えます。これは統計力学の基本的な方針そのものです (巨視的な現象は微視的な振る舞いの平均値で十分実験結果を再現できると考える)。
確率を導入するために、正準方程式に従って実現する粒子の状態に対する分布を考えます。粒子の状態は p, q で指定されるので、その状態を表す分布関数は時間 t での 6N 次元位相空間において f (p, q, t) と与えます。p, q は正準方程式に従っていて、(p1 , p2 , . . . , pN ) を p、(q1 , q2 , . . . , qN ) を q としています。f (p, q, t) はあらゆる時間に対して全位相空間積分によって

1N !(2πℏ)3N∫∫f (p, q, t)dpdq =f (p, q, t)dΓ = 1と規格化されています。
分布関数がどのように時間依存しているのか見ていきます。リウヴィルの定理からの関係 (「情報エントロピーとの関係」参照)

f (p, q, t) = f (p′ , q ′ , t′ )において、t′ を t から微小時間 dt 経過したものだとすると、p′ , q ′ は p, q から t + dt 経過したものなのでp′ = p +dp
dqdt , q ′ = q + dt
dtdtとなります。そうすると

1f (p, q, t) = f (p +dpdqdt, q + dt, t + dt)
dtdt= f (p, q, t) +∂f (p, q, t) dp∂f (p, q, t) dq∂f (p, q, t)
dt +dt +dt
∂pdt
∂qdt
∂t

= f (p, q, t) +N∑∂f (p, q, t) dpi∂f (p, q, t) dqi∂f (p, q, t)(·+·)dt +dt
∂pdt
∂qdt
∂tiii=1

最後の行は省略している内積がどうなっているのか分かりやすくするために書いています (粒子数を N = 1 にすれば和がなくなるのでもっと分かりやすくなります)。等式を成立させるには第二項と第三項が 0 になっている必要があるのでN∑dpi ∂fdqi ∂f∂f(·+·)+=0dt ∂pi
dt ∂qi∂ti=1(1)q, p は正準方程式に従っている必要があるので、正準方程式を使うことで分布関数の時間微分はN∑∂fdpi ∂fdqi ∂f=(−·−·)∂tdt∂pdt ∂qiii=1=N∑dH ∂f∂H ∂f(·−·)
dqi ∂pi∂pi ∂qii=1= {H, f }となります。この方程式のことをリウヴィル (Liouville) 方程式と言います。最後の { } はポアソン括弧{F, G} =N∑∂F ∂G∂F ∂G(·−·)∂q∂p∂piii ∂qii=1を表します。
リウヴィル方程式は分布関数に対する基本的な方程式になっています。リウヴィル方程式が解ければ各粒子がどういう状態になっているかを表す分布関数が分かり、適当な平均値を求めることができます。後で見るように平衡状態に対して外部から影響を与えたときの反応もリウヴィル方程式から分かります。
また、ここまでの話からなんとなく分かるように、分布関数を特別に選んだ (目的に合うようにした) ものが、ミクロカノニカル、カノニカル、グランドカノニカルです。なので、統計力学はリウヴィル方程式の話から入っていくほうが正統的です (分かりやすいかどうかは別にして)。
ハミルトニアンと分布関数が時間にあらわに依存していないとき (∂f (p, q)/∂t = 0)、リウヴィル方程式はN∑∂H ∂fdH ∂f·−·)=0{H, f } =(dq∂p∂piii ∂qii=1(2)となり、解析力学の話から分布関数が保存量になることを意味します (解析力学の「ハミルトン形式」参照)。つまり、系の運動は分布関数 f が一定になる領域に制限されます (f 一定面を動く)。そして、分布関数が時間依存しないことから、平衡状態の分布関数を表します (位相空間全体にわたって時間独立)。また、f の p, q 微分が 0 だとすれば ∂f /∂t = 0 になるので、位相空間が一様だとすれば平衡状態になります。

2さらに、分布関数がある保存量 S(p, q) (t にあらわに依存しない。∂S/∂t = 0) のみに依存するなら、分布関数はリウヴィル方程式 (2) を満たすと言えます。これは、S が保存量であるために系の運動は S が一定となる領域に制限され、それに伴って系の運動は分布関数が一定 (分布関数が S のみに依存するため) となる領域に制限されるからです。このようにして平衡状態を考えるとき分布関数は保存量のみの関数にすればいいという話になり、ミクロカノニカルアンサンブルが出てきます。
リウヴィル方程式は∂f= iLf , Lf = −i{H, f }∂t(3)という形で書かれることがあり、演算子 L はリウヴィル演算子と呼ばれます。もしくは∂f= Lf , Lf = {H, f }∂tと定義される場合もあります。ここでは (3) の表記を使います。このように書いたとき、リウヴィル方程式の解を形式的にきれいに書くことが出来ます。t = 0 での分布関数を f (p, q, 0) とし、L が時間に依存していないなら解は

f (p, q, t) = eiLt f (p, q, 0)と書けます。解になっていることは、時間微分∂f (p, q, t) = iLeiLt f (p, q, 0) = iLf (p, q, t)∂tから確かめられます。また、関数 A が時間に p, q を通してしか依存していないときのポアソン括弧を使った正準方程式 (p = p(t), q = q(t))

dA(p, q)= −{H, A}dt
をリウヴィル演算子を使って書けば

dA(p, q)= −iLA(p, q)dt
なので、L が時間に依存していないなら、正準方程式に従う A は t = 0 で A(p(t), q(t)) が A(p(0), q(0)) だとすることでA(p(t), q(t)) = e−iLt A(p(0), q(0))
この e−iLt は時間発展演算子と呼ばれます。名前の通り関数 A(p(0), q(0)) を A(p(t), q(t)) に変化させるからです。
関数の時間発展とリウヴィルの定理と合わせることで、eiLt のかかる方向を変えることが出来るのでそれも示しておきます。関数 A(p(0), q(0)) と eiLt B(p(0), q(0)) の全位相空間での積分∫dpdqA(p(0), q(0))eiLt B(p(0), q(0))

I=∫=dpdqA(p(0), q(0))B(p(−t), q(−t))

を考えます。二行目へは時間発展演算子による書き換えを行っています。位相空間全体にわたった積分なので、積分変数 p, q を時間 t だけ動かした3

A : p(0) ⇒ p(t) , q(0) ⇒ q(t)B : p(−t) ⇒ p(0) , q(−t) ⇒ q(0)という書き換えをしても問題なく∫
I=dp(t)dq(t)A(p(t), q(t))B(p(0), q(0))

微小な位相体積はリウヴィルの定理から dp(t)dq(t) = dpdq なので∫
I=dpdqA(p(t), q(t))B(p(0), q(0))

そして、A を時間発展演算子で書くと∫
I=dpdqe−iLt A(p(0), q(0))B(p(0), q(0))となるので∫∫dpdqA(p(0), q(0))eiLt B(p(0), q(0)) =dpdqe−iLt A(p(0), q(0))B(p(0), q(0))(4)という関係があることが分かります。よって、右側の B にかかっていた eiLt を左側の A に作用させるときは、この関係に従います。
(1) が連続の方程式の形になっていることから、リウヴィル方程式は位相空間における連続の方程式だと見ることが出来ます。実際にそうなっていることも示しておきます。まず、3 次元空間における連続の方程式∂ρ(x, t)= −∇ · J∂tを持ち出します。ρ は密度、J はフラックス (カレント) で、流速度 dx/dt = v と密度によって J = vρ となっています。これを位相空間に持っていきます。位相空間は p と q で張られる 6N 次元空間ですが、簡単のために 1 粒子で px , qx のみによる 2 次元空間だとします。そこに 2 次元の四角形を用意し、そこに単位時間あたり px 方向からJ(px )∆qx の量の粒子が入ってくるとします。で、四角内部から J(px + ∆px )∆qx だけ出て行くとします。qx 方向でも同様に考えると、px , qx 方向から入ってくる量と出て行く量の単位時間当たりの差 (四角内部にたまる量) は(J(px + ∆px )∆qx )out − (J(px )∆qx )in ≃∂J(px )∆px ∆qx∂px(J(qx + ∆qx )∆px )out − (J(qx )∆px )in ≃∂J(qx )
∆px ∆qx∂qx

この 2 つの和になります (辺の長さは微小として ∆px , ∆qx としています)。微小な時間間隔 ∆t 経てば∂J(px )
∂J(qx )
∆px ∆qx ∆t +∆px ∆qx ∆t∂px∂qxだけ四角内に溜まります。これにともなって四角内の物質の密度 ρ(px , qx , t) も変化するので微小時間 ∆t の間に密度は4∂ρ(px , qx , t)∆t∂tと変化します。そして、これに面積 ∆px ∆qx をかけることで、物質の流入出のつりあいから位相空間での連続の方程式∂ρ(px , qx , t)∂J(px )
∂J(qx )
∆px ∆qx ∆t = −∆px ∆qx ∆t −∆px ∆qx ∆t
∂t
∂px
∂qx
∂J(px ) ∂J(qx )∂ρ(px , qx , t)
= −−∂t∂px∂qxが出てきます。これは px を x、J(px ) を Jx 、qx を y 、J(qx ) を Jy だと思えば 2 次元の連続の方程式です。J は
ρdx/dt なので、これを今の場合に当てはめればJ(px ) =dqxdpx
ρ , J(qx ) =ρdtdtそして、6N 次元に持っていくとN∂∂ρ ∑ (∂dpidqi )−=ρ) −ρ)·(·(∂t∂pdt∂qdtiii=1で、密度 ρ を確率分布である分布関数 f (p, q, t) だとすることでN∑(∂f∂∂dpidqi )=−f) −f)·(·(∂t∂pidt∂qidti=1

正準方程式を使うと

NN
∑
∑
∂f∂∂H
∂
∂H
=·(f) −·(f)∂t∂p∂q∂q∂piiiii=1
i=1=N∑( ∂Hi=1=∂qiN∑( ∂Hi=1∂qi·∂f∂ ∂H∂H ∂f∂ ∂H )+(·)f −·−(·)f∂pi∂pi ∂qi∂pi ∂qi∂qi ∂pi·∂f∂H ∂f )−·∂pi
∂pi ∂qi= {H, f }となり、連続の方程式を求める方法でリウヴィル方程式が出てきます。また、(1) から、リウヴィル方程式はNdf∂fdp ∂fdq ∂f∂f ∑ dpi ∂fdqi ∂f=++=+·+·)=0(dt∂tdt ∂pdt ∂q∂tdt∂pdt ∂qiii=1となっていることも分かり、これを流体の概念に合わせれば、流線に沿った分布関数は時間変化しないことを表しています。

5

今度は量子論的なリウヴィル方程式を出します。最初に対応が見やすいように c 数に持っていった形を出します。量子論での演算子 A の統計平均は、量子力学での状態 ψi とその状態が現れる確率 Pi によって∑⟨A⟩ =Pi ⟨ψi |A|ψi ⟩iで与えられます (混合状態)。これを位置表示の波動関数にもっていくと∫⟨A⟩ =d3 xd3 x′∑Pi ⟨ψi |x′ ⟩⟨x′ |A|x⟩⟨x|ψi ⟩i∫=d3 xd3 x′∑

Pi ψi∗ (x′ )A(x′ , x)ψ(x)i∫=

d3 xd3 x′ A(x′ , x)ρ(x, x′ )となります。⟨x′ |A|x⟩ を x′ , x を行列成分と見ることで A(x′ , x) は行列です。本質的に c 数なので交換関係を気にする必要はないですが、行列の積を成分で書く場合 (AB = Aij Bjk ) と同じように A(x′ , x)ρ(x, x′ ) のように書かれます。ここでの

ρ(x, x′ ) =∑Pi ψi (x)ψi∗ (x′ )iは位置表示での密度行列です。見て分かるように、密度行列が確率分布を表していることから、古典的な場合での分布関数に対応しています。ちなみに、密度行列は確率 Pi と波動関数の規格化から∫∫
d3 xρ(x, x) =d3 x∑

Pi ψi (x)ψi∗ (x) = 1i∑
(Pi = 1 ,∫d3 xψi (x)ψi∗ (x) = 1)

iと規格化されています。というわけで、密度行列に対する時間微分をしたものが量子論的な場合になります。そうすると波動関数の時間微分が出てくるので、シュレーディンガー方程式iℏ∂ψi (x, t) = Hψi (x, t)∂tを持ち出します。H はハミルトニアン演算子です。ハミルトニアン演算子も同じように位置表示に直すと∫Hψ(x, t) = H⟨x|ψ; t⟩ = ⟨x|H|ψ; t⟩ =3 ′′′d x ⟨x|H|x ⟩⟨x |ψ; t⟩ =∫d3 x′ H(x, x′ )ψ(x′ , t)

位置の状態 ⟨x| にハミルトニアン演算子が作用することで波動関数にかかる運動量演算子 −iℏ∇ が出てきます。これによって∂iℏ ψi (x, t) =∂t∫

d3 x′ H(x, x′ )ψ(x′ , t)

これを踏まえて時間依存させた行列密度 (波動関数が時間依存する) の時間微分をしてみると

6

ρ(x, x′ , t) =∑Pi ψ(x, t)ψi∗ (x′ , t)iiℏ

∑ ( ∂ψi (x, t)∂ρ(x, x′ , t)∂ψ ∗ (x′ , t) )Pi iℏ=
ψi∗ (x′ , t) + ψi (x, t)iℏ i
∂t∂t
∂ti
=

∑Pi(∫

d3 yH(x, y)ψi (y, t)ψi∗ (x′ , t) − ψi (x, t)∫)
d3 yH ∗ (x′ , y)ψi∗ (y, t)

i

∫
=()
d3 y H(x, y)ρ(y, x′ , t) − ρ(x, y, t)H ∗ (x′ , y)

Pi は時刻 t での確率なので時間依存していません。ハミルトニアンはエルミート演算子なので H ∗ (x′ , y) =H(y, x′ ) (位置を行列成分と見て転置させている) としてiℏ∂ρ(x′ , x, t)
=∂t

∫

()
d3 y H(x, y)ρ(y, x′ , t) − ρ(x, y, t)H(y, x′ )これが量子論的なリウヴィル方程式です。右辺が交換関係の形になっているので、演算子の形では交換関係で書けるのではないかと予想できます。実際に密度行列のブラケット表記 (「分配関数と経路積分」参照)

ρ(t) =∑Pi |ψi ; t⟩⟨ψi ; t|iを時間微分するとiℏ∑∑∂ρ∂∂=Pi iℏ |ψi ; t⟩⟨ψi ; t| +Pi |ψi ; t⟩iℏ ⟨ψi ; t|
∂t
∂t
∂tii=

∑Pi H|ψi ; t⟩⟨ψi ; t| −i∑

Pi |ψi ; t⟩⟨ψi ; t|Hi= Hρ − ρHとなっています。これはハミルトニアン演算子 H の式iℏ∂|ψi ; t⟩ = H|ψi ; t⟩∂tを使っているだけです（ハミルトニアン演算子はエルミート）。よって、演算子に対するリウヴィル方程式は交換関係 [A, B] = AB − BA によってiℏ∂ρ= [H, ρ]∂tとなっていることが分かり、ポアソン括弧が交換関係になるという古典論と量子論の対応そのままになっています。
最後にリウヴィル方程式を使った線形応答の話を見ておきます。線形応答は系に与える影響とその応答が線形性を持っていると言うことです。ようは線形の方程式に従っている系のことです。まず、時間独立なハミルトニアンH0 (p, q) を用意します。これに時間依存する項 Ht (p, q) を加えることを考えます (時間依存性を添え字の t で表しています)。これは元々は時間独立なハミルトニアン H0 (p, q) によって記述されている状況に、外部からの影響として時間依存する Ht (p, q) を加えたということです。H0 (p, q) の系に影響 Ht (p, q) を与える系の存在自体による影響はないものとします (単に外部からの影響が加わっただけだとする)。なので、全ハミルトニアンは7

H = H0 (p, q) + Ht (p, q)となります。そして、Ht (p, q) は無限大の過去では(t = −∞)Ht (p, q) = 0となっているとします。これは始めの状態はハミルトニアン H0 (p, q) のみで記述されていることを意味します。
この状況でリウヴィル方程式を考えます。リウヴィル方程式は時間依存部分があっても変更されないので∂f= {H, f } = {H0 + Ht , f }∂t
初期条件として時間無限大の過去で Ht がいない平衡状態になっているというのを与えます。つまり、t = −∞ でf はカノニカルアンサンブルでの分布関数 f0 になっているとします (もしくはグランドカノニカル)。この初期条
件は Ht が十分小さいときに有効だということに注意してください。Ht が大きいと系は強い非平衡状態になるた
めに、t = −∞ で平衡状態だったという設定は無意味になります (過去に平衡状態であったことは強い非平衡状態
ではどうでもよくなる)。さらに、カノニカルアンサンブルとしていことから分かるように熱浴が存在しているんですが、その取り扱いを無視していることにも気を付けてください。
ハミルトニアン H0 でのリウヴィル演算子 L0 と Ht でのリウヴィル演算子 L1 によってリウヴィル方程式は∂f= i(L0 + L1 )f∂tと書けます。f を新しく g(p, q, t) を導入してf (p, q, t) = eitL0 g(p, q, t)と置き換えてみると (L0 が t に作用しないことを分かりやすくするために tL0 としています)、g(p, q, t) に対する式として∂f= i(L0 + L1 )f∂tiL0 eitL0 g + eitL0∂g= i(L0 + L1 )eitL0 g∂t∂g= ie−itL0 L1 eitL0 g∂tL0 は演算子なので L1 と交換させていません。ちなみに {H0 , f0 } = 0 なので eitL0 f0 = f0 であるために、g(t =−∞) = f0 という条件がかかっています。両辺を積分して (p, q の依存性は省いて書きます)∫ tg(t) = i′′dt′ e−it L0 L1 eit L0 g(t′ )積分範囲を −∞ ∼ t として、t = −∞ で g0 = g(t = −∞) になるとすれば8∫ tg(t) = g0 + i′−∞′dt′ e−it L0 L1 eit L0 g(t′ )∫ teitL0 g(t) = eitL0 g0 + ieitL0∫ tf (t) = f0 + i= f0 +−∞−∞′′−∞∫ t′dt′ e−it L0 L1 eit L0 g(t′ )dt′ ei(t−t )L0 L1 f (t′ )′dt′ ei(t−t )L0 {Ht′ , f (t′ )}となります。最後に (3) を使ってポアソン括弧にしています。これは最初の近似として f (t) = f0 を与えて右辺第二項に入れれば∫ tf (t) = f0 + i−∞′dt′ ei(t−t )L0 L1 f0 = f0 +∫ t′−∞

dt′ ei(t−t )L0 {Ht′ , f0 }という 1 次近似の形になり、この段階までを考慮するのが線形応答です。この結果でてきた f (t) をまた第二項に入れるということを繰り返えして出てくる展開式を考慮したのが非線形応答となります (Ht が 1 次じゃなくなるから)。

Ht の影響を受けた状況での A(p, q) の統計平均は 1 次近似までで∫⟨A⟩ =∫dΓA(p, q)f (p, q, t) =∫ tdΓA(p, q)f0 (p, q) +∫=∫′dt−∞′dΓA(p, q)ei(t−t )L0 {Ht′ , f0 (p, q)}−∞∫ tdΓA(p, q)f0 (p, q) +dt′∫′dΓ{Ht′ , f0 (p, q)}e−i(t−t )L0 A(p, q)となります。最後の行へは (4) を使っています。
ここでの話は具体的な量を使っていないので分かりづらいですが、状況を簡単に言えば、平衡状態に時間依存する寄与を与え、その反応を統計平均を通して見るというものです。これを線形性を維持する範囲内で行うのが線形応答理論です。与える寄与としては電場がよく現れます。分布関数ではないですが、かなり簡単な線形応答の例を有限温度の場の理論の「グリーン関数とスペクトル関数」の補足に載せています。
量子論でも同じように求めることができます。状況設定は同じで、出発点が演算子形式でのリウヴィル方程式iℏ∂ρ(t)= [H, ρ] = [H0 + Ht , ρ]∂tになるというだけです。ρ をρ(t) = e−iH0 t/ℏ ρ′ (t)eiH0 t/ℏ = U −1 (t)ρ′ (t)U (t)として入れると左辺は9∂ρi∂ρ′ iH0 t/ℏ
i= − H 0 e−iH0 t/ℏ ρ′ eiH0 t/ℏ + e−iH0 t/ℏ
e+ e−iH0 t/ℏ ρ′ H0 eiH0 t/ℏ∂tℏ∂tℏ= −∂ρ′iiH0 U −1 ρ′ U + U −1
U + U −1 ρ′ H0 Uℏ∂tℏ= −i −1∂ρ′i
U H0 ρ′ U + U −1U + U −1 ρ′ H0 Uℏ∂tℏ

= −i −1∂ρ′U [H0 , ρ′ ]U + U −1Uℏ∂t

H0 と U −1 = e−iH0 t/ℏ は交換します。右辺は[H0 + Ht , U −1 ρ′ U ] = [H0 , U −1 ρ′ U ] + [Ht , U −1 ρ′ U ]= U −1 [H0 , ρ′ ]U + Ht U −1 ρ′ U − U −1 ρ′ U Ht
= U −1 [H0 , ρ′ ]U + U −1 U Ht U −1 ρ′ U − U −1 ρ′ U Ht U −1 U= U −1 [H0 , ρ′ ]U + U −1 Ht′ ρ′ U − U −1 ρ′ Ht′ U
なのでリウヴィル方程式はU −1 [H0 , ρ′ ]U + iℏU −1∂ρ′U = U −1 [H0 , ρ′ ]U + U −1 Ht′ ρ′ U − U −1 ρ′ Ht′ U∂tiℏ∂ρ′= Ht′ ρ′ − ρ′ Ht′∂t= [Ht′ , ρ′ ]これを t 積分して、積分範囲を −∞ ∼ t として −∞ での解を ρ0 とすれば1
ρ (t) = ρ0 +iℏ′∫ t−∞dt′ [Ht′′ , ρ′ (t′ )]dt′となります。ρ に戻せば

U −1 ρ′ (t)U = U −1 ρ0 U +.ρ(t) = ρ0 +1iℏ

1= ρ0 +iℏ

1iℏ

∫ t−∞

∫ t−∞

∫ t−∞

dt′ U −1 (t)[U (t′ )Ht′ U −1 (t′ ), U (t′ )ρ(t′ )U −1 (t′ )]U (t)dt′dt′ U −1 (t)U (t′ )[Ht′ , ρ(t′ )]U −1 (t′ )U (t)dt′′′dt′ e−iH0 (t−t )/ℏ [Ht′ , ρ(t′ )]eiH0 (t−t )/ℏ dt′

1 次近似として ρ0 を第二項に入れれば
1
ρ(t) = ρ0 +iℏ∫ t−∞′′dt′ e−iH0 (t−t )/ℏ [Ht′ , ρ0 ]eiH0 (t−t )/ℏ dt′10古典論のときと比べて、(iℏ)−1 がいて、ポアソン括弧が交換関係になり、exp[∼] が両側にいるようになっています。

11