振動子と秩序変数

r(t)=0.171r(t) = 0.171
φ(t)=−0.653\varphi(t) = -0.653

秩序変数 r(t) の時間変化

r(t)r(t)
tt
r(t)r(t) 数値解 r=1−Kc/Kr = \sqrt{1 - K_c/K} 解析解

コーシー分布 g(ω)

g(ω; ω0, γ)g(\omega;\,\omega_0,\,\gamma)
ω\omega
g(ω)g(\omega)
NN
KK
ω0\omega_0
γ\gamma
δt\delta t

蔵本モデルの概要

NN 体の振動子からなる蔵本モデルは次の微分方程式で表される。

dθi(t)dt=ωi+KN∑j=1Nsin⁡(θj(t)−θi(t)),i=1,2,...,N(1)\frac{d\theta_i(t)}{dt} = \omega_i + \frac{K}{N}\sum_{j=1}^N \sin{( \theta_j(t) - \theta_i(t) ) }, \qquad i=1,2, ..., N \tag{1}

θi,ωi\theta_i, \omega_i は ii 番目の振動子の位相と自然振動数、KK は結合強度を意味する。

自然振動数 ωi\omega_i は確率密度関数 g(ω)g(\omega) に従う。 今回のシミュレーションでは、コーシー分布 g(ω;ω0,γ)g(\omega; \omega_0, \gamma) を考える。 コーシー分布の確率密度関数は、

g(ω;ω0,γ)=1πγ(ω−ω0)2+γ2g(\omega; \omega_0, \gamma)= \frac{1}{\pi} \frac{\gamma}{(\omega - \omega_0)^2 + \gamma^2}

で記述される。

次に秩序変数 zz を定義する。

z=r(t)e−1φ(t)=1N∑i=1Ne−1θi(t)r(t)=(∑i=1Ncos⁡θi(t))2+(∑i=1Nsin⁡θi(t))2φ(t)=arctan⁡(∑i=1Nsin⁡θi(t)∑i=1Ncos⁡θi(t))(2)\begin{aligned} z &= r(t)e^{\sqrt{-1}\varphi(t)} = \frac{1}{N}\sum_{i=1}^N e^{\sqrt{-1} \theta_i(t)} \tag{2} \\ r(t) &= \sqrt{\left(\sum_{i=1}^N \cos{\theta_i(t)}\right)^2 + \left( \sum_{i=1}^N \sin{\theta_i(t)} \right)^2} \\ \varphi(t) &= \arctan{\left( \frac{\sum_{i=1}^N \sin{\theta_i(t)}}{\sum_{i=1}^N \cos{\theta_i(t)}} \right)} \end{aligned}

r(t)r(t) の大きさを見ることで、振動子が同期しているかどうかが分かる。 r(t)r(t) が 1 に近いほど同期しており、0 に近いほど同期していない事を表している。

秩序変数を使って、(1) 式を次のように書き直すことができる。

dθi(t)dt=ωi+K×r(t)sin⁡(φ(t)−θi(t))(3)\frac{d\theta_i(t)}{dt} = \omega_i + K \times r(t) \sin{\left( \varphi(t) - \theta_i(t) \right) } \tag{3}

この導出は、方程式

r(t)e−1(φ(t)−θi(t))=1N∑j=1Ne−1(θj(t)−θi(t))r(t)e^{\sqrt{-1} (\varphi(t) - \theta_i(t))}=\frac{1}{N}\sum_{j=1}^N e^{\sqrt{-1}(\theta_j(t) - \theta_i(t))}

の虚部だけ考えると、

r(t)sin⁡(φ(t)−θi(t))=1N∑j=1Nsin⁡(θj(t)−θi(t))r(t)\sin{(\varphi(t) - \theta_i(t))} = \frac{1}{N}\sum_{j=1}^N \sin{(\theta_j(t) - \theta_i(t))}

が得られるので、これを (1) 式に代入すれば良い。

g(ω)g(\omega) がコーシー分布のときに、N→∞N\rightarrow \infty の r(t)r(t) の値は解析的に求まることが知られている。 結果だけ以下に書いておく。

r=1.0−KcK(Kc<K)Kc=2γ\begin{aligned} r &=\sqrt{1.0 - \frac{K_c}{K}} \qquad (K_c < K) \\ K_c &= 2\gamma \end{aligned}
導出

連続極限 N→∞N \to \infty では、振動子の分布を密度関数 f(θ,ω,t)f(\theta, \omega, t) で記述する。秩序変数は

r e−1φ=∫−∞∞ ⁣ ⁣∫02πe−1θ f(θ,ω,t) dθ dωr\,e^{\sqrt{-1}\varphi} = \int_{-\infty}^{\infty} \!\!\int_0^{2\pi} e^{\sqrt{-1}\theta}\, f(\theta,\omega,t)\,d\theta\,d\omega

と書ける。定常状態では振動子を2つのグループに分ける。

同期グループ (∣ω−ω0∣≤Kr|\omega - \omega_0| \le Kr): 位相が固定点に収束する。(3) 式の右辺を 0 とおくと

θ∗=φ+arcsin⁡ ⁣(ω−ω0Kr)\theta^* = \varphi + \arcsin\!\left(\frac{\omega - \omega_0}{Kr}\right)

非同期グループ (∣ω−ω0∣>Kr|\omega - \omega_0| > Kr): 位相がドリフトし続ける。定常分布は確率保存から

f(θ,ω)=C(ω)∣ω−ω0−Krsin⁡(θ−φ)∣f(\theta,\omega) = \frac{C(\omega)}{|\omega - \omega_0 - Kr\sin(\theta - \varphi)|}

秩序変数の自己無撞着方程式は、(2) 式の実部をとって

r=∫∣ω−ω0∣≤Krcos⁡(θ∗(ω)−φ) g(ω) dωr = \int_{|\omega - \omega_0| \le Kr} \cos(\theta^*(\omega) - \varphi)\, g(\omega)\,d\omega

ここで非同期グループの寄与は対称性から消え、同期グループのみが残る。 ω=ω0+Krsin⁡α\omega = \omega_0 + Kr\sin\alpha と変数変換すると

r=∫−π/2π/2cos⁡2 ⁣α    g(ω0+Krsin⁡α)  dαr = \int_{-\pi/2}^{\pi/2} \cos^2\!\alpha \;\; g(\omega_0 + Kr\sin\alpha)\;d\alpha

コーシー分布 g(ω;ω0,γ)g(\omega; \omega_0, \gamma) を代入し、Kr≫γKr \gg \gamma の近似で積分を評価すると

1=g(ω0)⋅π2  ⟹  Kc=2π g(ω0)=2γ1 = g(\omega_0)\cdot\frac{\pi}{2} \implies K_c = \frac{2}{\pi\, g(\omega_0)} = 2\gamma

r≠0r \neq 0 の解が存在する条件は K>Kc=2γK > K_c = 2\gamma であり、このとき

r=1−KcKr = \sqrt{1 - \frac{K_c}{K}}

が得られる。

数値シミュレーション

常微分方程式の数値解法の一つにオイラー法がある。 (1) 式に前進オイラー法を適用すると、

θit=θit−1+δt×{ωi+KN∑j=1Nsin⁡(θjt−1−θit−1)}\theta_{i}^{t} = \theta_i^{t-1} + \delta t \times \left\{ \omega_i + \frac{K}{N}\sum_{j=1}^N \sin{(\theta_j^{t-1} - \theta_i^{t-1})} \right\}

が得られる。 θit\theta_i^t は時刻 tt における ii 番目の振動子の位相、δt\delta t は時間幅である。これを愚直に実装すると、1ステップあたりの計算量は O(N2)O(N^2) になる。

(3) 式に前進オイラー法を適用すると、

θit=θit−1+δt×{ωi+K×r(t)sin⁡(φ(t)−θi(t))}\theta_{i}^{t} = \theta_i^{t-1} + \delta t \times \left\{ \omega_i + K \times r(t) \sin{\left( \varphi(t) - \theta_i(t) \right) } \right\}

が得られる。この場合、1ステップあたりの計算量は O(N)O(N) になる。

参考