高次相互作用を持った全結合イジング模型の相転移

概要

この記事では全結合イジング模型において、相互作用を一般の多体相互作用に拡張した模型の相転移について考えます。相互作用は一様で強磁性的なものを仮定します。このモデルについては検索してもあまりヒットしないのですが、モデルの性質自体はすでによく知られており^[1]、相互作用の次数を$pp$として$p>2p>2$で一次転移を示すことが分かっています。この記事ではこのことを確認します。

ちなみにどうしてこのような模型に興味を持ったかというと、自分はJijに入社した当初、高次制約なし二値最適化問題^[2]に対してシミュレーテッドアニーリングを行うソルバーを開発^[3]していたのですが、その際にソルバーが正しい挙動をしているか確かめるために厳密解の分かる模型が欲しくなり、ここで議論するような模型を考えたのがきっかけです。

モデル

一般の場合

高次相互作用を持った全結合イジング模型は以下で定義されます。

$$Ep=−1(N−1p−1)∑i1=1N∑i2=i1+1N⋯∑ip=ip−1+1Nσi1σi2⋯σip E_{p}=-\frac{1}{\binom{N-1}{p-1}}\sum^{N}_{i_1=1}\sum^{N}_{i_2=i_1+1}\cdots\sum^{N}_{i_p=i_{p-1}+1}\sigma_{i_1}\sigma_{i_2}\cdots\sigma_{i_p}$$

ここで、$NN$は全スピン数、$σij∈{−1,+1}\sigma_{i_j}\in\{-1, +1\}$はスピン変数、$pp$は相互作用の次数であり、

$$(N−1p−1)=(N−1)(N−2)⋯(N−p+1)(p−1)(p−2)⋯1 \binom{N-1}{p-1}=\frac{(N-1)(N-2)\cdots(N-p+1)}{(p-1)(p-2)\cdots 1}$$

は相互作用の規格化定数です。この規格化定数は自由エネルギーがスピン数$NN$に対して示量的になるように決めればよく、いくらかの任意性があります。ここでは以降の計算が簡単になるように$(N−1p−1)\binom{N-1}{p-1}$としており、例えば$Np−1N^{p-1}$としても$NN$が大きい極限では温度が適当にスケールされるだけで基本的な性質は変わりません。ちなみに$(N−1p−1)\binom{N-1}{p-1}$はある一つのスピンと相互作用しているスピンの総数になっています。

3体相互作用の場合

この記事では主に$p=3p=3$のモデルを考えます。ハミルトニアンは以下のとおりです。

$$Ep=3=−2(N−1)(N−2)∑i=1N∑j=i+1N∑k=j+1Nσiσjσk E_{p=3}=-\frac{2}{(N-1)(N-2)}\sum^{N}_{i=1}\sum^{N}_{j=i+1}\sum^{N}_{k=j+1}\sigma_{i}\sigma_{j}\sigma_{k}$$

ここで、便利のために以下の全スピンの和を定義しておきます。

$$Ψ=∑i=1Nσi \Psi=\sum^{N}_{i=1}\sigma_{i}$$

この量の３乗を計算すると、

$$Ψ3=∑i=1N∑j=1N∑k=1Nσiσjσk=∑i=1Nσi3+3∑i=1N∑j=1i≠jNσi2σj+∑i=1N∑j=1N∑k=1i≠j≠kNσiσjσk=∑i=1Nσi3+3(∑i=1N∑j=1Nσi2σj−∑i=1Nσi3)+6∑i=1N∑j=i+1N∑k=j+1Nσiσjσk=Ψ+3(NΨ−Ψ)+6∑i=1N∑j=i+1N∑k=j+1Nσiσjσk=(3N−2)Ψ+6∑i=1N∑j=i+1N∑k=j+1Nσiσjσk \Psi^3=\sum^{N}_{i=1}\sum^{N}_{j=1}\sum^{N}_{k=1}\sigma_{i}\sigma_{j}\sigma_{k}\\ =\sum^{N}_{i=1}\sigma^3_{i} + 3\sum^{N}_{i=1}\sum^{N}_{\substack{j=1 \\ i \neq j}}\sigma^2_{i}\sigma_{j}+\sum^{N}_{i=1}\sum^{N}_{j=1}\sum^{N}_{\substack{k=1 \\ i \neq j \neq k}}\sigma_{i}\sigma_{j}\sigma_{k}\\ =\sum^{N}_{i=1}\sigma^3_{i}+3\left(\sum^{N}_{i=1}\sum^{N}_{j=1}\sigma^2_{i}\sigma_{j} - \sum^{N}_{i=1}\sigma^3_{i}\right) + 6\sum^{N}_{i=1}\sum^{N}_{j=i+1}\sum^{N}_{k=j+1}\sigma_{i}\sigma_{j}\sigma_{k}\\ =\Psi+3(N\Psi - \Psi)+6\sum^{N}_{i=1}\sum^{N}_{j=i+1}\sum^{N}_{k=j+1}\sigma_{i}\sigma_{j}\sigma_{k}\\ =(3N-2)\Psi+6\sum^{N}_{i=1}\sum^{N}_{j=i+1}\sum^{N}_{k=j+1}\sigma_{i}\sigma_{j}\sigma_{k}$$

となります。ここで$σi2=1\sigma^2_{i}=1$を使いました。上記関係式を用いるとハミルトニアンを$Ψ\Psi$の関数として表現できて、

$$Ep=3=−13(N−1)(N−2){Ψ3−(3N−2)Ψ} E_{p=3}=-\frac{1}{3(N-1)(N-2)}\left\{\Psi^3-(3N-2)\Psi\right\}$$

となります。

分配関数の計算

上の式を分配関数の定義に代入すると^[4]、

$$Z=∑σ1=±1∑σ2=±1⋯∑σN=±1exp⁡(−βEp=3)=∑σ1=±1∑σ2=±1⋯∑σN=±1exp⁡[β3(N−1)(N−2){Ψ3−(3N−2)Ψ}] Z = \sum_{\sigma_{1}=\pm 1}\sum_{\sigma_{2}=\pm 1}\cdots\sum_{\sigma_{N}=\pm 1}\exp\left(-\beta E_{p=3}\right)\\ =\sum_{\sigma_{1}=\pm 1}\sum_{\sigma_{2}=\pm 1}\cdots\sum_{\sigma_{N}=\pm 1}\exp\left[\frac{\beta}{3(N-1)(N-2)}\{\Psi^3-(3N-2)\Psi\}\right]$$

となります。ここで全てのスピンについての和、$∑σ1=±1⋯∑σN=±1\sum_{\sigma_{1}=\pm1}\cdots\sum_{\sigma_{N}=\pm1}$について、個別に和を取るのではなく、全スピン$Ψ\Psi$で和をとるように書き換えることができます。

$$Z=∑Ψ=−NNW(Ψ)exp⁡[β3(N−1)(N−2){Ψ3−(3N−2)Ψ}] Z=\sum^{N}_{\Psi=-N}W(\Psi)\exp\left[\frac{\beta}{3(N-1)(N-2)}\{\Psi^3-(3N-2)\Psi\}\right]$$

ここで新しく現れた$W(Ψ)W(\Psi)$は全スピンの値が$Ψ\Psi$となる場合の数であり以下のようになります。

$$W(Ψ)={N!(N+Ψ2)!(N−Ψ2)!N±Ψが0以上の偶数の時0上記以外 W(\Psi)=\left\{ \begin{array}{ll} \frac{N!}{\left(\frac{N+\Psi}{2}\right)!\left(\frac{N-\Psi}{2}\right)!} & N\pm\Psiが0以上の偶数の時 \\ 0 & 上記以外 \end{array} \right.$$

$NN$が十分大きいとしてスターリングの公式$log⁡N!≃Nlog⁡N−N\log N!\simeq N\log N - N$を適用すると、$N±ΨN\pm\Psi$が0以上の偶数の時、

$$log⁡W(Ψ)=log⁡N!−log⁡(N+Ψ2)!−log⁡(N−Ψ2)!≃Nlog⁡N−N−(N+Ψ2log⁡N+Ψ2−N+Ψ2)−(N−Ψ2log⁡N−Ψ2−N−Ψ2)=−N2{(1+ΨN)log⁡(1+ΨN)+(1−ΨN)log⁡(1−ΨN)−2log⁡2} \log W(\Psi)=\log N! - \log\left(\frac{N+\Psi}{2}\right)!-\log\left(\frac{N-\Psi}{2}\right)!\\ \simeq N\log N - N - \left(\frac{N+\Psi}{2}\log\frac{N+\Psi}{2} - \frac{N+\Psi}{2}\right) - \left(\frac{N-\Psi}{2}\log\frac{N-\Psi}{2} - \frac{N-\Psi}{2}\right)\\ =-\frac{N}{2}\left\{\left(1+\frac{\Psi}{N}\right)\log\left(1+\frac{\Psi}{N}\right)+\left(1-\frac{\Psi}{N}\right)\log\left(1-\frac{\Psi}{N}\right)-2\log 2\right\}$$

となります。ここで新たにスピンあたりの磁化

$$ψ=ΨN \psi=\frac{\Psi}{N}$$

および、

$$σ(ψ)=log⁡2−12{(1+ψ)log⁡(1+ψ)+(1−ψ)log⁡(1−ψ)} \sigma(\psi)=\log2 - \frac{1}{2}\left\{(1+\psi)\log(1+\psi)+(1-\psi)\log(1-\psi)\right\}$$

を定義すると、

$$log⁡W(Ψ)≃Nσ(ψ) \log W(\Psi) \simeq N\sigma(\psi)$$

となるので、これの$log⁡\log$を外して、

$$W(Ψ)≃exp⁡(Nσ(ψ)) W(\Psi)\simeq \exp\left(N\sigma(\psi)\right)$$

という関係式を分配関数の表式に代入して$NN$が十分大きいとすると、

$$Z=∑Ψ=−NNW(Ψ)exp⁡[β3(N−1)(N−2){Ψ3−(3N−2)Ψ}]≃N2∫−11exp⁡[β3(N−1)(N−2){Ψ3−(3N−2)Ψ}+Nσ(ψ)]dψ=∫−11exp⁡N[β3N(N−1)(N−2){Ψ3−(3N−2)Ψ}+σ(ψ)+1Nlog⁡N2]dψ≃∫−11exp⁡N[βψ33+σ(ψ)]dψ Z=\sum^{N}_{\Psi=-N}W(\Psi)\exp\left[\frac{\beta}{3(N-1)(N-2)}\{\Psi^3-(3N-2)\Psi\}\right]\\ \simeq \frac{N}{2}\int^{1}_{-1}\exp\left[\frac{\beta}{3(N-1)(N-2)}\{\Psi^3-(3N-2)\Psi\}+N\sigma(\psi)\right]d\psi\\ =\int^{1}_{-1}\exp N\left[\frac{\beta}{3N(N-1)(N-2)}\{\Psi^3-(3N-2)\Psi\}+\sigma(\psi)+\frac{1}{N}\log\frac{N}{2}\right]d\psi\\ \simeq\int^{1}_{-1}\exp N\left[\frac{\beta\psi^3}{3}+\sigma(\psi)\right]d\psi$$

と評価することができます。ここで、$NN$が十分大きいとして

$$∑Ψ=−NN→N2∫−11dψ \sum^{N}_{\Psi=-N}\rightarrow \frac{N}{2}\int^{1}_{-1}d\psi$$

と、和を積分で書き換えています。

さて、次はこの積分を計算します。$σ(ψ)\sigma(\psi)$の定義を代入してそのまま計算することは困難なので、ラプラスの方法を使って評価します。つまり、上記積分値に寄与するのは指数関数の中が最大値を取るときなので、

$$f(ψ)=βψ33+σ(ψ) f(\psi)=\frac{\beta \psi^3}{3}+\sigma(\psi)$$

と、指数関数の中で$NN$に比例する項をまとめると、

$$Z≃∫−11exp⁡Nf(ψ)∼exp⁡(Nf(ψ∗))2πN∣f′′(ψ∗)∣ Z\simeq\int^{1}_{-1}\exp Nf(\psi)\sim \exp\left(Nf(\psi^*)\right)\sqrt{\frac{2\pi}{N|f''(\psi^*)|}}$$

と形式的に書くことができます。ここで、

$$f(ψ∗)=max⁡−1≤ψ≤1{f(ψ)} f(\psi^*)=\max_{-1\leq \psi \leq 1}\{f(\psi)\}$$

であり、$ψ∗\psi^*$は関数$f(ψ)f(\psi)$の最大値を与える量です(ただし、この表記では最大値を与える$ψ\psi$が唯一であると仮定しています)。

磁化の計算

分配関数の計算過程で出てきた$ψ∗\psi^*$という量が実はこの系の磁化に対応する量であることがわかります^[5]。
よって、$f(ψ)f(\psi)$の最大値を与える$ψ\psi$の$β\beta$依存性が重要になります。それが分かれば任意の温度$T=1/βT=1/\beta$での磁化の値が分かることになります。
関数$f(ψ)f(\psi)$の定義はその定義域を$−1≤ψ≤1-1\leq\psi\leq1$として、

$$f(ψ)=βψ33+log⁡2−12{(1+ψ)log⁡(1+ψ)+(1−ψ)log⁡(1−ψ)} f(\psi)=\frac{\beta \psi^3}{3}+\log2 - \frac{1}{2}\left\{(1+\psi)\log(1+\psi)+(1-\psi)\log(1-\psi)\right\}$$

でした。最大値を与える$ψ\psi$の候補は以下の方程式で与えられます。

$$f′(ψ)=βψ2+12log⁡1−ψ1+ψ=0⇒ψ=tanh⁡(βψ2) f'(\psi)=\beta \psi^2+\frac{1}{2}\log\frac{1-\psi}{1+\psi}=0\\ \Rightarrow \psi=\tanh(\beta \psi^2)$$

これは平均場近似による自己整合方程式とよく似ています。実際、ハミルトニアン$Ep=3E_{p=3}$において$σi→m+δσi\sigma_{i}\rightarrow m+\delta\sigma_{i}$という置き換えを行い$δσi\delta\sigma_{i}$の二次以上の項を落とすと上記方程式が得られます。
通常の２体相互作用を持ったイジング模型における自己整合方程式は$ψ=tanh⁡(βψ)\psi=\tanh(\beta \psi)$となるので、$ψ\psi$の次数が一つ大きくなったことが違いとなっています。

ここから先は実際に$f(ψ)f(\psi)$のグラフを書いてみて磁化$m=ψ∗m=\psi^*$の挙動を確認していきます。
先程、$f(ψ)f(\psi)$の定義域は$−1≤ψ≤1-1\leq\psi\leq1$と述べましたが、今考えている系は3体相互作用を持っているため、スピン反転対称性$σ→−σ\sigma\rightarrow -\sigma$を持たず、低温では明らかに磁化は正の値を取ることが分かります。よって、以下では$f(ψ)f(\psi)$の定義域を$0≤ψ≤10\leq\psi\leq1$に制限します。
以下に$β=1.5,1.7,1.8,2.3\beta=1.5,1.7,1.8,2.3$という典型的な逆温度での関数形を示しました。また、念のためソースコードも載せておきました。

import matplotlib.pyplot as plt
import numpy as np

beta=2.3
x = np.arange(0, 1.1, 0.1)
y1 = x
y2 = np.tanh(beta*x*x)
plt.plot(x, y1, label=r"$y=\psi$", color="black")
plt.plot(x, y2, label=r"$y=\tanh(\beta\psi^2),\beta=$"+f"{beta}", color="red")
plt.xlabel(r'$\psi$')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.legend()

ここから、$ψ=tanh⁡(βψ2)\psi=\tanh(\beta \psi^2)$という方程式は$β≲1.7\beta\lesssim 1.7$で$ψ=0\psi=0$、$β∼1.7\beta\sim 1.7$で$ψ=0\psi=0$と$ψ∼0.8\psi\sim 0.8$という２つの解を持ち、$β≳1.7\beta\gtrsim 1.7$では３つの解を持つことが分かります。$β≳1.7\beta\gtrsim 1.7$における$ψ=0\psi=0$と$β\beta$の増大にともなって単調減少する２つの解は、非物理的なためこれ以降は考察しません。
さて、上記の考察からこの系の磁化は

• $β≲1.7\beta\lesssim 1.7$ではゼロ
• $β≳1.7\beta\gtrsim 1.7$で$ψ∼0.8\psi\sim 0.8$という値に不連続に変化し、以後$β\beta$の増大にともなって増加していく

という振る舞いをすることが分かりました。したがってこの系は一次転移を示すことが示唆されます。
しかし、ここで一つ注意しなければいけない点があります。というのも$ψ=tanh⁡(βψ2)\psi=\tanh(\beta \psi^2)$という方程式の解は、$f(ψ)f(\psi)$の最大値を与える単なる候補であるということです。そこで、先程と同じ$β\beta$に対して今度は$f(ψ)f(\psi)$のグラフを書いてみます。

import matplotlib.pyplot as plt
import numpy as np

beta=2.3
x = np.arange(0, 1, 0.001)
y = beta*x*x*x/3+np.log(2)-0.5*((1+x)*np.log(1+x)+(1-x)*np.log(1-x))
plt.plot(x, y, label=r"$y=f(\psi),\beta=$"+f"{beta}", color="black")
plt.xlabel(r'$\psi$')
plt.xlim([0, 1])
plt.ylim([0.4, 0.8])
plt.legend()

結果は一目瞭然です。自己整合方程式は$β∼1.8\beta\sim 1.8$では明らかに有限の解を持っていましたが、もとの関数$f(ψ)f(\psi)$は$β=1.8\beta=1.8$では$ψ=0\psi=0$で最大値を取っていることが分かります。つまり、このときの磁化はまだゼロだということです。一方で$β=2.3\beta=2.3$では$f(ψ)f(\psi)$の最大値は有限の$ψ\psi$が与えていることが分かるので、この系が一次転移を示すことは間違いなさそうです。先程の自己整合方程式の解による議論は転移温度$Tc=1/βcT_{\text{c}}=1/\beta_{\text{c}}$を過大評価していたことが分かりました^[6]。

最後に、磁化$mm$の温度$TT$依存性のグラフを書いてみます。

import matplotlib.pyplot as plt
import numpy as np

temperature_range = np.arange(0.1, 1.0, 0.001)
magnetization = []
for T in temperature_range:
    beta = 1/T
    x = np.arange(0, 1, 0.001)
    y = beta*x*x*x/3+np.log(2)-0.5*((1+x)*np.log(1+x)+(1-x)*np.log(1-x))
    magnetization.append(x[np.argmax(y)])
plt.scatter(temperature_range, magnetization, label="magnetization", color="black",s=10)
plt.xlabel(r'$T$')
plt.ylabel(r'$m$')
plt.legend()

およそ$T=0.5T=0.5$くらいで磁化が不連続に飛んでいることが分かります。確かに１次転移です。

一般の$pp$体相互作用の場合の磁化

ここまでは$p=3p=3$の模型での議論でしたが、一般の$pp$体相互作用の模型でもほとんど同様の議論をそのまま行うことができます。その結果、ハミルトニアン

$$Ep=−1(N−1p−1)∑i1=1N∑i2=i1+1N⋯∑ip=ip−1+1Nσi1σi2⋯σip E_{p}=-\frac{1}{\binom{N-1}{p-1}}\sum^{N}_{i_1=1}\sum^{N}_{i_2=i_1+1}\cdots\sum^{N}_{i_p=i_{p-1}+1}\sigma_{i_1}\sigma_{i_2}\cdots\sigma_{i_p}$$

で定まる系の磁化$mm$の逆温度依存性は

$$f(ψ)=βψpp+log⁡2−12{(1+ψ)log⁡(1+ψ)+(1−ψ)log⁡(1−ψ)} f(\psi)=\frac{\beta \psi^p}{p}+\log2 - \frac{1}{2}\left\{(1+\psi)\log(1+\psi)+(1-\psi)\log(1-\psi)\right\}$$

という関数の最大値を与える$ψ\psi$となることが分かります。また、自己整合方程式は

$$ψ=tanh⁡(βψp−1) \psi=\tanh(\beta \psi^{p-1})$$

となります。ただし、$p=3p=3$の時と同様に磁化の値として採用する自己整合方程式の解は慎重に選ぶ必要があります。

最後に、グラフが重なってしまって少し見にくいですが、いくつかの$pp$についての磁化$mm$のグラフを載せておきます。

import matplotlib.pyplot as plt
import numpy as np

p_list=[2, 3, 4]
temperature_range = np.arange(0.1, 1.0, 0.001)
magnetization = []
for p in p_list:
    temp_mag = []
    for T in temperature_range:
        beta = 1/T
        x = np.arange(0, 1, 0.001)
        y = beta*x**p/p+np.log(2)-0.5*((1+x)*np.log(1+x)+(1-x)*np.log(1-x))
        temp_mag.append(x[np.argmax(y)])
    magnetization.append(temp_mag)
for mag, p in zip(magnetization, p_list):
    plt.scatter(temperature_range, mag, label=r"magnetization for $p=$"+f"{p}",s=10)
plt.xlabel(r'$T$')
plt.ylabel(r'$m$')
plt.legend()

まとめ

一様かつ強磁性的な$pp$相互作用を持った全結合イジング模型の相転移について解析しました。その結果、$p>2p>2$では一次転移を示すことが分かりました^[7]。

今後について

モンテカルロシミュレーションでここで行った解析結果が再現されるか確認してみたいと思います^[8]。

2023/12/22追記
モンテカルロシミュレーションの記事を書きました。

最後に

＼Rustエンジニア・数理最適化エンジニア募集中！／
株式会社Jijでは、数学や物理学のバックグラウンドを活かし、量子計算と数理最適化のフロンティアで活躍するRustエンジニア、数理最適化エンジニアを募集しています！
詳細は下記のリンクからご覧ください。皆さんのご応募をお待ちしております！
Rustエンジニア: https://open.talentio.com/r/1/c/j-ij.com/pages/51062
数理最適化エンジニア: https://open.talentio.com/r/1/c/j-ij.com/pages/75132

脚注

1. 例えば、このノートに記載があります。

2. HUBOまたはHOBO(Higher-Order Unconstraint Binary Optimization)、PUBO(Polynomial Unconstraint Binary Optimizaion)などと呼ばれます。

3. OpenJijという名のOSSとして公開されています。

4. 田崎 統計力学IIより。分配関数の計算はこの本の第11章に準拠、というかほぼそのままです。

5. ハミルトニアン$Ep=3E_{p=3}$に磁場の項を入れて本節と同様の議論を展開すれば自由エネルギーが磁場依存性を持つので偏微分して直接確かめることができます。

6. これはつまり、$σi→m+δσi\sigma_{i}\rightarrow m+\delta\sigma_{i}$と置いて$δσi\delta\sigma_{i}$の二次以上の項を無視するような典型的な平均場近似を行うと、$p>2p>2$では転移温度を過大評価してしまうことを示しています。

7. 正確には$p=3,4,5p=3,4,5$に場合について、数値的に磁化の値をプロットして、一次転移の振る舞いを確認したということになります。

8. そもそも、自分が実装したモンテカルロシミュレーションのコードが正しいか確認するためにこのような解析を行ったのでした。