$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$
この記事は物工/計数 Advent Calendar 2023 17日目の記事です。
応物アドカレなので専門分野の記事もあった方がいいと思い、今回は物理系の人たちには馴染み深いイジングモデルを題材にします1。全体的に物工B3の授業内容（特に統計力学第二）と関わりの深い内容なので、物工の雰囲気を知りたい人にはぜひ読んでいただけると幸いです。
イジングモデルは磁性体の相転移を記述するモデルとして導入されました。ハミルトニアンは
H = -\sum_{<i,j>} J_{ij}S_iS_j-\sum_i h_iS_i
です。ここで$S_i$は$d$次元立方格子上で定義されたスピン変数（値は$\pm1$）であり、第1項の和は最近接サイトについて取ります。以下では、相互作用および磁場の等方性を仮定します（$J_{ij}=J, h_i=h$）。
横磁場イジングモデル
早速ですが、本稿の主題である横磁場イジングモデルについて説明します。横磁場イジングモデルのハミルトニアンは
H = -\sum_{<i,j>} J\sigma_i^z\sigma_j^z -\Gamma \sum_i \sigma_i^x
です。先ほどと違って磁場が$x$方向にかかっています。対応してスピン変数は$\sigma^x$となり、これは$\sigma^z$と交換しません。
定性的に横磁場イジングモデルの基底状態について考えてみます。$\Gamma$が（$J$に比べて）小さければ通常のイジングモデルと同じ振る舞いを示すはずなので、低温ではスピンが$z$方向に揃います。横磁場を大きくしていくとx方向にスピンが揃うので、絶対零度であってもz方向の対称性が回復する相転移が起きると期待されます。私たちが日常的に目にする相転移現象はエネルギーとエントロピーの拮抗という描像で理解されますが、この場合の相転移は演算子の非可換性に起因するものであり、（絶対零度で起きていることからも分かるように）先述の相転移とはメカニズムのレベルで異なるものです。演算子の非可換性は量子論に由来するものなので、このような相転移は特に「量子相転移」と呼ばれ、物性物理のホットトピックの一つとなっています。2
以下では、特に1次元横磁場イジングモデルの基底状態について高橋 和孝、西森 秀稔 著『相転移・臨界現象とくりこみ群』の議論を参考に考察を行い、最後に数値実験で確認します。実は1次元横磁場イジングモデルには厳密解が存在し、自由エネルギー密度は熱力学極限で
f = -\frac{1}{\beta}\int_{-\pi}^{\pi} \frac{dk}{2\pi}\ln(2\cosh\beta\sqrt{\Gamma^2+J^2+2\Gamma J \cos k})
となります。絶対零度極限$\beta\rightarrow \infty$を取り、磁化率を計算すると
\chi = -\frac{\partial^{2}f}{\partial^{2}\Gamma}=\int_{-\pi}^{\pi} \frac{dk}{2\pi} \frac{\Gamma^2 \sin^2 k}{(\Gamma^2+J^2+2\Gamma J \cos k)^{3/2}}
ですが、これは$\Gamma=J$のとき発散し、ここが相転移点であると分かります。
しかし以下で見るように、相転移点を知るだけであれば厳密解は必要ありません。ここでは、$d$次元横磁場イジングモデルが$d+1$次元イジングモデルにマップできることを示した上で、具体例として2次元イジングモデルのKramers-Wannier双対性を使った臨界温度の議論から元の1次元横磁場イジングモデルの量子相転移点が計算できることを示します。
ちなみに、物工の3Aの演習では横磁場イジングの基底状態の厳密解を求めさせられる回がある（少なくとも今年はあった）ので、2年生の方は楽しみにしておいてください！3
横磁場イジングと古典イジングの対応
※以下、横磁場イジングと区別するため、通常のイジングモデルを「古典イジングモデル」と呼ぶことがあります（演算子の非可換性が問題にならないという意味で）。
横磁場イジングモデルと通常のイジングモデルに何らかの対応関係があると仮定してみます。すなわち、$H$を横磁場イジングのハミルトニアン、$H_c(\boldsymbol{\sigma})$を古典イジングモデルのハミルトニアンとして、
Z=\text{Tr}[e^{-\beta H}] = \sum_{\boldsymbol{\sigma}}e^{-\beta H_c(\boldsymbol{\sigma})}
という関係が成り立つとしましょう。目標は、古典イジングモデルのハミルトニアンの結合定数$J,h$を横磁場イジングモデルの結合定数$J,\Gamma$で表すことです。
この式の左辺を直接評価するのが難しいのは、横磁場イジングモデルのハミルトニアンが$\sigma^z$を含む項$H_z$と$\sigma^x$を含む項$H_x$の和として書かれているためです。両者は交換しないので、指数関数についての等式
e^{-\beta H} = e^{-\beta H_z}e^{-\beta H_x} \;\;\;(?)
は成立しません。$e^{-\beta H}$を素直に計算するには系のサイズを$N$として疎行列ではない$2^N\times2^N$行列の指数関数を計算せねばならず、大変です。
この手の非可換な行列の和の指数関数を評価するのに便利な公式として、鈴木-Trotter公式
e^{-\beta(H_z+H_x)}=\lim_{M\rightarrow\infty}(e^{-\frac{\beta}{M}H_z}e^{-\frac{\beta}{M}H_x})^M
があります。この式は
( 右辺 )=\lim_{M\rightarrow\infty}(1-\frac{H_z}{M}-\frac{H_x}{M}+O(M^{-2}))^M = ( 左辺 )
と考えると納得できます。この公式により、
\text{Tr}[e^{-H}] = \lim_{M\rightarrow\infty}\sum_{\boldsymbol{\sigma}_1} \bra{\boldsymbol{\sigma}_1}(e^{-\frac{\beta}{M}H_z}e^{-\frac{\beta}{M}H_x})^M\ket{\boldsymbol{\sigma}_1}
を得ます。ここで$\text{Tr}$をとる基底として$\otimes \sigma_i^z$の固有状態を選んでおきます。
さらに計算を進めるため、$M$個ある$e^{-\frac{\beta}{M}H_z}e^{-\frac{\beta}{M}H_x}$の間に完全系
\sum_{\boldsymbol{\sigma_l}} \ket{\boldsymbol{\sigma_l}}\bra{\boldsymbol{\sigma_l}} = 1\;\;\;(l=2,3,\cdots ,M)
を挿入します。すると
\text{Tr}[e^{-H}] = \lim_{M\rightarrow\infty}\sum_{\boldsymbol{\sigma}_1, \boldsymbol{\sigma}_2,\cdots, \boldsymbol{\sigma}_M} \bra{\boldsymbol{\sigma}_1}e^{-\frac{\beta}{M}H_z}e^{-\frac{\beta}{M}H_x}\ket{\boldsymbol{\sigma}_2}
\bra{\boldsymbol{\sigma}_2}e^{-\frac{\beta}{M}H_z}e^{-\frac{\beta}{M}H_x}\ket{\boldsymbol{\sigma}_3}
\cdots
\bra{\boldsymbol{\sigma}_M}e^{-\frac{\beta}{M}H_z}e^{-\frac{\beta}{M}H_x}\ket{\boldsymbol{\sigma}_1}
ですが、$\ket{\sigma_l}$は$H_z$の固有状態なので、
\text{Tr}[e^{-H}] =
\lim_{M\rightarrow\infty}\sum_{\boldsymbol{\sigma}_1, \boldsymbol{\sigma}_2,\cdots, \boldsymbol{\sigma}_M}
\bra{\boldsymbol{\sigma}_1}e^{-\frac{\beta}{M}H_x}\ket{\boldsymbol{\sigma}_2}e^{-\frac{\beta}{M}H_z(\boldsymbol{\sigma}_1)}
\bra{\boldsymbol{\sigma}_2}e^{-\frac{\beta}{M}H_x}\ket{\boldsymbol{\sigma}_3}e^{-\frac{\beta}{M}H_z(\boldsymbol{\sigma}_2)}
\cdots
\bra{\boldsymbol{\sigma}_M}e^{-\frac{\beta}{M}H_x}\ket{\boldsymbol{\sigma}_1}e^{-\frac{\beta}{M}H_z(\boldsymbol{\sigma}_M)}
となります。明らかに
e^{-\frac{\beta}{M}H_z(\boldsymbol{\sigma}_l)} = e^{\frac{\beta J}{M}\sum_{i}\sigma_{i,l}\sigma_{i+1,l}}
です。次に$H_x$の行列要素を考えます。各スピンが独立であることから$i$番目のスピンについてのみ考えると、等式4
\bra{\sigma_{l}}e^{h \sigma^x_i }\ket{\sigma_{l+1}} =\bra{\sigma_{l}}(\cosh h+ \sinh h\sigma^x_i)\ket{\sigma_{l+1}}
=A(h)e^{K(h)\sigma_{i,l}\sigma_{i,l+1}}
が成り立つような$A(h),K(h)$があれば古典イジング型相互作用に書き直せます（ここで$h=\frac{\beta \Gamma}{M}$）。そのような$A(h),K(h)$を探すと5、
K(h) = -\frac{1}{2}\log \tanh(h), \; A(h)^2 = \frac{1}{2}\sinh (2h)
を得ます。以上により、横磁場イジングモデルの分配関数は
Z = \lim_{M\rightarrow\infty}A(\frac{\beta \Gamma}{M})^{MN}\sum_{\boldsymbol{\sigma}}e^{\sum_{i,l}(\frac{\beta J}{M}\sigma_{i,l}\sigma_{i+1,l}+K(\frac{\beta \Gamma}{M})\sigma_{i,l}\sigma_{i,l+1})}
と書き表すことができます。これは非等方的な（$x$方向と$y$方向で結合定数が異なる）2次元イジングモデルの分配関数の形をしているので、両者は対応関係にあることが確かめられました。
2次元古典イジングモデルの解析
2次元正方格子イジングモデルは、Kramers-Wannier双対性を用いることで厳密な臨界温度が（厳密解を経由せずに）計算できる例として知られています。ここではその議論を簡単に紹介した後、2次元イジングモデルの臨界温度を決定する方程式に前節で求めた結合定数の変換則を代入することで、1次元横磁場イジングモデルの量子相転移点を求めます。
零磁場での2次元イジングモデルの分配関数は
Z = \sum_{\boldsymbol{\sigma}}e^{\beta J \sum_{<i,j>} \sigma^z_i\sigma^z_j} = \sum_{\boldsymbol{\sigma}}\prod_{<i,j>}(\cosh \beta J + \sigma^z_i\sigma^z_j \sinh \beta J)
= (\cosh \beta J)^{2N} \sum_{\boldsymbol{\sigma}}\prod_{<i,j>} (1 + \sigma^z_i\sigma^z_j \tanh \beta J)
と書けます。なお、簡単のため系は等方的であるとしました（横磁場イジングとの対応を考える上では非等方的な場合に関心がありますが、本質的には変わりません）。ここで右辺を$u=\tanh(\beta J)$で展開してみましょう。これは$\beta$が小さいとき良い展開になっているので、高温展開と呼ばれます。パウリ行列の性質$\text{Tr}[ \sigma^z] = 0$および$\text{Tr}[ {\sigma^z}^2] = \text{Tr}[ 1] = 2$から、和に寄与するのは全てのサイト$i$について$\sigma^z_i$のべきが偶数乗となっているときであり、これは和に寄与するサイトを格子に沿って結ぶと閉曲線を描くことと同値です。従って、$K=\beta J$、$n$本の辺からなる閉グラフの数を$\Omega_n$とすると
Z = (\cosh K)^{2N} 2^N \sum_{n} \Omega_n u^n
が成り立ちます。
一方、同じ分配関数を次のように表現することもできます。すなわち、格子全体で隣同士スピンが揃っていないペアが$\nu$個あるとすると、残っているペアは$2N-\nu$個ありますから、そのような配置パターンの総数を$\omega_\nu$として
Z=\sum_{\nu}\omega_\nu e^{\beta J (2N-\nu-\nu)} = e^{2NK}\sum_{\nu}\omega_\nu e^{-2K\nu}
です。これは$K$（つまり$\beta$）が大きいとき良い展開になっているので、低温展開と呼ばれます。
ここで2つの展開の対応関係について考えてみます。元の正方格子の全ての正方形の中心にサイトを作ると再び正方格子が得られますが（双対格子と言います）、そのとき高温展開での閉グラフは双対格子ではスピンの向きが異なるペアにちょうど対応します6。これは、$\omega_n = \Omega_n$、すなわち
Z=e^{2NK}\sum_{n}\Omega_n e^{-2Kn}
であることを示しています。
ここで
Z(K)=(\cosh K)^{2N} 2^N \sum_{n} \Omega_n u^n, \tilde{Z(K)} = e^{2NK}\sum_{n}\Omega_n e^{-2Kn}
とおき、変換
e^{-2\tilde{K}} =\tanh K
を考えます。すると
Z(K)=(\cosh K)^{2N} 2^N e^{-2N\tilde{K}}\tilde{Z}(\tilde{K})=(\sinh(2K))^{N}\tilde{Z}(\tilde{K})
が成り立ち、
\sinh(2\tilde{K}) = \frac{1}{2}(\frac{1}{\tanh K}-\tanh K)=\frac{1}{\sinh(2K)}
を用いると
\frac{Z(K)}{{\sinh(2K)}^{N/2} }=\frac{\tilde{Z}(\tilde{K})}{{\sinh(2\tilde{K})}^{N/2} }
とまとめられます。ここで、$Z(K)$と$\tilde{Z(K)}$は同じ系の分配関数なので一致するはずであり、それを考慮すると
\frac{Z(K)}{{\sinh(2K)}^{N/2}}=\frac{Z(\tilde{K})}{{\sinh(2\tilde{K})}^{N/2} }
です。このような低温展開と高温展開の双対性がKramers-Wannier双対性です。$\beta$を変化させていくと$K,\tilde{K}$も変化していきますが、ある$K$で$Z$に特異性が現れるとすると$\tilde{K}$も特異点になるはずです。相転移点がただ一つだけあるとすれば$K=\tilde{K}$が成り立つ必要があり、このことから相転移点を与える式は
e^{-2K} =\tanh K
です（あるいは$\sinh(2K) = \frac{1}{\sinh(2K)}$から方程式$\sinh(2K)=1$を導いても良いです）。これを解くとよく知られた2次元イジングモデルの臨界温度$\beta_c = \frac{K}{J} = \frac{1}{J}\ln(1+\sqrt{2})$を得ます。
今回考えたいのは非等方的な場合です。その場合の細かい計算は2021年度のPhysics Labのpdf に詳述されているので結果だけを紹介します。$x,y$方向の結合定数を$K,K'$とすると、双対格子への変換の式は
e^{-2\tilde{K'}} = \tanh K, e^{-2\tilde{K}} = \tanh K'
と書けます（つまり2方向の結合定数が混ざるような形になります）。相転移は$K=\tilde{K},K'=\tilde{K'}$が成り立つときに起きます。転移点を決める方程式は、先ほど等方的な場合で$\sinh(2\tilde{K}) = \frac{1}{\sinh(2K)}$を導いたのと同様にして
\sinh(2K) = \sinh(2\tilde{K})=\frac{1}{\sinh(2K')}
つまり
\sinh (2K) \sinh(2K')=1
となります。
それでは、前節の結果を代入して横磁場イジングモデルの量子相転移点を求めましょう。$K,K'$に$\frac{\beta J}{M},-\frac{1}{2}\log \tanh(\frac{\beta \Gamma}{M})$を代入してみます。一見複雑になりそうですが、驚いたことに（？）
\Gamma = J
を得ます。これが求めたかった相転移点に他なりません。
数値計算
最後に数値計算して確かめてみましょう。以下のpythonコードを実行すると、$\Gamma/J$を変化させていった時の磁化の2乗の期待値が計算できます7。なお、学科同期のM君がProcessPoolExecutorを用いた並列実行化、およびscipyの疎行列パッケージを用いた高速化を施してくれたため、そちらを掲載します8。
import matplotlib.pyplot as plt
import numpy as np
import math
from math import sqrt
from scipy.linalg import expm
from scipy.sparse import csr_array, bsr_array, csc_array
from sympy import init_printing, ones, matrix2numpy
from sympy.physics.quantum import represent,TensorProduct
from sympy.physics.quantum.gate import IdentityGate,X,Z
from sympy.matrices.sparse import SparseMatrix
from sympy import Matrix
from concurrent.futures import ProcessPoolExecutor
def toSparse(mat):#疎行列化する
if not isinstance(mat, np.ndarray):
mat = matrix2numpy(mat, dtype=np.float16)
return csr_array(mat)
def expH(t, M, N, J, G):
#初期状態にe^(-tH)を作用させた状態とその状態に収束するまでのエネルギー変化を返す関数
exphz = expmdiagonal(t/M * J * hz(N).toarray())
exphx = toSparse(expmH_x(t, M, N, G))
psi0 = toSparse(np.ones((2**N,1)))
ene = []
for j in range(M):
psi0 = exphz @ exphx @ psi0
norm = sqrt((psi0.T @ psi0)[0,0])
psi0 /= norm
ene.append(energy(psi0, N, J, G))
return psi0, ene
def hz(N):#ハミルトニアンの2体相互作用項
hz = np.zeros((2**N, 2**N))
if N > 2:
for i in range(N - 1):
a = represent(Z(i) * Z(i + 1), nqubits=N)
hz += a
a = represent(Z(N - 1) * Z(0), nqubits=N)
hz += a
return toSparse(hz)
def hx(N):#ハミルトニアンの横磁場項
hx = np.zeros((2**N, 2**N))
for i in range(N):
a = represent(X(i), nqubits=N)
hx += a
return hx
def expmdiagonal(A):#対角行列の指数関数を求める
N = len(A[0])
B = Matrix.zeros(N)
for i in range(N):
B[i, i] = math.exp(A[i, i])
return B
def expmH_x(t, M, N, G):#横磁場部分の指数関数
a = math.cosh(t*G/M)
b = math.sinh(t*G/M)
mat = Matrix([[a,b],[b,a]])
A = mat
for i in range(N-1):
A = TensorProduct(A,mat)
return np.array(A,dtype=np.float16)
def sumspin(N):#スピンの総和を計算する
s = np.zeros((2**N, 2**N))
for i in range(N):
a = represent(Z(i), nqubits=N)
s += a / N
return toSparse(s)
def MM(t, M, N, J, G):#磁化の2乗の期待値を求める
psi = expH(t, M, N, J, G)[0]
result = psi.T @ sumspin(N) @ sumspin(N) @ psi
return result[0, 0]
def energy(psi, N, J, G):#エネルギーの期待値を求める
hamiltonian = -J * hz(N) - G * hx(N)
mat = psi.T @ hamiltonian @ psi
return mat[0, 0]
#エネルギーが収束することの確認
N = 10
ene_result=expH(N,10*N,N,1,1)[1]
plt.plot(ene_result)
plt.xlabel('step number')
plt.ylabel('eneygy')
#磁化の2乗のプロット
if __name__ == "__main__":
init_printing() # ベクトルや行列を綺麗に表示するため
N = 10
Gamma = [0.001, 0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 3.0, 4.0, 5.0]
g_len = len(Gamma)
one_list = np.ones(g_len, dtype=np.int8)
N_list = one_list*N
N10_list = N_list*10
with ProcessPoolExecutor(max_workers=None) as executer:
res = executer.map(MM, N_list, N10_list, N_list, one_list, Gamma)
plt.plot(Gamma, list(res))
plt.title(f"N={N}")
plt.xlabel('$\Gamma /J$') # 理想的にはΓ=Jで相転移
plt.ylabel('Magnetic moment^2')
plt.show()
$N=10$の場合の計算結果を示します。$\Gamma/J=1$の付近で激しく減少しており、$N\rightarrow\infty$では$\Gamma/J=1$で磁化率が発散する振る舞いが想像できますね。これより横磁場イジングモデルは2次転移を起こすことがわかります（これは古典イジングと同じ結果です）。
最後に
イジングモデルの最初に触れたのはB1の秋学期だったと思いますが、そのシンプルな見た目と裏腹に面白い物理的性質がたくさんあるのがとても非自明で惹かれたのを覚えています。物性物理には他にも見かけ以上に奥が深いモデルがたくさんあります。物理って面白いですね〜。
参考文献
-
高橋 和孝、西森 秀稔『相転移・臨界現象とくりこみ群』丸善出版, 2017
-
西森 秀稔、大関真之『量子アニーリングの基礎』共立出版, 2018
-
統計力学第二（求先生）の講義ノート
-
Physics Lab 2021 物性班解説記事 相転移と臨界現象~Ising 模型を題材に~ (https://event.phys.s.u-tokyo.ac.jp/physlab2021/pdf/condmat_Ising.pdf)
-
理物開講「計算機実験2」（藤堂先生）講義資料（実装方法に関するヒントがあったので参考にした）
-
イジングなんて知ってるよ、という人ももう少しお付き合いください。 ↩
-
物工で言うと、例えば芝内・橋本研が研究しているそうです。量子相転移は非常に広くて奥深いトピックであり、これはほんの一例にすぎないことを付記しておきます。 ↩
-
凄まじい計算量です ↩
-
$\sigma^x$のTaylor展開を偶数・奇数の項に分け、$(\sigma^x)^2=1$を使います。 ↩
-
$\sigma_l=\sigma_{l+1}$とおくと$\cosh h = A(h)e^{K(h)}$、$\sigma_l=-\sigma_{l+1}$とおくと$\sinh h = A(h)e^{-K(h)}$です。これを$A(h), K(h)$について解きます。 ↩
-
境界条件によっては一対一に対応しない場合がありますが、自由エネルギー密度などの熱力学量には影響しません。 ↩
-
2乗をとっているのは、たとえ強磁性状態であっても+1と-1に揃った状態が打ち消しあって磁化の期待値が0になってしまうのを回避するためです。 ↩
-
M君以外にもデバッグに協力してくれた学科同期が何人もいます。本当にありがとうございます。 ↩