数学・物理通信
5 巻 12 号

編集

2016 年 2 月

新関章三・矢野2016 年 2 月 5 日

忠1

目次 (Contents)
1. ニュートリノ振動
小出義夫

2

世戸憲治

11

世戸憲治

21

2. 自然界で目にする波動現象

3. 自然界で目にする波動現象 (2)4. 編集後記
新関章三，矢野忠

27

Yoshio KOIDE

2

Kenji SETO

11

Kenji SETO

21

Shozo NIIZEKI and Tadashi YANO

27

1. Neutrino Oscillation

2. Wave Phenomena Appearing in Nature

3. Wave Phenomena Appearing in Nature (2)4. Editorial Commentsニュートリノ振動
小出義夫 1

Neutrino Oscillation
Yoshio Koide 2
1

はじめに2015 年度ノーベル物理学賞を東京大学宇宙線研究所の梶田隆章さんが受賞した．それで世間的にもニュート

リノが質量をもつことの実験的な証明である，ニュートリノ振動に社会的にも関心が起きたと思う．ところが
インターネットで調べてみてもそのニュートリノ振動についての本質的な説明があまりないように思われる．
それで物理とか数式がわかる人のためにニュートリノ振動について説明をされている，小出義夫氏の解説をこの『数学・物理通信』に転載することをお願いしたところ，快くご了承を頂いたので，ここに掲載する．
この解説はもともと小出氏の著書『物理学は何を発見してきたか？』(開成出版，2006）の付録 B で『量子力学的現象の例』というタイトルで述べられたものである．旧稿を転載することに同意して下さり，編者の求めに応じて加筆訂正をして下さった著者の小出義夫氏に深く感謝します．
なお，転載するにあたり，文章をいくつかの節を分けて読みやすいように変更を加えた．またスーパーカミ
オカンデの概念図と実際に小出氏が撮影された写真がそれぞれ一葉づつあったが，残念ながらここではこれらを割愛させてもらった．
なお，本文の脚注になっていた部分を付録 1, 2 にまわし，付録 4 として量子力学の基本法則を述べた．編集者がもともとあった付録も書き変えたので，付録の文責は編集者にある．
（編集者：矢野）2

ニュートリノとは
自然界には，ニュートリノと呼ばれる奇妙な粒子の存在が知られている．それは，３種類，νe ，νµ ，ντ が存在していて，電気的に中性で，質量は他の素粒子と比べて極端に小さく 3 ，また，物質との相互作用も極端に弱い（従って，物質中をほとんど自由に貫通してしまう）．これら３つのニュートリノ，νe ，νµ ，ντ は，弱い相互作用と呼ばれる素粒子反応において，荷電レプトンと呼ばれる粒子，電子 e− ，ミュー粒子 µ− ，タウ粒子 τ − とそれぞれ対になって関与するので，エレクトロンニュートリノ νe ，ミューオンニュートリノ νµ ，タウ
ニュートリノ ντ とそれぞれ名前が付いている．

1 静岡県立大学名誉教授
2 koide@kuno-g.phys.sci.osaka-u.ac.jp
3 各ニュートリノの質量は，現在でも正確には分かっていない．３つのニュートリノの内で，一番重いものでも 0.05 × 10−9 GeV（GeV はエネルギーの単位）くらいと推定されている．これに対して，ニュートリノを別とすれば，一番軽い素粒子である電子ですら，0.5×10−3 GeVの質量を持つ．ニュートリノはこれよりも 100 万分の１の軽さであるということになる．

2

これらニュートリノは，自然界の色々な場所で発生する．例えば，太陽の放射するエネルギーは，核融合反応によるものであるが，このとき，エレクトロンニュートリノ νe が発生する．ところが，地球上で νe を観測す
ると，太陽で核融合反応が起きて νe が放出されたと考えて計算した理論値とは合わず，予想される半分程度の数の νe しか観測されない．
（これは「太陽ニュートリノ問題」と呼ばれる．
）また，宇宙からいろいろの粒子が地球に降り注いでくるが，それが地球の大気（地表から 10km 程度上空）と衝突するとき，ミューオンニュー
トリノ νµ を発生させる．しかし，これも地表で観測すると 4 ，予想されるより半分程度の νµ しか観測されない．
（これは「大気ニュートリノ問題」と呼ばれる．
）
これらの現象は，以下に紹介する「ニュートリノ振動」と呼ばれる現象であることが現在分かっている．ニュートリノ振動は，このような自然現象としてだけではなく，人工的に発生させたニュートリノを用いての実験からも，その存在が確かめられている．この「ニュートリノ振動」とは，どのようなメカニズムで起こるのか，
次節で，大気ニュートリノ問題の場合に例をとって，解説しよう．

3

ニュートリノ振動—大気ニュートリノの場合

3.1

直観的な説明大気ニュートリノ問題は，簡単に言ってしまえば，大気中（上空約 10km）で発生したミューオンニュート
リノ νµ が，観測所に届くまでにタウニュートリノ ντ へと姿を変えてしまうことが原因である．ただし，この言い方は誤解を与えるので正しくない．この現象は，例えば 100 個の νµ が生まれて，その内のいくつかが，νµ → ντ へと姿を変えるために，結果として νµ の数が減ったということではない．それなら，生まれた νµ が，距離を走るにつれて，どんどん減って行くはずである．ところが，「ニュートリノ振動」の名前からも想像できるように，ある場所で測ると νµ の数は大きく減少して観測されるが，さらに離れた場所で観測すると今度は全く減少がないと観測され，そして更に距離をとって観測すると再び減少が観測され，· · · と，まさに「振動」的であることが観測される．この現象は，古典物理学的な意味で，νµ という粒子が別の粒子 ντ へと姿を変えるのではなく，量子力学的に，ニュートリノの状態が，発生したときは |νµ ⟩ であったものが，しばらく時間がたつと |ντ ⟩ になり，更に時間がたつとまた |νµ ⟩ になり，· · · と，「状態」の振動を起こしていることが原因なのだ．以下，このことを，量子力学の数式を用いて説明しよう．

3.2

量子力学による説明

ニュートリノの生成や検出には，素粒子反応を利用するので（即ち，荷電レプトン e− や µ− ，τ − を関与させての反応），その粒子状態は，|νe ⟩ か |νµ ⟩ か，あるいは |ντ ⟩ かに限られる．これら |νe ⟩，|νµ ⟩，|ντ ⟩ を「弱い相互作用の固有状態」と呼ぶ．
（量子力学における「状態」という概念は，初心者でなくとも，大変理解が難しい．そこで，通常の数値を表す記号とは区別して，「状態」を表すために |(· · · )⟩ という記号がよく用いられる．
）
4 正確には，「地表近くの地下」と言うべきであろう．例えば，日本には，岐阜県神岡街に「スーパーカミオカンデ」という巨大なニュートリノ観測装置が存在するが，これは，地下 1000m のところに設置されている．本当の意味で地表に観測装置を置くと，色々な宇宙からの粒子がすべて入り混ざって観測されてしまい，ニュートリノはそんなに数多く観測される粒子ではないので，ニュートリノを識別できなくなってしまう．地下に観測所を設けて置けば，ニュートリノ以外の粒子はすべて岩盤でシャットアウトできるので，岩盤を自由に貫通するニュートリノだけを観測できることになる．
（むろん，ほとんどのニュートリノは，その観測（検出）装置をも貫通して行ってしまう
ので，検出装置に引っかかるニュートリノは，そのほんの一部である）

3

一方，ニュートリノを観測していないとき（空間を自由に走っているとき），それはエネルギー固有状態にある．それを |ν1 ⟩，|ν2 ⟩，|ν3 ⟩ と書こう．これは時間に依存しない Schrödinger 方程式，

H|νi ⟩ = Ei |νi ⟩(i = 1, 2, 3)(3.1)

を満たす．ここで，Ei は，静止質量 mi を持つニュートリノ νi の相対論的エネルギーであり，Ei =√(mi c2 )2 + (cp)2(3.2)

で与えられる 5 ．
時刻 t におけるエネルギー固有状態 |νi ⟩ と，時刻 t = 0（そのときニュートリノが生まれた）におけるエネ
ルギー固有状態 |νi (0)⟩ とは，|νi (t)⟩ = e−i(Ei /ℏ)t |νi (0)⟩,(3.3)の関係がある 6 ．
むろん，偶然，弱い相互作用の固有状態 |νe ⟩，|νµ ⟩，|ντ ⟩ が，エネルギーの固有状態 |ν1 ⟩，|ν2 ⟩，|ν3 ⟩ と同一である（例えば，|νe ⟩ = |ν1 ⟩，|νµ ⟩ = |ν2 ⟩，|ντ ⟩ = |ν3 ⟩）ということもあり得るが，一般には別の状態と考えるべきである．大気ニュートリノ問題では，νν ↔ ντ の振動であるということが判明しているので，今 |νµ ⟩，|ντ ⟩ と |ν2 ⟩，|ν3 ⟩ の関係が，任意の時刻 t において，常に|νµ (t)⟩ = |ν2 (t)⟩ cos θ23 − |ν3 (t)⟩ sin θ23 ,(3.4)|ντ (t)⟩ = |ν2 (t)⟩ sin θ23 + |ν3 (t)⟩ cos θ23 ,(3.5)

であったとしよう 7 ．
（θ23 は，時刻 t に依存しないことに注意）(3.4),(3.5) を逆に解けば，|ν2 (t)⟩ = |νµ (t)⟩ cos θ23 + |ντ (t)⟩ sin θ23 ,(3.6)|ν3 (t)⟩ = −|νµ (t)⟩ sin θ23 + |ντ (t)⟩ cos θ23 ,(3.7)の関係と見てもよい．
今，ある時間 t に，ニュートリノの状態はどうなっているのかを，µ 粒子を関与させることによって検出することとしよう 8 ．この方法で検出できるのは |νµ (t)⟩ の状態のみである．この状態 |νµ (t)⟩ は，(3.4) から，|νµ (t)⟩ = |ν2 (t)⟩cosθ23 − |ν3 (t)⟩sinθ23 ,(3.8)

である．|ν2 (t)⟩，|ν3 (t)⟩ の時間変化は，(3.3) にしたがう：|ν2 (t)⟩ = e−i(E2 /ℏ)t |ν2 (0)⟩,(3.9)|ν3 (t)⟩ = e−i(E3 /ℏ)t |ν3 (0)⟩.

(3.10)

ニュートリノが発生した t = 0 では，ニュートリノ状態は |νµ (0)⟩ または |ντ (0)⟩ のはずである．したがって，
この |ν2 (0)⟩ および |ν3 (0)⟩ は，|νµ (0)⟩，|ντ (0)⟩ に戻してながめる必要がある．よって (3.8) は，(3.9),(3.10) を5 (3.2) の導出は付録 1 を参照せよ．
6 (3.3) の導出は付録 2 を参照せよ．
7 この関係は，３つのニュートリノについてのある統一的な法則に由来すると考えるが，ここでは単なる仮定として扱う．
8 ニュートリノを直接見ることはできない．電気的中性の粒子は見えない．検出可能の粒子は，荷電粒子のみである．

4代入し，さらに t = 0 とした (3.6),(3.7) を代入することによって，最終的に次のように書ける：|νµ (t)⟩ = e−i(E2 /ℏ)t |ν2 (0)⟩cosθ23 − e−i(E3 /ℏ)t |ν3 (0)⟩sinθ23= e−i(E2 /ℏ)t (|νµ (0)⟩cosθ23 + |ντ (0)⟩sinθ23 )cosθ23− e−i(E3 /ℏ)t (−|νµ (0)⟩sinθ23 + |ντ (0)⟩cosθ23 )sinθ23()= e−i(E2 /ℏ)t cos2 θ23 + e−i(E3 /ℏ)t sin2 θ23 |νµ (0)⟩()+ e−i(E2 /ℏ)t − e−i(E3 /ℏ)t sinθ23 cosθ23 |ντ (0)⟩= cµ (t)|νµ (0)⟩ + cτ (t)|ντ (0)⟩(3.11)

従って，t = 0 で発生した νµ ，即ち，|νµ (0)⟩ 状態が，t = t で νµ（即ち，νµ (t)）である確率 P (νµ (0) → νµ (t))
は，(3.11) から，νµ (0) → νµ (t) となる部分の係数

cµ (t) ≡ e−i(E2 /ℏ)t cos2 θ23 + e−i(E3 /ℏ)t sin2 θ23(3.12)

P (νµ (0) → νµ (t)) = |cµ (t)|2(3.13)の絶対値の２乗で与えられるので，P (νµ (0) → νµ (t)) = (e−i(E2 /ℏ)t cos2 θ23 + e−i(E3 /ℏ)t sin2 θ23 )(e+i(E2 /ℏ)t cos2 θ23 + e+i(E3 /ℏ)t sin2 θ23 )= cos4 θ23 + sin4 θ23 + (e−i(E2 −E3 )t/ℏ + e+i(E2 −E3 )t/ℏ )sin2 θ23 cos2 θ23
E2 − E3= cos4 θ23 + sin4 θ23 + 2cos
t sin2 θ23 cos2 θ23ℏ()
E2 − E3= cos4 θ23 + sin4 θ23 + 2 1 − 2sin2
t sin2 θ23 cos2 θ23
2ℏ
E
−
E
2
3
= (cos2 θ23 + sin2 θ23 )2 − 4sin2
t sin2 θ23 cos2 θ23
2ℏE2 − E3t,= 1 − sin2 2θ23 sin22ℏ

(3.14)を得る 9 (3.14) を見れば，P (νµ (0) → νµ (t)) が時間的に振動していることがわかるのだが，つぎの小節でそのことを具体的に検討して見よう．

3.3

具体的な検討もう少し，(3.14) を具体的に見るために，時間 t ではなく，ニュートリノの走った距離で，θ≡

E2 − E3t2ℏ

9 この式変形において，次の数学公式を用いた．

eiθ = cosθ + i sinθ ⇒ eiθ + e−iθ = 2cosθ.
cos2θ = 1 − 2sin2 θ,

sin2θ = 2sinθcosθ.

5(3.15)

の部分を見てみることにしよう．ニュートリノの質量はとびきり小さいので，ほぼ光速 c で走ると考えてよい
ので，走った距離 L とは

L = ctの関係がある．また Ei は (3.2) で与えられたが，ここでも mi が十分小さいことを考えると，近似公式(3.16)√
1+x≃

1 + 21 x を使って，√√
m22 c2m2 c2
E2 − E3 = cp 1 + 2 − cp 1 + 32
p
p
(
)
2 2
2 2
m c
m c
= cp 1 + 2 2 − 1 − 3 2
2p
2p
2
2
2
m − m3 3m − m23 4
= 2
c ≃ 2
c ,
2p
2E

(3.17)と書ける．ここで，最後の式変形は，m2 ，m3 が十分小さいので，E2 ≃ E3 ≃ cp = E と置いたことによる．
結局，(3.15) で定義した θ は，θ=

c3 m22 − m23
Lℏ
4E(3.18)と書ける．m2i を eV2 という単位で，E を GeV という単位で，また，L を km という単位とするとき，

θ = 1.27(∆m223 c4 )[eV2 ]L[km]
E[GeV](3.19)と書ける 10 ．ただし，ここで ∆m223 ≡ m22 − m23 と表した．
大気ニュートリノの場合は，E = 1 ∼ 10GeV と考えてよい．L は，上記で発生したニュートリノが地表に届くだけなら L = 10km と見なしてよいが，実は，ニュートリノは，地球の裏側で発生したものも地球の内部を貫通して観測所まで届く．地球の直径は，12,800km なので，このことを考慮すると L = 10 ∼ 13, 000km と考えるべきである．

P (νµ (0) → νµ (t)) = 1 − sin2 (2θ23 )sin2 θ,(3.20)

より，sin2 θ が最大のとき，即ち，θ = 12 π （および θ = 23 π, 52 π, · · · ）である．大気ニュートリノにおける E
および L の値に対して，θ = 21 π となる (m23 − m22 )c4 の値は，10−1 ∼ 10−4 eV2 程度であることが分かる．大気ニュートリノの観測データを詳しく分析することにより，現在，νµ → νµ 振動からの情報として，

sin2 2θ23 ≃ 1.00, (m23 − m22 )c4 ≃ 2.3 × 10−3 eV2 ,(3.21)ということが結論付けられている．
なお，同様にして，太陽ニュートリノのデータからは，νµ → νµ の結果として，

sin2 2θ12 ≃ 0.85,(m22 − m21 )c4 ≃ 7.5 × 10−5 eV2 ,という情報が得られている．

10 因子 1.27 の出てくる計算は付録 3 を見よ．

6

(3.22)

4

まとめまとめとして，もう一度，ニュートリノ振動について復習すると，ニュートリノ振動が起きる原因は，ニュートリノの生成および検出は，弱い相互作用の固有状態，|νe ⟩，|νµ ⟩，|ντ ⟩ として現れるのに対して，ニュートリノの観測を行っていない間（ニュートリノが生まれて，それが再び検出されるまでの間）は，エネルギーの固有状態 |ν1 (t)⟩，|ν2 (t)⟩，|ν3 (t)⟩ として，(3.3) に従ってその「状態」が時間的変化を行うということにある．このために，ニュートリノの状態が，(3.11) で示したように|νµ (t)⟩ = cµ (t)|νµ (0)⟩ + cτ (t)|ντ (0)⟩の形で変動するからである．注意すべきは，この現象は，量子力学的現象であり，|νµ (t)⟩ はただ１個のニュートリノについて記述された「状態」であるということである．決して，100 個の νµ など，統計学的な多数個の粒子集団を記述している「状態」ではない．だからこそ，時間とともに |cµ (t)|2 や |cτ (t)|2 の値が大きくなったり小さくなったり，振動を示すことになる．古典的物理学的発想でこれを解釈しようとすると νµ → ντ → νµ → ντ → · · ·などと粒子名が変化しているように見え，説明に困ることとなる．ニュートリノの種類を特定する実験では，
状態 |νe ⟩, |νµ ⟩, |ντ ⟩ を識別することになり，一方，観測をしていない間のニュートリノはエネルギー状態 |ν1 ⟩,|ν2 ⟩, |ν3 ⟩ によって状態の時間変化がなされている．それら固有状態は，別々であるという量子力学的「状態」を理解してはじめて説明可能となる．

5

付録

1

エネルギーと運動量の関係

特殊相対論において，4 次元空間のベクトルである空間座標 r および時間座標 ct に対応して,の運動量 p およびエネルギー E/c は４次元ベクトルを構成し，それらは

mi c2
Ei = √
1 − v 2 /c2

mi v
p= √
,
1 − v 2 /c2で与えられる．

(cmi v)2
1 − v 2 /c2()
2 2(mi c )v2
v2=
1− 2 + 2
1 − v 2 /c2cc(mi c2 )2=
1 − v 2 /c2(mi c2 )2 + (cp)2 = (mi c2 )2 +

= Ei2
を得る．

6

付録

2

エネルギー固有値状態

実は，エネルギー演算子 H は，

H = iℏ
7

∂
∂t

３次元空間

で与えられる．即ち，波動関数 ψ(t) の言葉で書けば，「状態」ψ(t) はiℏ∂
ψ(t) = Eψ(t)∂tを満たす．この ψ(t) はψ(t) = e−i(E/ℏ)t ψ(0)の形をしていればよいことは，実際に固有値方程式の左辺に代入してみることによって確認できる：

iℏ

∂ψ(t)
∂(e−i(E/ℏ)t )= iℏψ(0)∂t∂tE −i(E/ℏ)t= iℏψ(0)(−i )eℏ= Ee−i(E/ℏ)t ψ(0)= Eψ(t).

以上で (3.3) を確かめることはできたが，きちんと導くには次のようにする．そのために時間推進演算子 U (t)を次の式で導入しよう．

U (t) = e−(i/ℏ)Ht
(
)
(
)2
(
)3
(
)n
1
it
1
it
1
it
1
it
=1+
− H +
− H +
− H + ··· +− H+ ··· .
1!
ℏ
2!
ℏ
3!
ℏ
n!
ℏ
ここで H はハミルトン演算子である．(3.1) によれば

H|νi ⟩ = Ei |νi ⟩(i = 1, 2, 3)

H 2 |νi ⟩ = HEi |νi ⟩(i = 1, 2, 3)= Ei H|νi ⟩(i = 1, 2, 3)が成り立つから，同様に= Ei2 |νi ⟩(i = 1, 2, 3)

H 3 |νi ⟩ = Ei3 |νi ⟩(i = 1, 2, 3)

H 4 |νi ⟩ = Ei4 |νi ⟩(i = 1, 2, 3)

同様にして

···
であるから

U (t)|νi ⟩ = e−(it/ℏ)H |νi ⟩
[
]
(
)
(
)2
1
it
1
it
= 1+
− Ei +− Ei + · · · |νi ⟩
1!
ℏ
2!
ℏ
= e−(i/ℏ)Ei t |νi ⟩(i = 1, 2, 3)

8

が成り立つ．ここで粒子の始状態 ψ(0) はエネルギーの固有状態の重ね合わせで表すことができるから

ψ(t) = U (t)ψ(0)と表すことができる．
ここで，U (t) は t = 0 の状態から t = t の状態を導く演算子であるので時間推進演算子とよばれている．この時間推進演算子はユニタリー演算子であり，ハミルトニアン H はエルミート演算子である．

7

付録

3

定数の導出

(3.18) からθ=

c3 m22 − m23
Lℏ
4E(3.18)

であるが，これをつぎのように式を変形する．

1 ∆m2 c4
L, (∆m2 = m22 − m23 )4ℏc E∆m2 c4L= 1.27Eθ=

1である．ここで， 4ℏc= 1.27 である．

8

付録

4

量子力学の基本法則

量子力学を知っている人には余計な付録かもしれないが，量子力学を知らない方もおられると思うので念のために量子力学の基本法則についてまとめておく．
（量子力学の第一法則）
量子力学状態 ψ(x, t) は系の観測がなされない限り，時間に依存した Schrödinger 方程式iℏ∂
ψ(x, t) = Hψ(x, t)
∂t

(8.1)

によって時間的に変化する．ここで H はハミルトニアンである．

（量子力学の第二法則）
量子力学的状態は系の観測がなされるときには，状態の時間的変化は上の時間に依存した Schrödinger 方程式にしたがわず次のような確率的な変化をおこなう．すなわち，系の観測によって状態関数 ψ(x, t) は

ψ(x, t) =∞∑cn (t)un (x)

(8.2)

n=0

で表された一般の状態からその中のいずれかの一つの固有状態 un (x) へと遷移する．一般の重ね合わせの状態
からどの固有状態へと遷移するかは前もって予言することはできない．しかし，一般の重ね合わせの状態から特定の固有値状態への遷移確率はその固有状態の展開係数 cn (t) の絶対値の 2 乗 |cn (t)|2 に比例する．
量子力学での状態の時間的変化は第一法則と第二法則の２種類があることに気をつけてほしい．

9

また，状態関数 ψ(x, t) = T (t)u(x) が時間の関数部分 T (t) と空間の関数部分 u(x) とに変数分離できる場合には ψ(x, t) を時間に依存した Schrödinger 方程式に代入すればT (t) = Ae−i(E/ℏ)t ,A : 定数

Hu(x) = Eu(x)が得られる．ここで Hu(x) = Eu(x) を時間に依存しない Schrödinger 方程式という．ここで 1 次元の系のハ
ミルトニアン H はH=−ℏ 2 d2+ V (x),2m dx2V (x) :

ポテンシャル・エネルギー

(8.3)

である．

Schrödinger 流の量子化で座標表示を用いるとき，位置座標 x はそのままで，運動量 px は
px = −iℏ

d
dx

(8.4)

と x の微分演算子で表される．
直線上を運動する，質量 m の粒子のもつ古典力学での力学的エネルギーはH=

p2x+ V (x),2m

ポテンシャル・エネルギー

V (x) :

で表されるが，この式に上に与えられた運動量を表す微分演算子 (8.4) を代入すれば量子力学でのハミルトン演算子 (8.3) が得られる．

(2015.12.30)

10自然界で目にする波動現象
世戸 憲治 ∗

Wave Phenomena Appearing in Nature
Kenji Seto∗

1

はじめに
自然界で目にする波動現象というと，第 1 に，海の波を思い描くが，その他に，(1)

砂地に発生する風紋 (砂紋) がある．これは砂が風によって飛ばされるためにできるもので，その原理的なことはすでに解明されている．

(2)

北海道に住む筆者にとって，冬の雪はあたりまえであるが，雪原ではどんなに風が吹いても風紋らしきものは発生しない．これは雪があまりに軽すぎて，かなりの距離を飛ばされてしまうからなのか？

そのか

わり，車が走る雪道では，信号の手前 50m 位から信号を過ぎて 50m 位のあいだに，波長 10cm 程度の波型の凹凸ができてしまう．これはおそらく，信号の前後では，車が停止するにも，発進するにも，雪面に水平方向の大きな力を与えてしまうためと考えられるが，詳しいことはわかっていない．

(3)

これは波動現象には属さないのかもしれないが，秋になると空にイワシ雲が発生する．何故，このようなきれいな形になるのかと，見とれてしまうが，これはベナール対流で起こると言われている．雲が個々に独立してウロコ状のものはそれで説明できるとしても，それが横方向に長くつながって波のような形になったものは，やはり，波動としての要因が関与しているものと考えられる．

(4)

草原に風が吹くと，草が波を打ってしな垂れていくのを見ることがある．不思議なのは，風が吹く横方向にその位相がきれいに揃っていることである．

(5)

高速道路を走っていると，車の密度が高いところと低いところが発生し，この粗密が波動として伝播することがある．

などが挙げられる．
次節では，(1) の「風紋現象」について，Werner モデルと言われるものについて解説する．この部分は，「複雑系の物理学 (上)，砂丘と風紋」(田口善弘，琉球大学集中講義, 2004) に依るところが多いが，この中の記述は分かり難い部分もあるので，ここでは，でき得る限り，分かり易く丁寧に解説することを心がけた．また，3 節
では，(5) の「高速道路における波動現象」について，筆者なりの見解を述べることにする．

∗ 北海学園大学名誉教授

E-mail: seto@pony.ocn.ne.jp

11

2

風紋現象
初めに，砂地面に発生する風紋現象について述べる．風紋は，砂地面に風が吹き，砂が飛ばされることによって起こる現象である．ただし，この砂の「飛び」は風紋ができるための第 1 要因ではあるが，位相の揃ったきれ
いな風紋はそれだけではできない．きれいな風紋ができるためには，砂が高いところから低い方へ崩れ落ちる効
果も考慮する必要がある．ここでは，説明のし易さから，この砂の「飛び」を先に導入し，「崩れ」の方は後で導入することにする．
砂地面に 2 次元座標 (x, y) を設定する．風は x 軸の正方向に向かって吹くものとし，x 軸方向に波が発生す
るものとする．ただし，初めは，y 方向には一様とし，x 軸のみの 1 次元座標として扱っていく．このとき，座標 x, 時刻 t における砂面の高さを h(x, t) とする．ただし，この h(x, t) は風紋ができていないときの平均の高さを基準として測ったものとし，h は風紋の波長に比べ十分に小さいものとする．また，この座標 x, 時刻 t で砂が飛ばされる距離 L(x, t) は高さが高いところから飛ばされる方が大きくなるものとし，

L(x, t) = L0 + bh(x, t)(2.1)と仮定する．ここに，L0 は風の強さによって決まる量，また，b は正の比例定数とする．さらに，座標 x, 時刻t において，単位長さ，単位時間あたりに砂が飛ぶ量を N (x, t), また，着地する量を M (x, t) とする．砂は距離
L だけ飛ぶのであるから，x に着地する砂は，座標
ξ = x − L(ξ, t)(2.2)

から出発していることになる．ここで，右辺の L が L(ξ, t) であって，L(x, t) ではないことに注意する．ただし，これ以降，砂が飛んでいる時間は，h(x, t) の時間変化に比べ小さいものとして無視することにする．この式
によって，ξ は x の関数 ξ(x) となる．いま，dξ を微小長さとし，区間 [ξ, ξ + dξ] から，単位時間に，飛ぶ砂の量は N (ξ, t)dξ であり，これが，区間 [x, x + dx] に着地するとすると，それは，M (x, t)dx に等しいことから，

M (x, t) = N (ξ, t)

dξ
dx

(2.3)

という関係式が成り立つ．
ここで，時間発展に関する方程式をたてる．dt を微小時間として，この時間内に区間 [x, x + dx] で増える砂の量 M (x, t)dxdt から，減少する砂の量 N (x, t)dxdt を引いたものが，この区間における高さの変化による砂[]の量 h(x, t + dt) − h(x, t) dx となるので，

[

]
[
]
h(x, t + dt) − h(x, t) dx = M (x, t) − N (x, t) dxdt

(2.4)

が成り立つ．両辺を dxdt で割り，書き直すと，∂h(x, t)= M (x, t) − N (x, t)
∂t

(2.5)

となる．これは，流体力学における連続の方程式に相当する．さらに，(2.3) 式を代入すると，∂h(x, t)dξ
= N (ξ, t)− N (x, t)
∂t
dx
12

(2.6)

となる．ここで，dξ/dx の部分は (2.2) 式を x で微分し，(2.1) 式を用いると，となるので，(2.6) 式は，

dξ∂h(ξ, t)
=1−b
dx
∂x

(2.7)

[
∂h(x, t)∂h(ξ, t) ]= N (ξ, t) 1 − b− N (x, t)
∂t
∂x

(2.8)

となる．ここで，大きな仮定，砂が単位時間，単位長さあたりに飛ぶ量 N は座標にも時間にも依らずに一定

N (ξ, t) = N (x, t) ≡ N(2.9)∂h(x, t)∂h(ξ, t)
= −C
∂t
∂x

(2.10)C = Nb

(2.11)と仮定すると，方程式は，となる．ここで，定数 C をと定義する．さらに，この h(ξ, t) の部分を，(2.2) (2.1) 式を用い，h が十分に小さいことから，()
h(ξ, t) = h x − L0 − bh(ξ, t), t ∼= h(x − L0 , t)と近似してしまうと，方程式は，∂h(x − L0 , t)∂h(x, t)
= −C
∂t
∂x

(2.12)(2.13)と線形化されたものになる．この方程式で，もし，右辺の h の中で，x − L0 の L0 がなければ，これは，1 次元の右進行波を表す波動方程式である．この −L0 の存在は何を意味するのだろうか．ここで，k を正の波数，νを振動数として，

h(x, t) = ei(kx−νt)(2.14)

とおいてみる．これをこの方程式に代入すると，[]
ν = Ck cos(kL0 ) − i sin(kL0 )(2.15)と，ν は複素数で求められる．これを (2.14) 式に戻すと，

h(x, t) = e−Ck sin(kL0 )t · eik[x−C cos(kL0 )t](2.16)となる．実際は，この式の実部，および，虚部の線形結合が解となる．いずれにしても，sin(kL0 ) > 0 のとき
は，時間がたつとゼロに収束し，風紋はできない．逆に，sin(kL0 ) < 0 のときは，時間がたつと発散してしまう解になる．もちろん，風紋ができるためには，sin(kL0 ) < 0 であればよいのだが，このままだと，いったん風紋
ができ始めるとそれがどんどん進行していき安定性に欠けたものになってしまう．
これを止めるには，砂が表面を崩れ落ち，高さの変化をやわらげる効果を入れる必要がでてくる．ここでは，単位時間に砂が崩れ落ちる量は，高さの勾配に比例するものとし，正の比例定数を D として，D∂h/∂xとする．区間 [x, x + dx] を考え，微小時間 dt 内に，点 x + dx において，崩れでこの区間に入ってくる量D[∂h/∂x]x+dx dt から，点 x で崩れで出ていく量 D[∂h/∂x]x dt を引いたものが，この区間内の dt 時間におけ13

る「崩れ」による変化量を表すことになる．前に，「飛び」だけで導いた (2.4) (2.5) 式の右辺にこの「崩れ」による変化量を付加すると，新しい方程式は，∂h(x, t)∂ 2 h(x, t)= M (x, t) − N (x, t) + D∂t∂x2(2.17)となる．これが，「飛び」と「崩れ」の両方を取り入れたときの砂の総量保存の式となる．この新しく付け加わった 2 階微分の項を，ここでは，散逸項と呼ぶことにする．これは，熱伝導方程式における熱が温度の高いところ
から低い方へ移動する散逸効果と同じものと考えてよい．
前述したように，この式の右辺 1，2 項目は (2.13) 式右辺のように近似されるので，解くべき方程式は，∂h(x, t)∂h(x − L0 , t)
∂ 2 h(x, t)
= −C
+D
∂t
∂x
∂x2

(2.18)となる．この方程式に (2.14) 式を代入すると，[]
ν = Ck cos(kL0 ) − ik C sin(kL0 ) + Dk

(2.19)と求められ，これを (2.14) 式に戻すと，

h(x, t) = e−k[C sin(kL0 )+Dk]t · eik[x−C cos(kL0 )t](2.20)となる．これは，C sin(kL0 ) < 0 で風紋が成長しすぎるのを，Dk > 0 で妨げる効果になっている．実際の風紋では，C sin(kL0 ) + Dk < 0 で風紋が成長し，ある程度成長すると，D の値も大きくなり，この値がゼロに近づ
いていき，成長が止まると考えられる．なお，風紋の移動速度 V は，

V = C cos(kL0 )(2.21)となる．この風紋の速度 V は，sin(kL0 ) < 0 のとき，正負両方，あるいは，ゼロの値をとり得るので，必ずしも，風が吹く方向に進むとは限らない．
ここで，風紋ができるための条件をより詳しく求めてみよう．(2.20) 式右辺の第 1 指数関数に現れる因子を

f = C sin(kL0 ) + Dk

(2.22)と置いたとき，f > 0 ならば風紋はできず，風紋ができるためには f < 0 でなければならない．その境界となる f = 0 としたとき，k について正根をもつためのパラメータに対する条件を考えてみる．パラメータはC, D, L0 の 3 個あるが，変数 kL0 について考えると，それは，明らかに，D/(CL0 ) の値で決まる．したがって，パラメータとしては，D のみを変化させることで十分である．D が十分大きいときは，正根がない．また，

D が十分小さいときは複数個の正根を持つ．したがって，D の値を大きいところから小さい方に動かしたとき，初めて根を持つときは重根になっている．そこで，2 本の方程式，

df
=0
dk

f = 0,(2.23)

を解くことで，パラメータに対する条件が求められる．
これら 2 本の式から，三角関数の部分を消去し，k を求めると，

√
k=

C2
1
− 2
D2
L0
14

(2.24)

となり，これを df /dk = 0 の式に代入すると，

1,cos(α) = − √1 + α2

√
ただし，

α=

C2 2
L −1
D2 0

(2.25)となる．この方程式は解析的には解けないが，数値的に，正の最小根を求めてみると，

α = 2.396287 · · · ,

したがって，√
C
L0 = 1 + α2 = 2.596573 · · ·D(2.26)となる．これは，風紋ができるための最小条件で，一般には，

C
L0 > 2.596573 · · ·D(2.27)

であればよい．これから，C, L0 が大きいほど，また，D が小さいほど，風紋ができやすいことになるが，C, L0
をあまり大きくすると，許される波数 k の範囲が広がり，かえって，風紋ができにくくなる．また，D を小さくしすぎると位相の揃ったきれいな波にはならない．
これまで，座標１次元で述べてきたが，ここで，座標を (x, y) の 2 次元とし，方程式 (2.18) において，右辺最後の散逸項の部分の 2 階微分を，2 次元の Laplacian に拡張した( ∂2∂h(x, y, t)∂h(x − L0 , y, t)
∂2 )
= −C
+D
+
h(x, y, t)
∂t
∂x
∂x2
∂y 2

(2.28)とする．この最後の項は，砂が四方八方に散らばる効果を表している．波の位相を横方向 ( y 方向) に揃えてくれるのはこの項のおかげである．以後の 4 節で，この方程式に基づいた数値解法例を示す．

3

高速道路における波動現象
高速道路を走っているとき，車がある程度混んでくると，密度の高い部分と低い部分が発生することがある．

しかも，この粗密が波動として伝播する．この現象が起こる理由を考えてみよう．ここでは，車の流れを流体力学になぞらえてみる．場所 x, 時刻 t における車の密度を ρ(x, t), また，その速度を v(x, t) とする．dx を微小距離として，区間 [x, x + dx] を考える．点 x において，微小時間 dt のうちに，この区間に入ってくる車の台数は

ρ(x, t)v(x, t)dt であり，また，点 x + dx において，この区間から出ていく車の台数は ρ(x + dx, t)v(x + dx, t)dt[]
である．その差が，この区間内の dt 時間経過後の車の台数変化 ρ(x, t + dt) − ρ(x, t) dx となるので，[]
ρ(x, t + dt) − ρ(x, t) dx = ρ(x, t)v(x, t)dt − ρ(x + dx, t)v(x + dx, t)dt(3.1)

が成り立つ．両辺を dxdt で割り，書き直すと，[]∂ ρ(x, t)v(x, t)∂ρ(x, t)
=−
∂t
∂x

(3.2)

となる．これは流体力学における連続の方程式そのものである．
ここで，第 1 段階として，車の速度 v(x, t) は密度 ρ(x, t) の関数としよう．密度が小さいときは速度を上げ，密度が大きくなると速度を下げる．しかし，高速道路には速度制限があるので，密度が小さいといっても無制限に速度を上げるわけにはいかない．ここでは，最も単純化したモデルとして，v は ρ の 2 次式[
( ρ(x, t) )2 ]v(x, t) = v0 1 −
ρ0
15

(3.3)

と仮定してみる．ここに，v0 は制限速度，ρ0 は最大に車が詰まってしまい走れなくなったときの密度である．
この式から，ρ → 0 で v → v0 , また，ρ → ρ0 で v → 0 となるので，一応，所要の条件は満たされる．しかし，実際の高速道路で起きている現象は，これほど単純ではない．運転していて，先が詰まってきたときは速度を下げ，先が空いてきたときは速度を上げるということもある．これは，車の速度が，密度勾配にも左右されることを意味する．いま，密度 ρ が小さいところと，大きいところの 2 か所があって，これら 2 か所における密度勾配∂ρ/∂x が同じであったとしよう．密度が小さいところでは，かなりのスピードをだしているので，先が詰まってきたときは，大きく減速しなければいけないが，もともと密度が大きいところでは，あまりスピードをだしていないので，先が詰まってきたときの減速量も小さいはずである．ここでは，これらの事情を考慮し，第 2 段階として，(3.3) 式に密度勾配に依存する項を付け加え，[( ρ(x, t) )2 ]E ∂ρ(x, t)v(x, t) = v0 1 −−ρ0ρ(x, t) ∂x(3.4)と修正する．ここに，E は正の定数とする．この式で，分母に ρ(x, t) を付けたのは，同じ密度勾配であっても密度が小さいときの方が，減速量，あるいは，加速量を大きくするためである．
ここで，第 3 段階として，さらなる拡張をする．この密度勾配に関する項は，各運転者の直前の車の位置で決まると言えるが，実際の運転では，運転者が自分の車の数台前までの車の位置に注意して運転しなければならな
い．これは密度のより高階の微分係数までをも考慮する必要がでてくるということである．ここでは，密度の 2
階微分 ∂ 2 ρ/∂x2 も取り入れて，この 2 階微分が正のときは先が詰まってくることが予測されるので速度を下げ，2 階微分が負のときは先が空いてくることが予測されるので速度を上げるということにする．式で表すと，(3.4)
式にさらなる項を付け，[( ρ(x, t) )2 ]E ∂ρ(x, t)
F ∂ 2 ρ(x, t)v(x, t) = v0 1 −−−ρ0ρ(x, t) ∂x
ρ(x, t) ∂x2(3.5)と拡張する．ここに，F は適当な正の定数とする．
この (3.5) 式を (3.2) 式に代入すると，(∂ρ(x, t)∂ [
ρ2 (x, t) )]∂ 2 ρ(x, t)∂ 3 ρ(x, t)= −v0
ρ(x, t) 1 −
+E
+F
2
2
∂t
∂x
ρ0
∂x
∂x3

(3.6)

となり，これが解くべき方程式となる．この式右辺の 1 項目は非線形項であり，2 項目は前節でも導入した散逸項で，密度の急激な変化をやわらげる効果を持つ．また，3 項目は分散項と呼ばれ，密度を分散させる効果を持つ．
この方程式を一般的に解くことは難しいが，E, F のどちらか一方がゼロの場合の定常波解は簡単に求めることができる．

E ̸= 0, F = 0 の場合
方程式は，

(
∂ [
ρ2 (x, t) )]∂ρ(x, t)∂ 2 ρ(x, t)= −v0
ρ(x, t) 1 −
+E
2
∂t
∂x
ρ0
∂x2

(3.7)

となり，この方程式は，Burgers 方程式 *1 の次数を 1 つ上げたもので，簡単にその定常波解を求めることができる．c を定数として，

ξ = x − ct

ρ = ρ(ξ),*1 Burgers 方程式は波動関数を u としたとき，∂u/∂t + u∂u/∂x − ν∂ 2 u/∂x2 = 0 と書かれる．

16

(3.8)

と仮定すると，方程式は，

−c

dρ
d[ (
ρ2 )]
d2 ρ
= −v0
ρ 1− 2 +E 2
dξ
dξ
ρ0
dξ

(3.9)

dρ (
v0 )
= v0 − c − 2 ρ 2 ρdξ
ρ0(3.10)となり，一度，ξ で積分して，Eとなる．ここで，ξ → −∞ で ρ → 0, dρ/dξ → 0 とし，積分定数をゼロとした．0 < c < v0 とし，変数分離してから，さらに積分すると，

aeαξ
ρ= √
,
1 + e2αξ√
ただし， a = ρ0v0 − c,v0α=v0 − cE(3.11)と求められる．この解は，ξ → −∞ で ρ → 0, であり，ξ → ∞ で ρ → a となるキンク解である．このときの車の速度 v(x, t) は (3.4) 式にこの解を代入すると，場所にも，時間にも依存しない一定値 v = c となる．もっとも，この v が一定値 c となることは，この解を (3.4) 式に代入するまでもなく，(3.8) 式を仮定した (3.2) 式を積分することで，簡単に求められる．実際，(3.10) 式は v = c とした (3.4) 式になっている．これから，このキンクは，車の動きとともに形を変えずに速度 c で進行することになる．

E = 0, F ̸= 0 の場合
方程式は，(∂ρ(x, t)∂ 3 ρ(x, t)∂ [ρ2 (x, t) )]
+
F
= −v0
ρ(x, t) 1 −
∂t
∂x
ρ02
∂x3

(3.12)となり，これは，変形 Korteweg-deVries(mKdV) 方程式と少し形は違うが本質的には同じものである *2 ．この方程式の定常波解を求めるため，(3.8) 式を仮定し，ξ 積分をすると，

F

d2 ρ
v0
= (v0 − c)ρ − 2 ρ3dξ 2
ρ0(3.13)となる．ここで，ξ → ±∞ で ρ → 0, dρ/dξ → 0, d2 ρ/dξ 2 → 0 を仮定する．さらに，両辺に 2(dρ/dξ) を掛けて ξ で積分すると，( dρ )2dξ)
β2 (= 2 b2 − ρ2 ρ2 ,b√
ただし， b = 2a = ρ0√2(v0 − c)
,
v0

√
β=

v0 − cF(3.14)となる．これから dρ/dξ を求め，さらに，変数分離して積分すると，

ρ=bcosh(βξ)(3.15)と，これはソリトン型の解が求められる．このときの車の速度 v(x, t) は，キンク解のときと同様に，一定値

v = c となるが，E = 0 とした (3.5) 式にこの解を代入することでも確かめられる．したがって，このソリトン
は車の流れといっしょに速度 c で移動することになる．なお，このときのソリトンは通常のものとは違い，c の値によって，ソリトンの高さが高いほど移動速度は遅くなり，高さが低いほど移動速度は速くなる．

*2 mKdV 方程式は，波動関数を u としたとき，∂u/∂t + 6u2 ∂u/∂x + ∂ 3 u/∂x3 = 0 と書かれる．

17

4

数値解法によるグラフ表示

4.1

風紋の場合

初めに，2 節で述べた風紋の数値解法例を示す．方程式 (2.28) を差分方程式に変換するため，∆t を微小時間とし，この ∆t を単位として測った時間 t/∆t を改めて整数値の t とする．同様に，∆x, (=)∆y を微小距離とし，これを単位として測った距離 x/∆x, y/∆y を改めて，整数値の x, y とする．さらに，この方程式における距離 L0 ，および，係数 C, D を

L0 /∆x → L0 ,

C

∆t
→ C,
∆x

D

∆t
→ D
∆x2

(4.1)

と改めて置き直すことにする．この置き換えで，微分方程式 (2.28) は差分方程式に，]
C[h(x, y, t + 1) − h(x, y, t) = − h(x − L0 + 1, y, t) − h(x − L0 − 1, y, t)
2
[
]
[
]
+ D h(x + 1, y, t) − 2h(x, y, t) + h(x − 1, y, t) + D h(x, y + 1, t) − 2h(x, y, t) + h(x, y − 1, t) (4.2)と変換されるので，以後，この方程式を数値的に解くことにする．
ここでは，点 (x, y) を 200 × 200 の格子点上にとり，その境界では x 方向，y 方向共に周期境界条件を課すことにする．また，パラメータ L0 , C, D については，風紋が最もきれいにできる

L0 = 15,

C = 0.01,

D = 0.01(4.3)と設定した場合を求めてみる．その上で，t = 0 の初期値として，風紋の高さ h を各点ごとに乱数で与えておき，方程式 (4.2) に従って，時刻 t = 2300 となるまで求める．そのときの高さ h を以下の図 1 に示す．

図1

風紋

この図は，右方向に x 軸，斜め上方向に y 軸，上方向に高さ h をとって 3 次元的に示したものである．この図から，波の位相が y 方向にかなりの程度揃っているのがわかる．また，ところどころで，位相の切り替えが起こっているのも自然らしく，いかにも風紋という感じが見てとれる．このまま，時間を経過させると，高さはどんどん高くなり最終的には発散してしまうが，実際の風紋では前に述べたように高さが増えるにしたがい，係数D の値も増加し，ある程度成長したところで止まるようになっているものと考えられる．
18

4.2

高速道路の場合

つぎに，高速道路における波動現象を表す (3.6) 式の数値計算によるグラフを示す．風紋のときと同様に，∆t, ∆x を微小時間，微小距離とし，これらを単位として，時間，距離を測ることにし，係数 v0 , E, F を，v0

∆t
→ v0 ,
∆x

E

∆t
→ E,
∆x2

F

∆t
→ F
∆x3

(4.4)

と改めて置き直すことにする．この変換で，方程式は差分方程式

(
(
v0 [
ρ2 (x + 1, t) )
ρ2 (x − 1, t) )]ρ(x + 1, t) 1 −
−
ρ(x
−
1,
t)
1
−
2
ρ02
ρ02
[
] F[
]
+ E ρ(x + 1, t) − 2ρ(x, t) + ρ(x − 1, t) +ρ(x + 2, t) − 2ρ(x + 1, t) + 2ρ(x − 1, t) − ρ(x − 2, t) (4.5)
2ρ(x, t + 1) − ρ(x, t) = −に変換される．
この式に基づいて数値計算したものを，図 2，図 3 に示す．

図2

高速道路

図3

E = 0.01, F = 0.001

高速道路

E = 0.001, F = 0.01

これらの図では，パラメータ v0 , ρ0 , E, F の値を，それぞれ，
図2:

v0 = 0.01, ρ0 = 10, E = 0.01,

図3:v0 = 0.01, ρ0 = 10, E = 0.001, F = 0.01

F = 0.001(4.6)とし，図 2 と図 3 では，v0 , ρ0 の値は同じであるが，E, F の値が，それぞれ，10 倍異なるようにとってある．
図は前と同様，右横方向に x 軸，斜め上方向に t 軸，上方向に ρ 軸をとり立体的に示したものである．なお，x軸に関しては，0 から 200 までとし，その両端を自由境界とした．また，時間については，20000 ステップまで
計算し，その間，200 ステップごとに，計 101 本の曲線を描いている．初期値に関しては，図 2，図 3 共に，

ρ(x, 0) =6)
cosh (x − 40)/10
(

(4.7)

と設定した．図 2 からは，時間の経過とともに，密度 ρ の形はキンク解に近づいていくのがわかる．また，図 3
では，初期状態で設定したものがこの系での完全なソリトンではないために，時間発展の途中で，それが 3 個のソリトンに分裂しているのが見て取れる．先に述べたようにこの系では，高さが低いソリトンの方が速度が速い
ために，低い方が x 軸上を先に進んでいくのがわかる．

19

5

おわりに
ここでは，
「はじめに」のところで述べた (1)「風紋現象」と (5)「高速道路における波動現象」の 2 つの現象を扱ってみた．このうち，「風紋現象」については，Werner モデルと言われるものを紹介しただけであるが，こ
こでは，自己流ではあるが，でき得る限り分かり易い解説となるよう努めた．また，「高速道路における波動現象」についても，流体力学的に扱ったものが存在するが，ここでは，キンク解とソリトン解に結び付ける筆者の見解を披露してみた．もちろん，これは，実測データに基づくものではないので，その可否は何とも言えないが，理論的モデルとしては興味あるものになったと考えている．
しかし，雪国に住む筆者として，本当にやりたかったのは，(2) の「雪道における凹凸現象」である．初めは，これも風紋現象と同じようなものかと考えていたが，そう簡単ではないことがわかる．車が発進するときは，確
かに雪が後方に飛ばされるが，停止するときはほとんど雪が飛ばされることはないからである．しかも，どちらかと言うと，信号の手前で停止するところの方が，信号を過ぎた後よりも，凹凸が激しいように思えるからである．また，夏の暑いときにアスファルト路面が柔らかくなったときも，雪道ほどではないが，凹凸が発生することがある．この場合は，アスファルトはまったく飛んではいないので，風紋とは違う要因で起こっていることは確かである．雪道のある程度固まった雪は，流体なのか，弾性体なのか，あるいは，弾塑性体なのか，どう扱ってよいか，その辺のことからしてよくわからない．しかし，いつの日か，解決したい問題と考えている．

[ 謝辞 ]
今回もこの原稿を書くにあたって，京都大学名誉教授の中西襄先生にご精読いただき，大変有益なコメントをたくさんいただきました．ここに，謹んで感謝いたします．

20

自然界で目にする波動現象 (2)世戸 憲治 ∗

Wave Phenomena Appearing in Nature (2)Kenji Seto∗

1

はじめに
前回の「自然界で目にする波動現象」では，「風紋現象」と「高速道路における波動現象」を扱った．この風紋現象については，インターネットで検索してもたくさんの論文を見ることができるが，その中では，「砂丘現象」
についても述べられているものが多い．砂丘などというものは，どうせ，不定形だろうと考えていたが，どうやらそうではないらしい．その典型的なものは，バルハン砂丘と言われるもので，その模式図を，インターネット
から引用して図 1 に示す *1 ．

図1

バルハン砂丘の模式図

この図では，風は左から右に吹いており，三日月型の砂丘が形成されている．この砂丘の風上に面したところ
は傾斜がなだらかでその角度はおよそ 15 度，また，風下に面したほうは傾斜が急で，その角度はおよそ 30 度から 35 度とのことである．このバルハン砂丘の大きさは，大きなものでは，その高さは 30m 位，また，基底部の幅はおよそ 370m 位に達するということである．ちなみに，このバルハン砂丘というのは，1881 年にロシアの博物学者 Alexander von Middendorf によって，タジキスタンの砂漠で発見され，報告されたものである．このバルハンという命名も元はロシア語に由来するもので，英語では barchan dune, 日本語では三日月型砂丘と呼ばれている．初めに，このバルハン砂丘の写真を見たときは，何故にこのようなきれいな形になるのか不思議であったし，しかも，この形が数理モデルで再現できるということに驚いてしまった．今回も，新しいことはな
いのだが，自分でもこのきれいな形を再現してみようと思いたったしだいである．
∗ 北海学園大学名誉教授

E-mail: seto@pony.ocn.ne.jp*1 Barchan-Wikipedia21

2

方程式の導入
今回の砂丘を表す方程式は，前回の風紋を表す方程式をほんの少し修正しただけのもので，ほとんどの部分は同じものである．しかし，今回の独立性のため，敢えて，同じ部分も繰り返すことにする．
風紋，あるいは，砂丘は，砂地面に風が吹き，砂が飛ばされることによって起こる現象である．ただし，この砂の「飛び」は風紋や砂丘ができるための第 1 要因ではあるが，それだけでは，きれいなものはできない．きれいな風紋や砂丘ができるためには，砂が高いところから低い方へ崩れ落ちる効果も考慮する必要がある．
砂地面に 2 次元座標 (x, y) を設定する．風は x 軸の正方向に向かって吹くものとし，x 軸方向に波が発生す
るものとする．ただし，初めは，y 方向には一様とし，x 軸のみの 1 次元座標として扱っていく．このとき，座標 x, 時刻 t における砂面の高さを h(x, t) とする．ただし，この h(x, t) は風紋や砂丘ができていないときの平均の高さを基準として測ったものとし，h は，風紋の波長，あるいは，砂丘の大きさに比べ十分に小さいものとする．
初めに，砂の「飛び」の方を考える．この座標 x, 時刻 t で砂が飛ばされる距離 L(x, t) は高さが高いところ
から飛ばされる方が大きくなるものとし，

L(x, t) = L0 + bh(x, t)(2.1)と仮定する．ここに，L0 は風の強さによって決まる量，また，b は正の比例定数とする．さらに，座標 x, 時刻t において，単位長さ，単位時間あたりに砂が飛ぶ量を N (x, t), また，着地する量を M (x, t) とする．砂は距離
L だけ飛ぶのであるから，x に着地する砂は，座標
ξ = x − L(ξ, t)(2.2)

から出発していることになる．ここで，右辺の L が L(ξ, t) であって，L(x, t) ではないことに注意する．ただし，これ以降，砂が飛んでいる時間は，h(x, t) の時間変化に比べ小さいものとして無視することにする．この式
によって，ξ は x の関数 ξ(x) となる．いま，dξ を微小長さとし，区間 [ξ, ξ + dξ] から，単位時間に，飛ぶ砂の量は N (ξ, t)dξ であり，これが，区間 [x, x + dx] に着地するとすると，それは，M (x, t)dx に等しいことから，

M (x, t) = N (ξ, t)

dξ
dx

(2.3)

という関係式が成り立つ．
一方，砂の「崩れ」については，単位時間に砂が崩れ落ちる量は，高さの勾配に比例するものとし，正の比例
定数を D として，D∂h/∂x とする．
ここで，時間発展に関する方程式をたてる．初めに，砂の「飛び」の方を考慮すると，dt を微小時間として，
この時間内に区間 [x, x + dx] で増える砂の量 M (x, t)dxdt から，減少する砂の量 N (x, t)dxdt を引いたものが，この区間における「飛び」による変化量を表す．また，「崩れ」の方を考慮すると，この区間 [x, x + dx] で，微小時間 dt 内に，点 x + dx において，崩れでこの区間に入ってくる量 D[∂h/∂x]x+dx dt から，点 x で崩れで出ていく量 D[∂h/∂x]x dt を引いたものが，この区間内の dt 時間における「崩れ」による変化量を表す．これ[]ら，2 つの変化量を加えたものが，高さの変化による砂の量 h(x, t + dt) − h(x, t) dx となるので，

[ ∂h
[
]
[
]
∂h ]
h(x, t + dt) − h(x, t) dx = M (x, t) − N (x, t) dxdt + D−
dt∂x x+dx ∂x x
22

(2.4)

が成り立つ．両辺を dxdt で割り，書き直すと，∂h(x, t)∂ 2 h(x, t)= M (x, t) − N (x, t) + D
∂t
∂x2

(2.5)

となる．これが，「飛び」と「崩れ」の両方を取り入れたときの砂の総量保存の式となる．この式右辺の 2 階微分の項を，ここでは，散逸項と呼ぶことにする．これは，熱伝導方程式における熱が温度の高いところから低い方へ移動する散逸効果と同じものと考えてよい．
さらに，この式に，(2.3) 式を代入し，M (x, t) を消去すると，dξ∂ 2 h(x, t)∂h(x, t)= N (ξ, t)− N (x, t) + D
∂t
dx
∂x2

(2.6)

となる．ここで，dξ/dx の部分は (2.2) 式を x で微分し，(2.1) 式を用いると，dξ∂h(ξ, t)
=1−b
dx
∂x

(2.7)

[
∂h(x, t)∂h(ξ, t) ]∂ 2 h(x, t)= N (ξ, t) 1 − b− N (x, t) + D
∂t
∂x
∂x2

(2.8)

となるので，(2.6) 式は，となる．
ここまでは，前回のものとまったく同じである．前回の風紋のときは，ここで，単位時間，単位長さあたりに飛ぶ砂の量 N (x, t) を，座標にも，時間にも依存しない一定値と仮定した．今回は，これを止めて，風上に面したところでは砂は一定量のものが飛ぶが，風下に面したところでは，砂は飛ばないと仮定する．風は x 軸の正方向に向かって吹くので，言い換えると，高さの勾配が正のときは一定値 N0 ，勾配が負のところではゼロと仮定する．式で表すと，()
N (x, t) = N0 θ ∂h(x, t)/∂x(2.9)となる．ここに，θ は単位階段関数である．
この (2.9) 式を (2.8) 式に代入し，解こうとしたとき，ξ の存在が問題となる．これは，h を十分に小さいとした近似を用いるしか方法はないであろう．まず，(2.8) 式に含まれる h(ξ, t) のところは，(2.2) (2.1) 式を用いて，()
h(ξ, t) = h x − L(ξ, t), t ∼= h(x − L0 , t)

(2.10)と近似する．また，N (ξ, t) の部分については，初め，N (ξ, t) ∼= N (x − L0 , t) と近似して数値的に解いてみたが，これでは近似が強すぎて，砂の総量が時間とともに減少してしまうことがわかった．これは着地する砂の量が過少評価になるためらしい．そこで，近似をもう 1 段進めて，()
N (ξ, t) ∼= N x − L0 − bh(x − L0 , t), t(2.11)とする．この変数の中身を ξ 微分するところは，x 微分に置き換える．以上から，方程式は，()[∂h(x − L0 , t) ]∂h(x, t)∂ 2 h(x, t)= N x − L0 − bh(x − L0 , t), t 1 − b− N (x, t) + D
∂t
∂x
∂x2
となる．

23

(2.12)

これまで，座標１次元で述べてきたが，ここで，座標を (x, y) の 2 次元とし，方程式 (2.12) において，右辺最後の散逸項の部分の 2 階微分を，2 次元の Laplacian に拡張した

( ∂2
(
)[
∂h(x, y, t)∂h(x − L0 , y, t) ]∂2 )= N x − L0 − bh(x − L0 , y, t), y, t 1 − b− N (x, y, t) + D+ 2 h(x, y, t)
2
∂t
∂x
∂x
∂y
(2.13)とする．この最後の項は，砂が四方八方に散らばる効果を表しており，バルハン砂丘全体の形を形成するのに重()

要な働きをする．この方程式に含まれる N x − L0 − bh(x − L0 , y, t), y, t の項は (2.9) 式に代入すると，未知関数の中に未知関数をを含む複雑な方程式となり，もはや，解析的に解くことは不可能であろう．つぎの 3 節
で，この方程式に基づいた数値解法例を示す．

3

数値解法例
方程式 (2.13) を差分方程式に変換するため，∆t を微小時間とし，この ∆t を単位として測った時間 t/∆t を改めて整数値の t とする．同様に，∆x, (=)∆y を微小距離とし，これを単位として測った距離 x/∆x, y/∆yを改めて，整数値の x, y とする．さらに，この方程式における距離 L0 ，および，係数 N0 , b, D を

L0 /∆x → L0 ,N0 ∆t → N0 ,b/∆x → b,

D

∆t
→ D
∆x2

(3.1)

と改めて置き直すことにする．この置き換えで，微分方程式 (2.13) は前進型の差分方程式に，

(
)[
]]
b[
h(x, y, t + 1) − h(x, y, t) = N x − L0 − bh(x − L0 , y, t), y, t 1 − h(x − L0 + 1, y, t) − h(x − L0 − 1, y, t)
[
]
[ 2
]
− N (x, y, t) + D h(x + 1, y, t) − 2h(x, y, t) + h(x − 1, y, t) + D h(x, y + 1, t) − 2h(x, y, t) + h(x, y − 1, t)(3.2)と変換されるので，以後，この方程式を数値的に解くことにする．
ここでは，パラメータの値を，

L0 = 15,

N0 = 0.01,

b = 0.5,

D = 0.02(3.3)とおき，数値解法による h(x, y, t) の時間発展の様子を，以下の図 2-1 から図 2-6 までに示す．これらの図では，x, y を 200 × 200 の格子点上にとり，横方向に x 軸，斜め上方向に y 軸，上方向に h 軸をとって立体的に描いたものである．なお，x, y の境界は周期境界とした．

図 2-1

図 2-2

t=0
24

t = 5, 000

図 2-1 の t = 0 ステップでは，初期値 h(x, y, 0) を，−0.5 から +0.5 までの一様乱数で与えてある．これが時間発展して t = 5, 000 ステップの図 2-2 では，前回と同じような風紋ができている．

図 2-3

図 2-4

t = 10, 000

t = 15, 000

図 2-3 の t = 10, 000 ステップでは，グラフの中央部に 3 個のバルハン砂丘の卵とでもいうものができているのがわかる．ここではこれら 3 個に，左から，A, B, C という名前を付け，以後，これら 3 個のバルハン砂丘
がどのように時間発展するかを見ていく．t = 15, 000 ステップの図 2-4 では，これら，A, B, C のバルハン砂丘はいずれも x 軸の正方向に動いているのがわかる．ただし，図 2-3 で一番右にあった C は周期境界のため図2-4 では，一回りして左側にきている．また，A と C は成長して大きくなっているが，B はそれほど大きくはなっていない．

図 2-5

図 2-6

t = 20, 000

t = 25, 000

t = 20, 000 ステップの図 2-5 では B と C が合体してしまい 1 個のバルハン砂丘になっている，と言うより実際は，B は C ほど成長しなかったので，C が B を取り込んでしまった形になっている．この C が図 1 で示したバルハン砂丘の模式図に最も近い形になっている．t = 25, 000 ステップの図 2-6 では，これら A, B, C の他にたくさんのバルハン砂丘ができて，それらが，横方向 (y 方向) に連結してしまい曲がりくねった山脈のよ
うな形になっているのがわかる．このようなものをメガバルハンというが，一般に目にする砂丘はこの形なのであろう．
なお，これらの図では，x, y 軸のスケールと h 軸のスケールは見やすさのために変えてある．h 軸のスケー

25

ルは x, y 軸のスケールの 10/3 = 3.333 · · · 倍になっていることを断わっておく．

4

おわりに筆者としては，でき得る限り，解析的に解くことを，信条にしているが，今回のものは数値計算がメインになってしまった．数値計算をしているときに，同じパラメータでやっているのに，回を重ねるごとに，出てくるパターンが違っていることに気が付いた．初めは，パソコンが壊れたものと思い込んでしまいどうしようかと迷ってしまった．よく考えてみると，乱数を使っているので，回を重ねて実行すると，乱数を作るための種が異なってきて，違うパターンになることに気が付いた．今回は，これら何度か実行したうちの最もバルハン砂丘の模式図に近いものを採用することにした．初期値に乱数を使ってしまうとどうしても複数個のバルハン砂丘ができてしまい，最終的には，それらが連結したものになってしまう．初期値をうまくとると，図 1 に示すような 1個だけのバルハン砂丘を作ることができるはず．しかし，これはかなり難しい問題で，いろいろと試してはみたが，いまのところ，すべて失敗に終わってしまった．もっとも，今回の数値計算では，t = 25, 000 ステップまで計算するのに，およそ，40 分の時間を要する．その間，じっと，パソコンを眺めているのも結構つらいものがある．もっと高速のパソコンがあったらと恨めしくもなるしだいである．

[ 謝辞 ]
今回もこの原稿を書くにあたって，京都大学名誉教授の中西襄先生にご精読いただきました．ここに，謹んで感謝いたします．

26

編集後記
新年あけましておめでとうございます．
いつものことながら，新年を迎えたと思ったら早くも節分は目前です．今年も， 継続は力なり の精神で編集委員は頑張って仕事に励むつもりですのでよろしくお願い申し上げます．
今回の投稿は 3 編です．物理学の分野で東京大学の教授が昨年度のノーベル物理学賞を取った．内容
は，今まではニュートリノは質量ゼロが定説だったのが実は質量はあったという話である．
編集者（新関）は素粒子については全くの素人であるが，今号の解説論文のように詳しく分りやすい
解説があれば，分かったような気がする．
本年も健康には十分注意しましょう．
病院の患者癒やすやハンドベル（クリスマス）（2016 年 2 月１日，新関章三）
明けましておめでとうございます．今年も読者と投稿者の皆様によろしくおつきあいをお願いします．

12 月 30 日に 5 巻 11 号を発行した．その編集後記で「2015 年度のノーベル物理学賞がニュートリノ
振動の発見に与えられたので，その解説論文を今年中に発行したかったのだが，著者のご都合で来年回しとなった」と書いたら，すぐに著者の小出さんからメールが来て，すでに 12 月 8 日に原稿を送ったのだが，着かなかったかとの問い合わせがあった．
あわてて原稿が不着だったことを伝えるとその日のうちに原稿が再送されて来た．年が明けて早々に今度は世戸さんから 2 編の論文の投稿があった．これで 5 巻 12 号の原稿が十分にできたので昨年の積み残しということで 5 巻 12 号を発行する．しかし，編集者としては今回の措置は特例中の特例である．
今回の編集後記は新関さんの順番だが，5 巻 11 号の編集後記の経緯を説明をした．
メールの届かないことが最近多発しており，原稿を投稿されても私からの受取りの返事が行かなかったら，確認のメールかまたは再度の原稿の送付をお願いしたい．特に私が不在とかの理由がなければ即日または少なくとも 3 日以内には返事をするつもりなので．
（矢野

27

忠）