球の衝突と密度変化だけから起こる相転移(牛を球とみなして何が嬉しいのか)【論文紹介】#15

例えば、ペットボトルに小石を入れて振る。
この中で起こる現象をよく見ると実は未解明の現象がある。
パラメータとして、小石の数・大きさ、振る振幅・振動数を変えると、小石のほとんどがペットボトルの振動と同期して振動する状態と、同期せずばらばらの振動数でぶつかり合う状態があり、パラメータ空間中での状態の遷移が非常に普遍的な相転移を示すのである。

今週紹介する論文はこちら。
Interplay between an absorbing phase transition and synchronization in a driven granular system
外力のある粉粒系における吸収状態転移と同期現象の相互作用
2401.12817.pdf (arxiv.org)

吸収状態転移(APTs: Absorbing phase transitions)は、凝縮系、伝染病、地震、生態学、化学反応など、非平衡系に広く見られる。
一番わかりやすい例としては、感染速度s、治癒速度tの感染症があったとすると、s>tなら感染は広がっていくが、s<tならそのうち感染者数がゼロになり、それ以降感染者は発生しなくなる(元の状態には戻らない)という吸収状態になる。他にもこの文献に割と詳しく書いてある。ja (jst.go.jp)
吸収状態転移は、先行研究[7-9, 15, 19, 22]では転移は連続的に起こり、[10, 26, 27]では転移は不連続に起こることが分かっており、流体力学や弾性に起因する相互作用がこの効果の原因ではないかと議論されている[28]。
本研究では、パラメータを変えることによって転移が連続になったり不連続になったりする興味深い系を考える。

一辺 L の正方形の2枚の板を幅hだけ開けて平行にし、その間に直径 σ、質量 m の球をたくさん入れる。2枚の板は、z方向(板に垂直の方向)に振動数 f で振動できるものとする。
球と板や球同士の衝突は完全弾性衝突ではなく、z方向の運動エネルギーは板から供給されるが、水平方向の速度はだんだん失っていき、次第に上下動だけになる(図1左)。あるいは球同士の衝突で垂直方向の速度が水平方向の速度になることもある(図1右)。

図1　シミュレーションする状況
粒は上下の壁との衝突時に、接線方向の摩擦力により水平方向のエネルギーを失い、次第に上下動だけになる(左)。
粒同士の衝突時は、板の振動から得られる垂直エネルギーが粒子の速度の水平成分に伝達される(右)。

球の半径 σ = 2.5 mm、振動数 f = 53 Hz、質量 m = $${6.54 × 10^{-5}}$$ kgはこれ以降固定する。シミュレーションされた球の数 Nは、$${10^3}$$から$${2×10^4}$$の範囲である。エネルギー単位(1個の球が持つおおよそのエネルギー)は $${m(σf)^2 = 1.15 × 10^{-6}}$$ J で与えられる。水平方向の平均運動エネルギーTは粒体温度(granular temperature)と呼ばれ、次のように書ける。
T=$${\frac{1}{T}\int_{0}^{T}\frac{1}{2}m＜v(t)_x^2+v(t)_y^2＞dt}$$
難しそうな見た目の式だがなんてことはない、運動エネルギー($${\frac{1}{2}mv^2}$$)の水平成分を全部足して時間$${T}$$で割って、運動エネルギーの時間平均を出しているだけだ。
また、板の面積$${L^2}$$に対する球全部の断面積$${\frac{Nπσ^2}{4}}$$の割合を、充填率$${\phi = \frac{πNσ2}{4L^2}}$$と定義する。粒子数 N を増やすと$${\phi}$$も増える。
下グラフは、板の間隔 h と球の直径 σ の割合を(a)h/σ=1.88と(b)h/σ=1.63にし、振幅Aと球の直径σの割合を0.105～0.121(黄緑～青)、球の充填率$${\phi}$$(横軸)を0.05～0.45まで変えて粉体温度T(縦軸)が最終的に落ち着いた(stable stateの)値Tssをプロットしたものである(各点において、球の初期条件はランダムで、Tが落ち着くまでシミュレーションを行った)。

図2　充填率(横軸)の変化による粉体温度(縦軸)の転移
板の間隔hと球の直径σの割合を(a)h/σ=1.88と(b)h/σ=1.63、
振幅Aと球の直径σの割合を0.105～0.121(黄緑～青)で変化させた

まず、ぱっと見でわかる通り、(a), (b)ともに充填率$${\phi}$$を変えると、粉体温度(水平方向の運動エネルギー)がゼロになる相と有限の値を取る相で相転移が起こっている。特に、粉体温度ゼロの相は、ランダムな初期状態からしばらくは球同士が衝突して水平方向の速度があるが、一度粉体温度ゼロになれば、球は全て完全な上下動しかしなくなるので、球同士の衝突はなくなってその後いくら振り続けても変化はない吸収状態である。有限温度の相では、充填率が大きいためにいつまで経ってもどこかで球同士の衝突が起こっていて水平方向の速度はゼロにならない。この変化は充填率とともにだんだんと起こるだけでなく、ある充填率で粉体温度が突然大きくなるという不連続な転移も起こりうるというのは驚くべきことだ。
しかも、不連続な転移になる条件は謎めいている。板の間隔hが大きめの(a)では、振幅Aを大きくしていくと不連続な転移から連続な転移に変化していくのに対し、板の間隔hが小さめの(b)では、逆に振幅Aを大きくしていくと連続な転移から不連続な転移に変化していく。
吸収状態転移は、充填率$${\phi}$$を小さくすると起こり、その転移が連続か不連続かは、板の間隔hと振幅Aとに依るようだが(a)と(b)で見事に逆の結果になり、どうやらhとAだけでは連続か不連続かの議論はできなさそうだ。

グラフ(a),(b)の違いを解明するため、球の振動の同期具合を調べてみる。粉体温度ゼロの相では、すべての球が板と同期して振動しているが、そうでない場合は、ばらばらに非周期的に動く場合と、同じ周期で動くが位相が異なる場合が考えられる。
振動の同期具合 s は次のように定義する。
$${s = \lim_{t_0,t→∞} \frac{1}{t} \int^{t_0+t}_{t_0}|⟨exp( iθt')⟩| dt' }$$
これも見た目難しそうだがどうということはない。
各粒子の振動の位相差をθで表し、半径1の円周上の偏角θのところに点があるとみなし、複数の点の位置の平均が中心からどれだけ離れているかを表している。偏角θが均等にばらばらになっていればs(点の重心位置)は0に近くなるが、偏角θが偏っていればsは原点から離れて大きな値になる。(動画としてはこれの4:04~がわかりやすい蔵本予想とは何か？～天才が証明したもの～　＠第８回日曜数学会 - ニコニコ (nicovideo.jp))
同期具合 s を、振幅 A と板の幅 h を変えてみて色で示したのがグラフ(c)である。

図2(c)　振幅Aと板の間隔hを変化させたときの同期具合s
図2(a), (b)で使ったh, Aの値が右上の黄緑～青の線

図2(c)の興味深い点は、同期具合の低い暗い領域が斜めに横切っているために、同期具合 s が A のみに依るとも h のみに依るとも切り分けることができないところだ。
そして、図2(a), (b)で使ったh, Aの値が右上の黄緑～青の線で示されているが、図2(a), (b)と見比べると同期具合の大きい明るい領域のところで不連続な転移になっている。

同期具合がなぜ転移の連続性に関わるかというと、次のように考えられる。
振動が同期した球同士が衝突するというのは、真横から衝突するので水平方向の速度成分しか変化しない上に完全弾性衝突でないので水平方向の速度は確実に落ちる。
同期してない球同士は、真横ではなくz方向に角度がついて衝突するので、z方向の速度が水平方向の速度に転換されて、水平方向の速度は増えることもある。
つまり、振動が同期してない状態より、同期した状態になると、水平方向のエネルギー損失は大きくなるので、急激に水平方向の速度が落ちるようになるためである。

より本質的な理解を得るために、このモデルを球1個1個シミュレーションするのではなく、ここまでに見えた3つの主要なメカニズム、①壁からのエネルギー補給、②水平方向の速度の減速、③縦方向の振動の同期化を、統計的な粗視化モデルで考える。
i番目とj番目の球の衝突で、球の衝突前の速度$${\bm{v}_i}$$と衝突後の速度$${\bm{v}'_i}$$を次のように考える。
$${\bm{v}’_i=\bm{v}_i+\frac{1+α}{2}(\bm{v}_{ij}\cdot\bm{σ}_{ij})\bm{σ}_{ij}+Δ_{ij}\bm{σ}_{ij}}$$
$${\bm{v}’_j=\bm{v}_i-\frac{1+α}{2}(\bm{v}_{ij}\cdot\bm{σ}_{ij})\bm{σ}_{ij}-Δ_{ij}\bm{σ}_{ij}}$$
$${\bm{v}_{ij}}$$はi番目とj番目の球の相対速度、$${\bm{σ}_{ij}}$$はその単位ベクトル。
α ははね返り係数(1が完全弾性衝突)で、2項目は水平方向のはね返りによる速度変化を表す。$${Δ_{ij}}$$＞0は速度注入で、3項目は垂直方向から水平方向へのエネルギーの移動(①)を表す。
板との衝突による減速(②)は、次のように粘性係数γを用いて表す。
$${\frac{d\bm{v}_i}{dt}=-γ\bm{v}_i}$$
振動の同期(③)については、前の衝突から一定の同期時間$${τ_s}$$以内なら$${Δ_{ij}=Δ}$$(>0)とし、$${τ_s}$$を超えたら$${Δ_{ij}=0}$$とする。つまり、ある程度の同期時間$${τ_s}$$経ったら勝手に振動は同期し、垂直方向から水平方向へのエネルギーの移動はなくなるものとする。特に、極限$${τ_s}$$ → ∞はいつまで経っても同期が起こらない系を表す。
このモデルで、粒子数 N = 20000、 τ∆/σ = 0.025、 α = 0.95、 τ γ = 0.01 でシミュレーションを行った結果が図3である。

図3　粗視化モデルによる理論値とシミュレーションの比較
左上は、同期がない場合。縦線は平均自由行程の議論から予測される臨界充填率。

球1個1個のシミュレーションで得られた値とよく一致する結果となった。
$${\hat{τ}}$$は平均自由行程で、$${τ_s}$$のそれに対する比で同期時間を表している。同期時間が有限の場合は不連続な転移となり、$${τ_s}$$ → ∞では転移は連続になっている。
また、同期時間が有限の場合は、球同士の散逸衝突によって吸収状態に達するが、同期時間が無限の場合は、粘性抗力(≒壁との衝突)によって吸収状態に達することも観測された。この2種類の転移の性質を有限サイズ解析で確認するため、次のようなエネルギー変化率Gを考える。
$${G(\phi,T)=\frac{\partial{T}}{\partial{t}}=\frac{ω(\phi, T)}{2}＜E'-E＞_{coll}-2γT}$$
1項目は衝突(collision)によるエネルギー変化、2項目は粘性抵抗によるエネルギーの損失を表す。この式を計算するにあたって、前の衝突から一定の同期時間$${τ_s}$$を超えるかどうかで$${Δ_{ij}}$$の値が変わるという部分で場合分けがあるために計算が煩雑になるが、$${Δ}$$を次のように近似することで計算は容易になる。
$${\bar{Δ}(\phi, T, τ_s)=\bar{Δ}_{ij}=Δ(1-\exp(-2ω(\phi,T)τ_s))}$$
括弧内は、同期してない球とポアソン分布的な衝突を仮定して、衝突する2つの粒子の少なくとも1つが過去に$${τ_s}$$より短い時間で衝突した確率を表す。前の衝突からの時間が$${τ_s}$$より長くなればΔの値はほぼゼロになることを場合分けなしで1つの式で書ける。
また、$${Δ_{ij}}$$は平均を取るとき次のように近似する。
$${＜Δ_{ij}g(\bm{v}_i,\bm{v}_j)＞\simeq\bar{Δ}＜g(\bm{v}_i,\bm{v}_j)＞}$$(gは任意の関数)
この近似は、粒子速度とその同期状態の間の相関を無視しているが、$${τ_s}$$ → ∞の極限では$${Δ_{ij}}$$の値は変わらなくなるので厳密になる。
球の速度分布をガウス分布と仮定すると、最終的にエネルギー変化率Gは次のようになる。
$${G(\phi,T)=\frac{ω(\phi, T)}{2}(m\bar{Δ}^2+α\bar{Δ}\sqrt{πmT}-T(1-α^2))-2γT}$$
ここから、グラフにプロットされる粉体温度の定常値は次のように求まる。
$${T_ss=(\frac{ε+\sqrt{ε^2+4mΔ^2(1-α)}}{2(1-α^2)})^2}$$
ただし、$${ε=αΔ\sqrt{πm}-4\frac{γ}{\bar{ω}}}$$, $${\bar{ω}=\frac{8\phiχ}{σ\sqrt{πm}}}$$, σはEnskog因子(導出の途中で用いる仮定)
この式は、球の充填率$${\phi}$$を変えても$${T_{ss}}$$は常に連続で不連続な転移は起こらない。これは、衝突の頻度を平衡と同じと仮定したことに起因する。この仮定は、平均自由時間が1/γに比べて小さい場合にのみ成り立つ。実際、ある臨界充填率$${\phi_c}$$以下では水平運動はゼロにならなければならないことを示すことができる。この$${\phi_c}$$において、粉体温度はゼロになり、したがって衝突後の粒子の水平方向の速度変化は∆に等しくなければならなず、衝突間の距離はΔ/γとなる。この距離と平均自由行程 l(φ) が等しくなるという条件で$${\phi_c}$$は求まる(その値が図3左上の縦線)。この値は、同期がない場合にシミュレー ションで観測された臨界充填率とよく一致する。同期がない場合($${\hat{τ}/τ_s}$$ = 0)は、エネルギー注入に関する近似が厳密になるため、挿入図に見られるような転移を除いて、非常にうまく一致する。同期時間が有限である場合、理論値は転移点をよく予測し、シミュレーションによる数値とよく一致する。立ち上がってすぐのところはずれが目立つが、これは近似によるものである(よい近似でも臨界点付近は収束が遅いために誤差が大きく出やすい)。

上記の通り、①壁からのエネルギー補給、②水平方向の速度の減速、③縦方向の振動の同期化を取り入れた理論は、シミュレーションとよく一致することから、吸収状態転移が不連続になるかどうかのメカニズムをかなり良く表していると言えるだろう。
そして、非常に単純なモデルから構築された理論であるため、ここに条件を加えることによって容易に、吸収状態転移の連続性の議論に関して広く応用できることが期待できる。

最後に
物理学者ジョークとして「牛を球体とみなす」というフレーズを聞いたことがあるかもしれない(柞刈湯葉氏の著書ではなく、その元ネタのほう球形の牛 - Wikipedia)。この1フレーズだけを見ると、物理学者が難しいこと考えすぎてわけのわからないことを言い出したかのように見えるが、本研究を見れば、その先に見える、球の衝突と密度変化だけから起こる相転移といった興味深い現象こそが物理学者の興味の対象だということがわかっていただけるだろう。