Vicsekモデルの初出論文を読む(簡潔で美しく深遠!)【論文紹介】#24
今回はまた、古典的な名論文を紹介したい。
魚や鳥や草食哺乳類の群れなど、自己駆動粒子系を記述するモデルの代表例であるヴィチェック(Vicsek)モデルの初出論文である。本論文にはVicsek modelとは書いてないが、筆頭著者のTamás Vicsekの名前からのちにそう呼ばれるようになった。
今や界隈では知らない人はいない超有名モデルだが、本記事では知らない人にも伝わるように書くのでご安心ください。(というか本論文自体も、当時Vicsekモデルが存在しない世界にVicsekモデルを知らしめるために書かれた論文なのでかなりわかりやすい。この論文紹介シリーズは、物理に興味のある高校生も視野に入れているが、本記事を読むだけでなく”初めての英語論文読破”にも本論文はお勧めしたい。)
今回紹介する論文はこちら。
Novel Type of Phase Transition in a System of Self-Driven Particles
自己駆動粒子系における新しいタイプの相転移
0611743 (arxiv.org)
(↑こちらは2006年11月29日にarXivで無料公開されたものだが、初出のPhysical Reviewでは1995年8月7日掲載である。Phys. Rev. Lett. 75, 1226 (1995) - Novel Type of Phase Transition in a System of Self-Driven Particles (aps.org))
まず、読む前から感じさせられるこの論文のすごさは、本文が8ページという短さで、今や5000以上の論文で引用されているという驚異的な質の高さである(最後にまとめられている図7ページと参考文献欄2ページを含めると計17ページ)。
しかもその内容は、プログラムを組めば誰でもパソコンで実験できる計算機実験のモデルについてである。論理演算の宇宙(ありうるすべての集合)の中からこんな素晴らしい宝石を掘り出したのかと驚く。
(Pythonでの実装はこちら↓。(他の言語でもググれば出てくると思う))
Pythonプログラミング（モンテカルロ・シミュレーション） (tohoku.ac.jp)
Vicsekモデルが対象とする自己駆動粒子系とは、魚や鳥や草食哺乳類の群れである。各個体は一定の速さで自走する粒子とみなし、近隣の個体と相互作用して向きを変えるとして、単純なルールで群れの動きを再現するのである。
動物は群れを成すと、天敵から逃れたり、餌を取るための移動が効率的になることが知られているが、そのことが理論的に詳しく理解できたのはVicsekモデルのおかげである。他にも、Vicsekモデルを応用した研究は今やまとめ切れないほどある。
シミュレーションモデルとしては次のようなものである。
1辺の長さLの正方形の周期境界条件のフィールドに一定の速さ v で自走する粒子がたくさんいる。各粒子は、距離 r 以内にいる他の粒子と時間 Δt の間隔で相互作用する。(下の図1a参照)
初期条件は N 個の粒子がフィールド上にランダムに配置され、各粒子の移動の向き θ ($${-\frac{\pi}{2}<θ<\frac{\pi}{2}}$$)はランダムだが、移動の速さは v で全員同じとする。
向き θ は、時間ステップ1回ごとに、
$${\langleθ(t)\rangle_r=\arctan\frac{\langle\sin(θ(t))\rangle_r}{\langle\cos(θ(t))\rangle_r}}$$
だけ変化する。ここで$${\langle\sin(θ(t))\rangle_r, \langle\cos(θ(t))\rangle_r}$$は、近隣 r の範囲内の粒子の向き θ のsin, cosを取ったものの平均値である。これは、近隣 r 以内にある粒子の向きの平均と同じ向きになろうとすることを表す。
ただし、ここには幅 η 以内のランダムなノイズ Δθ が入るとする。(このノイズは、磁性体を記述するイジングモデルでの温度(熱雑音)に相当する)
合わせると、時間ステップ1回ごと(t→t+1)の向き θ の変化は次のように書ける。(式番号は論文と同じ)
$${θ(t+1)=\langleθ(t)\rangle_r+Δθ}$$ , (2)
粒子には通し番号が降られ、i番目の粒子の位置は$${\bm{x}_i}$$と書き、その時間ステップ1回ごと(t→t+1)の変化は次のように書ける。
$${\bm{x}_i(t+1)=\bm{x}_i(t)+\bm{v}_i(t)・Δt}$$ , (1)
ここで$${\bm{v}_i(t)}$$は、一定の自走の速さ v に、先ほど見た向き θ を加味した速度ベクトルである。
つまりどの粒子も、時間ステップ1回ごとに$${\bm{v}_i(t)・Δt}$$だけ直線的に進んでは、向きを$${\langleθ(t)\rangle_r+Δθ}$$だけ変えるという、カクカクと折れ曲がった経路を描くことになる。(下の図1aがわかりやすい)
ちなみにこのシミュレーションモデルは、自走の速さ v をゼロにすると、磁性体におけるスピン(原子レベルの磁石)の向きの揃い具合をモデル化したXYモデルと同じになる。
また、v を無限大にすると、粒子の相対位置は毎回ランダムに変わる(向きが非常に近い2粒子でもΔtの間に十分遠くまで離れてしまう)ので、XYモデルの平均場近似と同じになる。
本論文では v を0.003～0.3まで変えてみたが、その範囲では結果はあまり変わらなかった。以下では、粒子の密度 ρ(=N/L^2)とノイズ幅 η の大きさを変えた結果を示す。まず図1a~dは、各粒子の速度ベクトルを矢印で表し、直近20ステップの軌跡も付けたスナップショットである。いずれも自走する速さは v=0.03、粒子数は N=300 で固定だがフィールドの1辺の長さLを変えて密度を変えた。変更したパラメータは以下の通り。
a: フィールドの1辺 L=7, ノイズ幅 η=2.0
b: フィールドの1辺 L=25, ノイズ幅 η=0.1
c: フィールドの1辺 L=7, ノイズ幅 η=2.0
d: フィールドの1辺 L=5, ノイズ幅 η=0.1
(a)密度が低くノイズが大きい場合 (b)密度が高くノイズが大きい場合
(c)密度が高くノイズが大きい場合 (d)密度が高くノイズが小さい場合
密度が高くノイズが小さい図1dでは見事に、ほとんどの粒子が同じ方向に移動するようになっている。これは、近隣と同じ向きになろうとする相互作用のおかげで起こる非自明な現象だ。初期状態はランダムでどの向きも対称(特別な向きはない)のに、自発的に対称性が破れて向きが揃って秩序的な運動をするようになった。
ここで、粒子の全体的な動きを定量的に評価するために、次のような値を考える。
$${v_a=\frac{1}{N・v}|\sum_{i=1}^N\bm{v}_i|}$$(3)
まず絶対値の中身は、全粒子の速度をベクトル的に足している。その足しあわされた結果のベクトルの向きを無視して大きさだけを見るのが絶対値記号。それをNで割ることで平均値とし、vで割ることで0～1の範囲内になるように規格化している。
これが何を意味するかというと、粒子の向きがランダムの場合は絶対値の中身がベクトル的に打ち消し合ってゼロになる。向きがそろっていれば、$${\bm{v}_i}$$はNの分足し上げられ、$${N・v}$$で割ると1に近い値になる。これは0なら向きがバラバラで、1なら完全に向きが揃っているという、”秩序具合”を表す値だと解釈できるので$${v_a}$$をオーダーパラメータと呼ぶ。イメージとしては↓この動画(3:56~)がわかりやすい。
オーダーパラメータを用いた解析としてまず、密度 ρ は一定でノイズ幅 η を変えて調べてみたのが図2左、ノイズ幅 η は一定で密度 ρ を変えてみたのが図2右。
粒子数を増やすと”解像度”が上がって相転移がシャープに見えてくる
右：密度ρ(横軸)を変えたときのオーダーパラメータ(縦軸)の変化
まずノイズ幅 η を小さくすると、向きが揃いやすくなってオーダーパラメータが大きくなるのが見て取れる。それも粒子数を多くすると(密度は変わらないように、粒子数に比例するようにフィールドの面積も大きくしている)、ある臨界点で急激に立ち上がるような相転移的な振る舞いが見えてくる。粒子数はいわば、ドット絵のピクセル数(≒解像度)とも解釈できるので、粒子数が多いほうがより現象を精細に見ているといえる。
この相転移的な振る舞いは、磁性体において、ある臨界温度より温度を下げると自発磁化が発生する(磁化した鉄は臨界温度より加熱すると磁化を失うといったほうがイメージしやすいかな)のと同様の現象である。↓これの図3のような。
https://event.phys.s.u-tokyo.ac.jp/physlab2021/pdf/condmat_Ising.pdf
密度 ρ を変えた場合も同様に(左右反転しているが)、密度を上げれば相互作用する数が増えることによって急速にオーダーパラメータが立ち上がるという相転移的な振る舞いがみられる。
図2のオーダーパラメータが急激に立ち上がる部分を、もう少し詳しく見てみたのが図3である。
図3のグラフは両対数で、$${v_a\propto(η_c(L)-η)^β}$$とした場合の臨界指数 β を求めるためのプロットとフィッティング直線である。$${η_c(L)}$$は立ち上がり始める η の臨界値で(これは図2左の、粒子数を変えたデータから推測される、粒子数無限大極限の推測値である)、そこを原点とした関数の形を推定する。
臨界指数 β は、もし0.5だったら、$${y=\sqrt{ax}}$$のお馴染みの無理関数の形になるが、それを両対数軸に書くと傾き0.5の直線になる。図3左では、臨界指数がいくつになるかわからないデータを両対数軸でプロットしてみて、直線でフィッティングしてみると傾きが0.45になったという図。その結果から、図2左の臨界点から立ち上がるところの曲線は$${v_a\propto(η_c(L)-η)^{0.45}}$$となる。これは$${y=\sqrt{ax}}$$の無理関数より少し立ち上がりが遅いということである。
同様に、密度ρを横軸に取った図3右からは、$${v_a\propto(ρ-ρ_c(L))^{0.35}}$$であり、ηよりも立ち上がりが遅いということがわかる。
さらに図3左から読み取れる興味深いこととして、系サイズが大きくなるにつれて、データが直線のフィッティングにきれいに乗る領域が左に広がっている。このことは、相互作用する範囲は半径 r の範囲内だけのはずだが、粒子が自走することで相対位置が変わることによって、実質的に相互作用できる範囲が広がっているためと解釈できる。
ノイズが少ない側でもフィッティングに乗るかどうかは、相互作用が r を超えた広範囲にわたっているかどうかにかかっている。粒子数が40だとそもそも広い範囲での相互作用の効果が見えず、フィッティングから外れてしまう。広範囲にわたっていないと、図1(b)のように、小さなクラスタ内では向きが揃っているが、別のクラスタでは揃っていないという状況になり、オーダーパラメータが大きくならない。
このことは、粒子が動かないとした場合のXYモデルには見られない新奇な現象である。
また、自走する速さ v をゼロにするとXYモデルと一致するなら、上記の η と ρ の肩の2つの臨界指数は(XYモデルの既存研究により)同じ値になるはずである。それが図3に示す通り一致しないのは、XYモデルとの唯一の違いである自走することによる効果である。
このモデルは単純であるにもかかわらず、ランダム状態での向きの対称性を自発的に破り、輸送がない状態から全体的な輸送を湧き起こす、運動論的相転移を含む、豊かで現実的なダイナミクスをもたらすことを示す。
以下は私の感想だが、本論文は、今後の応用が効きやすいように基本的な部分を見抜いてそこだけを鋭く簡潔に切り抜いて記述したという、現象の描写の上手さも際立っていると思う。
自走する速さ v をゼロにするとXYモデルのようなミクロスケールを記述するモデルに共通した部分もあり、動物の群れのようなマクロスケールを記述する可能性も広く持っているという、現実世界の非常に上手い抽象化をしたモデルとも言えるし、ミクロとマクロの共通点を見通しているような深遠な洞察も感じる。
少し前にTwitter(現X、以下Twitter)でこんな素敵な言葉を見た。
>理系の人達は眼をしっかりと見開いて、皆で摩訶不思議なこの世界をじっくりと観察し、その驚異に感嘆しまくっている
理系の人達は眼をしっかりと見開いて、皆で摩訶不思議なこの世界をじっくりと観察し、その驚異に感嘆しまくっているわけですが、このような主張をされる方は、一人で目を閉じていたいというお気持ちなのかなと、思われました。この宇宙の真実を深く堪能したいかどうかという差かなと。 https://t.co/Z1ahcuLCr4
— Masahiro Hotta (@hottaqu) August 14, 2024
この宇宙にはまだまだ、目を見張るような、将来に人類の至宝と呼ばれるような、素晴らしい理論の宝石がどこかに発見されず眠っていることだろう。
そんな、論理的宇宙の広大さを感じさせられる、簡潔で美しく深遠な論文だと思う。
(余談)
今回の論文にぴったりな音楽はこれだと思う。広大な論理的宇宙の中で輝く宝石のようなイメージ。
最後まで読んでいただきありがとうございます。気に入っていただけたら「スキ」ボタンを押したりフォローしたりしていただけますと私のモチベが上がります。内容についての質問や感想もお待ちしております。
また、研究機関には所属せずにやっておりますので、有料ジャーナルのアクセス料程度のサポートをいただけるとありがたいですし、今後の記事の質が上がるかもしれません。どうぞよろしくお願いいたします。