Statistical Mechanics: Theory and Molecular Simulation Chapter 16 - Critical phenomena

nori

2024年3月6日 21:47

Statistical Mechanics: Theory and Molecular SimulationはMark Tuckermanによって著された分子シミュレーションに関連した熱・統計物理学，及び量子力学の参考書になります。

現時点では（恐らく）和訳がないこともあり，日本での知名度はあまりないような気がしますが，被引用数が約1,600であることを考えると海外では高く認知されている参考書みたいです。

https://scholar.google.com/citations?user=w_7furwAAAAJ

本記事では，Chapter 16（Critical phenomena）の章末問題の和訳とその解答例を紹介します。解答例に間違いが見受けられた場合はお知らせいただけると助かります。

Problem 16.1

外部磁場が存在する場合の一次元イジングモデルのブロックスピン変換を考える。ブロックスピン変換は１つおきにスピンを足し合わせることにする，そのような工程はdecimationと呼ばれる。

a. この変換の変換演算子$${T}$$を書き下し，この変換の長さスケーリング因子（length scaling factor）が$${b=2}$$となることを示せ。

b. この変換のRG（Renormalization Group）方程式を導出し，固定点を見出せ。

c. $${(x,y)}$$平面にRGフローをスケッチせよ。固定点にはどのような性質があり，臨界点の存在について何を暗示するであろうか？

系のハミルトニアンが

$$
\begin{align*}
\mathcal{H}(\{\sigma\})&=-J\sum_{i=1}^N\sigma_i\sigma_{i+1}-H\sum_{i=1}^N\sigma_i-N (C/\beta)
\end{align*}\tag{1.1}
$$

で与えられるとする。ここで，$${C}$$は定数である。
逆温度$${\beta=\frac{1}{k_{\rm B}T}}$$を用いて還元ハミルトニアンを$${\Theta=\beta\mathcal{H}}$$，及び$${K=\beta J}$$，$${h=\beta H}$$とおくと，分配関数$${\Delta(N,h,T)}$$は(1.2)式となる。

$$
\begin{align*}
\Delta(N,h,T)&=\sum_{\sigma_1=\pm 1}\cdots\sum_{\sigma_N=\pm 1}\exp\left[K\sum_{i=1}^N\sigma_i\sigma_{i+1}+h\sum_{i=1}^N\sigma_i+NC\right]
\end{align*}\tag{1.2}
$$

a. 偶数番目のスピンを$${\sigma\rightarrow\sigma'}$$に繰り込み変換する場合，変換演算子$${T}$$は

$$
\begin{align*}
T(\sigma'_{i};\sigma_{2i-1},\sigma_{2i})&=\delta_{\sigma'_{i}\sigma_{2i}}
\end{align*}\tag{1.3}
$$

で与えられる。
隣り合うスピンを繰り込むため，長さスケーリング因子は$${b=2}$$となる。

b. 繰り込み変換後の還元ハミルトニアンを$${\Theta'(\{\sigma'\})}$$とおくと，

$$
\begin{align*}
&{\rm e}^{-\Theta'(\{\sigma'\})}\\
&=\sum_{\sigma_1}\cdots\sum_{\sigma_N}\left(\delta_{\sigma_1'\sigma_2}\delta_{\sigma_2'\sigma_4}\cdots\right){\rm e}^{K\sum_{i=1}^N\sigma_i\sigma_{i+1}+h\sum_{i=1}^N\sigma_i+NC}\\
&=\sum_{\sigma_1}\sum_{\sigma_3}\cdots{\rm e}^{K\sigma_1\sigma_{1}'+\frac{h}{2}(\sigma_1+\sigma_1')+C}{\rm e}^{K\sigma_1'\sigma_{3}+\frac{h}{2}(\sigma_1'+\sigma_3)+C}\cdots\\
&=2^{N/2}\prod_{i=1}^{N/2}{\rm e}^{\frac{h}{2}(\sigma_i'+\sigma_
{i+1'})+2C}\cosh\left(K(\sigma_i'+\sigma_{i+1}')+h\right)
\end{align*}\tag{1.4}
$$

となる。

$$
\begin{align*}
{\rm e}^{\frac{h}{2}(\sigma_i'+\sigma_{i+1}')+2C}\cosh\left(K(\sigma_i'+\sigma_{i+1}')+h\right)&={\rm e}^{K'\sigma_i'\sigma_{i+1}'+\frac{h'}{2}(\sigma_i'+\sigma_{i+1}')+C'}
\end{align*}\tag{1.5}
$$

を満たす$${K',h',C'}$$が求めるRG方程式に対応する。
$${(\sigma_i',\sigma_{i+1}')=(1,1),\ (1,-1),\ (-1,-1)}$$の場合をそれぞれ考えると，

$$
\begin{align*}
{\rm e}^{h+2C}\cosh\left(2K+h\right)&={\rm e}^{K'+h'+C'}\\
{\rm e}^{2C}\cosh h&={\rm e}^{-K'+C'}\\
{\rm e}^{-h+2C}\cosh\left(-2K+h\right)&={\rm e}^{K'-h'+C'}
\end{align*}\tag{1.6}
$$

$$
\begin{align*}
\therefore {\rm e}^{2h'}&={\rm e}^{2h}\frac{\cosh(2K+h)}{\cosh(2K-h)}\\
{\rm e}^{4K'}&=\frac{\cosh(2K+h)\cosh(2K-h)}{\cosh^2h}\\
{\rm e}^{4C'}&={\rm e}^{8C}\cosh(2K+h)\cosh(2K-h)\cosh^2h
\end{align*}\tag{1.7}
$$

$${x:={\rm e}^{-4K},\ y:={\rm e}^{-2h},\ z:={\rm e}^{-8C},\ 0\le \{x,y,z\}\le 1}$$と定義すると，最終的にRG方程式は(1.8)式に帰着する。

$$
\begin{align*}
x'&=\frac{x(1+y)^2}{(x+y)(1+xy)}\\
y'&=\frac{(x+y)y}{1+xy}\\
z'&=\frac{xy^2z^2}{(x+y)(1+xy)(1+y)^2}\\
\end{align*}\tag{1.8}
$$

$${x,y}$$の更新には$${z}$$が寄与しないため，$${x,y}$$のみに着目すればよい。
固定点は$${(x,y)=(0,0),\ (0,1),\ (1,0\le\forall y\le 1)}$$である。

c. RGフロー，及び各固定点の性質に関する説明は参考文献3に解答が与えらてれいるため省略する。

Problem 16.2

一般化された磁気系のモデルはPottsモデルと呼ばれ，各格子に配置されたスピン自由度が$${N}$$個の離散的な状態を取る。１次元で３つの状態を取るPottsモデルは以下のように定義される。番号$${i}$$が割り当てられた各サイトのスピンが$${1,\ 2\ ,3}$$の値を取りうる。ハミルトニアンは

$$
\begin{align*}
\mathcal{H}&=-J\sum_i\delta_{s_i,s_{i+1}}
\end{align*}
$$

で与えられる。
Problem 16.1と同じdecimationを用いてRG方程式を導出し，固定点を見出せ。このモデルには臨界点は存在するか？

$$
\begin{align*}
&{\rm e}^{-\Theta'(\{\sigma'\})}\\
&=\sum_{s_1}\cdots\sum_{s_N}\left(\delta_{s_1's_2}\delta_{s_2's_4}\cdots\right){\rm e}^{K\sum_{i=1}^N\delta_{s_i,s_{i+1}}+NC}\\
&=\sum_{s_1}\sum_{s_3}\cdots{\rm e}^{K\delta_{s_1,s_{1}'}+C}{\rm e}^{K\delta_{s_1',s_3}+C}{\rm e}^{K\delta_{s_3,s_2'}+C}{\rm e}^{K\delta_{s_2',s_5}+C}\cdots\\
\end{align*}\tag{2.1}
$$

例として，$${s_3}$$の和に着目すると，

$$
\begin{align*}
\sum_{s_3=1}^3{\rm e}^{K\delta_{s_1',s_3}+C}{\rm e}^{K\delta_{s_3,s_2'}+C}&=\begin{cases}{\rm e}^{2K+2C}+2{\rm e}^{2C} & {\rm if\ }s_1' = s_2'\\
2{\rm e}^{K+2C}+{\rm e}^{2C} & {\rm if\ }s_1' \neq s_2'\end{cases}
\end{align*}\tag{2.2}
$$

つまり，(2.2)式と整合するように$${{\rm e}^{K'\delta_{s_1',s_2'}+C'}}$$の$${K',\ C'}$$を決めることを考える。

$$
\begin{align*}
\therefore R(K)&=\ln\left(\frac{{\rm e}^{2K}+2}{2{\rm e}^{K}+1}\right)
\end{align*}\tag{2.3}
$$

固定点は$${K=0, \infty}$$になり，それぞれ高温極限，低温極限に対応する。

Problem 16.3

ファンデアワールス理論におけるマクスウェルの等面積則では$${T< T_c}$$において$${\partial P/\partial V >0}$$から生じる理論上の非物理的な振る舞いを修正することを試みる。手順の概略をFig. 16.20に示す。

1次相転移による体積の不連続な変化を表す共役線を導入する。共役線の箇所は，共役線と$${PV}$$曲線によって囲まれる２つの領域の面積が等しくなる（$${a_1=a_2}$$）ように決められる。(4.7.35)式のファンデアワールス方程式とマクスウェルの等面積則を用いることにより，平均場の臨界指数$${\beta}$$の値が$${1/2}$$となることを示せ。

$$
\begin{align*}
\frac{P}{k_{\rm B}T}&=\frac{\rho}{1-\rho b}-\frac{a\rho^2}{k_{\rm B}T}\\
\rho&=\frac{N}{V}
\end{align*}\tag{4.7.35}
$$

解答は参考文献4にそのものが掲載されているため，ここでは略。。。

Problem 16.4

外部磁場がない1/2スピンの１次元イジング模型において，分配関数は以下の形式で表すことができる。

$$
\begin{align*}
\Delta&={\rm Tr}\left(P^N\right)
\end{align*}
$$

pair cell変換と呼ばれるRG変換を考えた場合，$${\Delta}$$は以下のように再定式化される。

$$
\begin{align*}
\Delta&={\rm Tr}\left[\left(P^2\right)^{N/2}\right]
\end{align*}
$$

転送行列は$${P'=P^2}$$に再定義される。
この変換に対するRG方程式を見出し，期待される安定固定点が得られることを示せ。

$${\Theta_0=K\sum_{i=1}^{N}\sigma_i\sigma_{i+1}+NC}$$とした場合，

$$
\begin{align*}
P&=\begin{bmatrix}{\rm e}^{K+C}&{\rm e}^{-K+C}\\{\rm e}^{-K+C}&{\rm e}^{K+C}\end{bmatrix}\\
P^2&=2{\rm e}^{2C}\begin{bmatrix}\cosh 2K&1\\1&\cosh 2K\end{bmatrix}\\
&\Longrightarrow\begin{bmatrix}{\rm e}^{K'+C'}&{\rm e}^{-K'+C'}\\{\rm e}^{-K'+C'}&{\rm e}^{K'+C'}\end{bmatrix}
\end{align*}\tag{4.1}
$$

$$
\begin{align*}
\therefore {\rm e}^{2K'}&=\cosh 2K
\end{align*}\tag{4.2}
$$

$${{\rm e}^{2K}=\cosh 2K}$$のとき，$${K=0, \infty}$$となり，期待される固定点が得られる。

Problem 16.5

以下の仮定に基づく長いポリマー鎖の単純モデルを考える：

鎖の配座エネルギー$${E}$$は単に主鎖の二面角のみから決定される。
各主鎖は「trans」配座を表す$${t}$$，及び「gauche」配座を表す$${g^+,\ g^-}$$の３つの状態のいずれかのみを取る。
各配座は固有のエネルギーを持ち，また隣り合う二面角のみの配座に影響を及ぼす。ポリマーが$${N}$$個の原子からなる場合，$${N-3}$$個の二面角$${\phi_1,\cdots,\phi_{N-3}}$$が存在することになる。全エネルギー$${E(\phi_1,\cdots,\phi_{N-3})}$$は以下のようになる。

$$
\begin{align*}
E(\phi_1,\cdots,\phi_{N-3})&=\sum_{i=1}^{N-3}\varepsilon_1(\phi_i)+\sum_{i=2}^{N-3}\varepsilon_2(\phi_{i-1}, \phi_i)
\end{align*}
$$

2つのエネルギー$${\varepsilon_1,\ \varepsilon_2}$$は以下の値を取るものとする。

$$
\begin{align*}
\varepsilon_1(t)&=0\\
\varepsilon_1(g^+)&=\varepsilon_1(g^-)=\varepsilon\\
\varepsilon_2(g^+,g^-)&=\varepsilon_2(g^-,g^+)=\infty\\
\varepsilon_2(\phi_{i-1}, \phi_i)&=0\ \ \ \ \ \ {\rm for\ all\ other\ combinations}
\end{align*}
$$

a. この系のカノニカル分配関数を計算せよ。その際，$${\sigma:=\exp(-\beta\varepsilon)}$$の表記を利用せよ。

b. $${N\rightarrow\infty}$$の極限において，分配関数が

$$
\begin{align*}
\lim_{N\rightarrow\infty}\frac{1}{N}\ln Q&=\ln\chi\\
\chi&=\frac{1}{2}\left[(1+\sigma)+\sqrt{1+6\sigma+\sigma^2}\right]
\end{align*}
$$

と振る舞うことを示せ。

c. $${N\rightarrow\infty}$$の場合に全ての二面角がtrans配座を取る確率を求めよ。

d. $${N\rightarrow\infty}$$の場合に二面角がtrans配座とgauche配座を交互にを取る確率を求めよ。

a. 行列要素$${\lang \phi|P|\phi'\rang:={\rm e}^{-\beta(\frac{\varepsilon_1(\phi)+\varepsilon_1(\phi')}{2}+\varepsilon_2(\phi,\phi'))}}$$を有する行列$${\mathbf{P}}$$を導入する。

$$
\begin{align*}
\mathbf{P}&=\begin{bmatrix}1&\sigma^{1/2}&\sigma^{1/2}\\
\sigma^{1/2}&\sigma&0\\
\sigma^{1/2}&0&\sigma\end{bmatrix}
\end{align*}\tag{5.1}
$$

$${\mathbf{P}}$$の固有値を$${\lambda}$$とおくと，

$$
\begin{align*}
|\mathbf{P}-\lambda\mathbf{I}|&=\begin{vmatrix}1-\lambda&\sigma^{1/2}&\sigma^{1/2}\\
\sigma^{1/2}&\sigma-\lambda&0\\
\sigma^{1/2}&0&\sigma-\lambda\end{vmatrix}\\
&=(\sigma-\lambda)\{\lambda^2-(1+\sigma)\lambda-\sigma\}\\
&=0\\
\therefore \lambda&=\sigma,\ \sigma_{\pm}\\
\sigma_{\pm}&=\frac{1}{2}\left[(1+\sigma)\pm\sqrt{1+6\sigma+\sigma^2}\right]
\end{align*}\tag{5.2}
$$

また，各固有ベクトルを$${\mathbf{e}_{\sigma},\ \mathbf{e}_{\sigma_+},\ \mathbf{e}_{\sigma_-}}$$，及びそれらを列成分に持つ行列を$${\mathbf{V}=\begin{bmatrix}\mathbf{e}_{\sigma}&\mathbf{e}_{\sigma_+}&\mathbf{e}_{\sigma_-}\end{bmatrix}}$$とする。
以上を用いて分配関数$${Q(N,T)}$$を計算すると，

$$
\begin{align*}
&Q(N,T)\\
&=\sum_{\phi_1}\cdots\sum_{\phi_{N-3}}\exp\left(-\beta\left[\sum_{i=1}^{N-3}\varepsilon_1(\phi_i)+\sum_{i=2}^{N-3}\varepsilon_2(\phi_{i-1}, \phi_i)\right]\right)\\
&=\sum_{\phi_1}\cdots\sum_{\phi_{N-3}}{\rm e}^{-\beta\varepsilon_1(\phi_1)}\lang\phi_1|P|\phi_2\rang\cdots\lang\phi_{N-4}|P|\phi_{N-3}\rang{\rm e}^{-\beta\varepsilon_1(\phi_{N-3})}\\
&=\sum_{\phi_1}\sum_{\phi_{N-3}}{\rm e}^{-\beta\varepsilon(\phi_1)}\lang\phi_1|P^{N-4}|\phi_{N-3}\rang{\rm e}^{-\beta\varepsilon_1(\phi_{N-3})}
\end{align*}\tag{5.3}
$$

$${\mathbf{a}:=\begin{bmatrix}1&{\rm e}^{-\beta\varepsilon/2}&{\rm e}^{-\beta\varepsilon/2}\end{bmatrix}^{\rm T}}$$とおくと，(5.3)式は

$$
\begin{align*}
{\rm r.h.s}&=\mathbf{a}^{\rm T}\left(\sum_{\lambda=\sigma,\sigma_+,\sigma_-}\lambda^{N-4}\mathbf{e}_\lambda\mathbf{e}_\lambda^{\rm T}\right)\mathbf{a}\\
&=\sum_{\lambda=\sigma,\sigma_+,\sigma_-}\lambda^{N-4}\left\|\mathbf{e}_\lambda^{\rm T}\mathbf{a}\right\|^2
\end{align*}\tag{5.4}
$$

に帰着する。

b. $${\sigma={\rm e}^{-\beta}\varepsilon}$$であることを思い出すと，$${N\rightarrow\infty}$$において$${\sigma^N, \sigma_-^N\rightarrow 0}$$となる。
そのため，

$$
\begin{align*}
Q(N,T)&\simeq \sigma_+^{N-4}\left\|\mathbf{e}_{\sigma_+}^{\rm T}\mathbf{a}\right\|^2
\end{align*}\tag{5.5}
$$

$$
\begin{align*}
\therefore\lim_{N\rightarrow\infty}\frac{1}{N}\ln Q&=\lim_{N\rightarrow\infty}\frac{N-4}{N}\ln \sigma_++\frac{\left\|\mathbf{e}_{\sigma_+}^{\rm T}\mathbf{a}\right\|^2}{N}\\
&=\ln\sigma_+\\
&=\ln\left\{\frac{1}{2}\left[(1+\sigma)+\sqrt{1+6\sigma+\sigma^2}\right]\right\}
\end{align*}\tag{5.6}
$$

$$
\begin{align*}
P&\simeq\frac{1}{ \sigma_+^{N}\left\|\mathbf{e}_{\sigma_+}^{\rm T}\mathbf{a}\right\|^2}
\end{align*}\tag{5.7}
$$

$$
\begin{align*}
P&\simeq\frac{2\left(2{\rm e}^{-\beta\epsilon}\right)^{\frac{N}{2}}}{ \sigma_+^{N}\left\|\mathbf{e}_{\sigma_+}^{\rm T}\mathbf{a}\right\|^2}
\end{align*}\tag{5.8}
$$

Problem 16.6

各格子上のスピンが$${z}$$軸成分のみを有し，また隣り合うスピン間にのみ相互作用が存在するイジング模型のハミルトニアンは(16.3.3)式で与えられる。

$$
\begin{align*}
\mathcal{H}&=-\frac{1}{2}\sum_{<i,j>}J_{ij}\sigma_i\sigma_j-h\sum_{i}\sigma_i
\end{align*}\tag{16.3.3}
$$

各スピン変数$${\sigma_i}$$は$${-1,\ 0,\ 1}$$の３つの値を取ることができる。
平均場理論を用いて磁化の超越方程式を見出し，このモデルの臨界温度を決定せよ。臨界指数はどうなるであろうか？

$${J_{ij}=J,\ m:=\frac{1}{N}\lang\sum_i\sigma_i\rang}$$とおくと，平均場近似においてハミルトニアンは

$$
\begin{align*}
\mathcal{H}&\simeq\frac{NJm^2}{2}-(Jm+h)\sum_{i}\sigma_i
\end{align*}\tag{6.1}
$$

と近似される。
(6.1)式を用いて系の分配関数$${\Delta(N,h,T)}$$，１スピン当たりのギブス自由エネルギー$${g(h,T)}$$を計算する。

$$
\begin{align*}
\Delta(N,h,T)&={\rm e}^{-\beta NJm^2/2}\left(\prod_{i=1}^N\sum_{\sigma_i=0,\pm 1}{\rm e}^{\beta(Jm+h)\sigma_i}\right)\\
&={\rm e}^{-\beta NJm^2/2}\left(2\cosh(\beta(Jm+h))+1\right)^N\\
\end{align*}\tag{6.2}
$$

$$
\begin{align*}
g(h,T)&=-\frac{1}{N\beta}\ln\Delta(N,h,T)\\
&=\frac{Jm^2}{2}-\frac{1}{\beta}\ln\left[2\cosh(\beta(Jm+h))+1\right]\\
\end{align*}\tag{6.3}
$$

(6.3)式から$${h=0}$$での磁化率を求めると，

$$
\begin{align*}
m&=-\left.\frac{\partial g}{\partial h}\right|_{h=0}\\
&=\frac{2\sinh(\beta Jm)}{2\cosh(\beta Jm)+1}
\end{align*}\tag{6.4}
$$

臨界温度の場合，(6.4)式の右辺の$${m=0}$$の切片が$${1}$$になる。

$$
\begin{align*}
\therefore T_c&=\frac{2J}{3k_{\rm B}}
\end{align*}\tag{6.5}
$$

$${g(0,T)}$$を$${m}$$に関してテイラー展開すると，$${m,\ m^3}$$の項は$${0}$$となるため，スピン自由度が$${2}$$の場合と同様にして臨界指数$${\beta=1/2}$$となる。

Problem 16.7

1978年にH. J. MarisとL. J. Kadanoffが外部磁場が存在しない２次元イジング模型に対する繰り込み群の方法を導入した（Maris and Kadanoff, 1978; Chandler, 1987）。この問題ではこの方法によるRG方程式の導出を一歩ずつ検証していくことにする。

a. Fig. 16.13(a)で示されているようなラベル付けされたスピンが２次元周期格子上に存在する状況を考える。

Fig. 16.13: (a) ラベル付けされた２次元格子状のスピン。(b)decimation後のスピン格子。（参考文献１より転載）

最初の手順では，足し合わされる各スピンが単一のボルツマン因子のみに現れるよう分配関数の総和を分配することによって，格子状のスピン数が半分になるように足し合わせる。
結果として得られるスピン格子がFig. 16.13(b)に対応する場合，分配関数が以下の数式となることを示せ。

$$
\begin{align*}
Q&=\sum_{\rm remaining\ spins}\cdots\left[{\rm e}^{K(\sigma_1+\sigma_2+\sigma_3+\sigma_4)}+{\rm e}^{-K(\sigma_1+\sigma_2+\sigma_3+\sigma_4)}\right]\\
&\ \ \ \ \ \ \ \ \ \ \times\left[{\rm e}^{K(\sigma_2+\sigma_3+\sigma_7+\sigma_8)}+{\rm e}^{-K(\sigma_2+\sigma_3+\sigma_7+\sigma_8)}\right]\cdots
\end{align*}
$$

ここで，$${K=\beta J}$$である。

b. [ ]内の項を

$$
\begin{align*}
{\rm e}^{K(\sigma_1+\sigma_2+\sigma_3+\sigma_4)}+{\rm e}^{-K(\sigma_1+\sigma_2+\sigma_3+\sigma_4)}&=g(K){\rm e}^{K'(\sigma_1\sigma_2+\sigma_2\sigma_3+\sigma_3\sigma_4+\sigma_4\sigma_1)}
\end{align*}
$$

のように書き下せるかを考える。ここで，$${K'}$$は新しい結合定数であり，また$${g(K)}$$は以下のスピンの値を代入した際に成立することを要求することで決定される。

$$
\begin{align*}
\sigma_1&=&\sigma_2&=&\sigma_3&=&\sigma_4&=&\pm 1\\
\sigma_1&=&\sigma_2&=&\sigma_3&=&-\sigma_4&=&\pm 1\\
\sigma_1&=&\sigma_2&=&-\sigma_3&=&-\sigma_4&=&\pm 1\\
\sigma_1&=&-\sigma_2&=&\sigma_3&=&-\sigma_4&=&\pm 1\\
\end{align*}
$$

上記の手順では$${g(K),\ K'}$$を決定できないことを示せ。

c. 代わりに複数の結合定数$${K_1,\ K_2,\ K_3}$$を導入し，以下の方程式を満足する$${g(K),\ K_1,\ K_2,\ K_3}$$を探すことを考える。

$$
\begin{align*}
{\rm e}^{K(\sigma_1+\sigma_2+\sigma_3+\sigma_4)}+{\rm e}^{-K(\sigma_1+\sigma_2+\sigma_3+\sigma_4)}&=g(K){\rm e}^{\frac{K_1}{2}(\sigma_1\sigma_2+\sigma_2\sigma_3+\sigma_3\sigma_4+\sigma_4\sigma_1)+K_2(\sigma_1\sigma_3+\sigma_2\sigma_4)+K_3\sigma_1\sigma_2\sigma_3\sigma_4}
\end{align*}
$$

b.で言及したスピンの値を挿入することにより，$${g(K),\ K_1,\ K_2,\ K_3}$$を$${K}$$の関数として表せ。得られる分配関数の意味を考察せよ。

d. c.で得られた結果は厳密にはRGの手順にはならないが，$${K_2,\ K_3}$$が無視できるのであれば以下のRG方程式に帰着することを示せ。

$$
\begin{align*}
K_1&=\frac{1}{4}\ln\cosh(4K)
\end{align*}
$$

これは臨界点に繋がるであろうか？

e. d.の結果を改善するために，無視するのは$${K_3}$$のみにして$${K_2}$$は近似的に取り扱うことが提案される。分配関数の数式内で$${K_1}$$の係数を持つスピンの和を$${\Sigma_1}$$，$${K_2}$$の係数を持つスピンの和を$${\Sigma_2}$$としよう。以下の近似を考える。

$$
\begin{align*}
K_1\Sigma_1+K_2\Sigma_2&\simeq K'(K_1, K_2)\Sigma_{i,j}'\sigma_i\sigma_j
\end{align*}
$$

ここで，$${\Sigma_{i,j}'}$$は近隣点のみを対象に合計することを表す。この近似で得られる分配関数を求めよ。

f. スピン当たりの自由エネルギーを

$$
\begin{align*}
a(K)&=\frac{1}{N}\ln Q(K)
\end{align*}
$$

と定義する。
$${a(K)}$$が以下のRG方程式を満たすことを示せ。

$$
\begin{align*}
a(K)&=\frac{1}{2}\ln g(K)+\frac{1}{2}f(K')
\end{align*}
$$

g. $${K'\simeq K_1+K_2}$$を示せ。

h. $${K'}$$のRG方程式を導出し，臨界点の存在を予測するモデルとなっていることを示せ。また，$${K_c}$$の値を求めよ。

i. $${K=K_c}$$まわりで自由エネルギーをテイラー展開することにより，臨界指数$${\alpha}$$を計算し，オンサーガによる厳密解と比較せよ。

j. 最後に，得られた臨界温度とオンサーガの結果を比較せよ：

$$
\begin{align*}
\frac{J}{k_{\rm B}T_c}&=0.44069
\end{align*}
$$

k. １次元スピン格子上のスピンを一つ飛ばしで足し合わせることにより，磁場のない１次元イジング模型に対してMaris-Kadanofスキームの類似方法を考案せよ。どのようなRG方程式が得られるだろうか？得られた方程式によって期待される固定点が得られることを示せ。

a. 例として，$${\sigma_5}$$に関する和を実行すると，

$$
\begin{align*}
\sum_{\sigma_5=\pm 1}{\rm e}^{K(\sigma_1\sigma_5+\sigma_2\sigma_5+\sigma_5\sigma_3+\sigma_5\sigma_4)}&={\rm e}^{K(\sigma_1+\sigma_2+\sigma_3+\sigma_4)}+{\rm e}^{-K(\sigma_1+\sigma_2+\sigma_3+\sigma_4)}
\end{align*}\tag{7.1}
$$

が得られる。
$${\sigma_6}$$以降に関しても同様に和を実行することにより，与えられた表式に帰着する。

b. $${K', g(K)}$$は以下の方程式を満たす必要がある。

$$
\begin{align*}
2\cosh(4K)&=g(K){\rm e}^{4K'}\\
2\cosh(2K)&=g(K)\\
2&=g(K)\\
2&=g(K){\rm e}^{-4K'}\\
\end{align*}\tag{7.2}
$$

2つの未知数に対して４つの方程式が存在するため，一意に決定することができない。

$$
\begin{align*}
2\cosh(4K)&=g(K){\rm e}^{2K_1+2K_2+K_3}\\
2\cosh(2K)&=g(K){\rm e}^{-K_3}\\
2&=g(K){\rm e}^{-2K_2+K_3}\\
2&=g(K){\rm e}^{-2K_1+2K_2+K_3}\\
\end{align*}\tag{7.3}
$$

$$
\begin{align*}
\therefore K_1&=\frac{1}{4}\ln\cosh(4K)\\
K_2&=\frac{1}{8}\ln\cosh(4K)\\
K_3&=\frac{1}{8}\ln\cosh(4K)-\frac{1}{2}\ln\cosh(2K)\\
g(K)&=2\cosh^{1/2}(2K)\cosh^{1/8}(4K)
\end{align*}\tag{7.4}
$$

d. (7.4)式より，$${K_2,\ K_3}$$を無視できる場合はRG方程式が

$$
\begin{align*}
K_1&=\frac{1}{4}\ln\cosh(4K)\\
\end{align*}\tag{7.5}
$$

となるのは明らかである。
(7.5)式の固定点は$${K=0}$$のみであり，これは高温極限に対応する。
そのため，この近似では自発磁化の存在を予言できない。

$$
\begin{align*}
Q&=\sum_{\rm remaining\ spins}\cdots\left[{\rm e}^{K(\sigma_1+\sigma_2+\sigma_3+\sigma_4)}+{\rm e}^{-K(\sigma_1+\sigma_2+\sigma_3+\sigma_4)}\right]\cdots\\
&=g^{N/2}(K)\sum_{\rm remaining\ spins}\cdots{\rm e}^{\frac{K_1}{2}(\sigma_1\sigma_2+\sigma_2\sigma_3+\sigma_3\sigma_4+\sigma_4\sigma_1)+K_2(\sigma_1\sigma_3+\sigma_2\sigma_4)+K_3\sigma_1\sigma_2\sigma_3\sigma_4}\cdots\\
& \\
&\simeq g^{N/2}(K)\sum_{\rm remaining\ spins}\cdots{\rm e}^{\frac{K_1}{2}(\sigma_1\sigma_2+\sigma_2\sigma_3+\sigma_3\sigma_4+\sigma_4\sigma_1)+K_2(\sigma_1\sigma_3+\sigma_2\sigma_4)}\cdots\\
&\simeq g^{N/2}(K)\sum_{\rm remaining\ spins}\cdots{\rm e}^{K'(K_1, K_2)\Sigma_{i,j}'\sigma_i\sigma_j}\cdots\\
\end{align*}\tag{7.6}
$$

$$
\begin{align*}
\therefore Q(K,N)&\simeq g^{N/2}(K)Q(K',N/2)
\end{align*}\tag{7.7}
$$

$$
\begin{align*}
a(K)&=\frac{1}{N}\ln Q(K,N)\\
&=\frac{1}{2}\ln g(K)+\frac{1}{2}\frac{1}{N/2}\ln Q(K',N/2)\\
&=\frac{1}{2}\ln g(K)+\frac{1}{2}\ln f(K')\\
f(K')&:=\frac{1}{N/2}Q(K',N/2)=a(K')
\end{align*}\tag{7.8}
$$

g. 興味ある状態は隣り合うスピン同士の向きがほぼ揃っている状態であり，その場合$${\Sigma_1,\ \Sigma_2,\ \Sigma_{i,j}'\sigma_i\sigma_j\simeq N}$$が成立する。

$$
\begin{align*}
K_1\Sigma_1+K_2\Sigma_2&\simeq K'(K_1, K_2)\Sigma_{i,j}'\sigma_i\sigma_j\\
K_1N+K_2N&\simeq K'(K_1, K_2)N\\
\therefore K'(K_1, K_2)&\simeq K_1+K_2
\end{align*}\tag{7.9}
$$

h. (7.4)，(7.9)式より，

$$
\begin{align*}
K'&=\frac{3}{8}\ln\cosh(4K)
\end{align*}\tag{7.10}
$$

固定点を$${K_c}$$とおくと，

$$
\begin{align*}
K_c&=\frac{3}{8}\ln\cosh(4K_c)
\end{align*}\tag{7.11}
$$

(7.11)式の形状を図7.1に示す。$${0}$$以外の固定点$${K_c=0.50698}$$が存在し，臨界点を表すと期待される。

i. 参考文献5によると，

$$
\begin{align*}
\alpha&=2-\frac{\ln2}{\ln\left.\frac{\partial K'}{\partial K}\right|_{K=K_c}}\\
&=2-\frac{\ln2}{\ln\left[\frac{3}{2}\tanh(K_c)\right]}\\
&=2-\frac{\ln2}{\ln\left[\frac{3}{2}\tanh(0.50698)\right]}\\
&\simeq 0.131
\end{align*}\tag{7.12}
$$

になるらしい。。。

j. $${K_c=0.50698}$$より，得られた転移温度はオンサーガの厳密解より少し低い値を取る。近似を導入したことにより，本来持つ強相関の性質が少し弱まったことに由来すると考えられる。

$$
\begin{align*}
&Q(N,T)\\
&=\sum_{\sigma_1}\sum_{\sigma_2}\cdots\sum_{\sigma_N}{\rm e}^{K\sum_{i=1}^N\sigma_i\sigma_{i+1}}\\
&=\sum_{\sigma_1}\sum_{\sigma_3}\cdots\left[{\rm e}^{K(\sigma_1+\sigma_3)}+{\rm e}^{-K(\sigma_1+\sigma_3)}\right]\left[{\rm e}^{K(\sigma_3+\sigma_5)}+{\rm e}^{-K(\sigma_3+\sigma_5)}\right]\cdots\\
&=2^{N/2}\sum_{\sigma_1}\sum_{\sigma_3}\cdots\cosh(K(\sigma_1+\sigma_3))\cosh(K(\sigma_3+\sigma_5))\cdots
\end{align*}\tag{7.13}
$$

$${\cosh(K(\sigma_i+\sigma_{i+1}))={\rm e^{K_1\sigma_i\sigma_{i+1}+K_2}}}$$を満たす繰り込み変換を考えると，

$$
\begin{align*}
\cosh(2K)&={\rm e}^{K_1+K_2}\\
1&={\rm e}^{-K_1+K_2}\\
\therefore K_1&=K_2=\frac{1}{2}\ln\cosh(2K)
\end{align*}\tag{7.14}
$$

これは，Problem 16.4で得られたRG方程式と等しく，固定点は$${K_c=0,\ \infty}$$の２つである。

Problem 16.8

7.6節のWang-Landau法（Wang and Landau, 2001）は容易にスピン格子に適用できる。離散的なスピンの値を採用することにより，スピン格子の全エネルギー$${E}$$も同様に離散的となる。そのため，カノニカル分配関数を

$$
\begin{align*}
Q(\beta)&=\sum_{E}\Omega(E){\rm e}^{-\beta E}
\end{align*}
$$

のように表記することができる，ここで$${\Omega(E)}$$は状態密度である。
ランダムに選ばれたスピンの向きを反転させ，(7.6.3)式に従ってその反転が受容されるかを決定する試行を行う。

$$
\begin{align*}
A(E_2|E_1)&=\min\left[1,\frac{\Omega(E_1)}{\Omega(E_2)}\right]
\end{align*}
$$

試行後のエネルギー$${E_i\ (i=1\ {\rm or}\ 2)}$$の状態密度は$${\ln\Omega(E_i)\rightarrow\ln\Omega(E_i)+\ln f}$$に変換される。
外場がない場合の$${50\times 50}$$スピン格子を対象に，複数の温度に対する分配関数と１スピンの自由エネルギーを計算するためのWang-Landauモンテカロルコードを書き下せ。各スピンがランダムに選ばれた値を初期状態として用意し，格子に対して周期境界条件を課すこと。

pythonによるコード例を以下に示す。問題文に与えられた条件をそのまま用いた場合は収束性に難があったため，スピン格子のサイズを$${10\times 10}$$に縮小した。

import itertools
import numpy as np
import sys

np.random.seed(2023)

# Input parameters
lnf          = 1         # initial scaling factor
lattice_size = 10        # size of 2D spin lattice 
n_trial      = 10000000  # The max value of spin flip trial
n_grids      = 51        # grid points of histgram and DOS

def flipSpinOrientation(spin_lattice):
    spin_lattice_copy = spin_lattice.copy()
    i, j = np.random.choice(lattice_size, 2)
    spin_lattice_copy[i, j] *= -1

    return spin_lattice_copy
    
def calculateTotalEnergy(spin_lattice):
    E = 0
    for i in range(lattice_size):
        ii = 0 if i == lattice_size - 1 else i + 1
        E -= np.sum(spin_lattice[i, :] * spin_lattice[ii, :])
        E -= np.sum(spin_lattice[:, i] * spin_lattice[:, ii])

    return E

def judgeAcceptance(lnOmega1, lnOmega2):
    dlnOmega = lnOmega1 - lnOmega2
    if dlnOmega > 0:
        return True
    elif np.exp(dlnOmega) > np.random.rand():
        return True
    else:
        return False

# generate initial values of spin lattice
spin_lattice = np.random.rand(lattice_size ** 2).reshape((lattice_size,) * 2)
spin_lattice = np.where(spin_lattice >= 0.5, 1, -1)

# prepare initial density of states
lnOmega = np.zeros(n_grids)
grids = np.linspace(-2, 2, n_grids + 1)
# prepare histograme
histo = np.zeros(n_grids)

E1 = calculateTotalEnergy(spin_lattice) / (lattice_size ** 2)
grid_id1 = int((0.249 * E1 + 0.5) * n_grids)
lnOmega[grid_id1] += lnf
histo[grid_id1] += 1
for n in range(n_trial):
    spin_lattice_candidate = flipSpinOrientation(spin_lattice)
    E2 = calculateTotalEnergy(spin_lattice_candidate) / (lattice_size ** 2)
    grid_id2 = int((0.249 * E2 + 0.5) * n_grids)

    if judgeAcceptance(lnOmega[grid_id1], lnOmega[grid_id2]):
        spin_lattice = spin_lattice_candidate
        grid_id1 = grid_id2
        E1 = E2

    lnOmega[grid_id1] += lnf
    histo[grid_id1] += 1

    if np.min(histo) >= 0.8 * np.mean(histo):
        if np.max(lnOmega) > 600:
            lnOmega -= (np.max(lnOmega) - 600)
        else:
            lnOmega -= np.max(lnOmega)
        lnOmega += np.log(2) - lnOmega[0]

        if lnf < 1e-8:
            break
        else:
            print(f"lnf:{lnf} -> {lnf * 0.5} at {n} steps")
            lnf *= 0.5
            histo = np.zeros(n_grids)

# Output logarithm of density of states of reduced spin lattice
for i in range(n_grids):
    print(0.5 * (grids[i] + grids[i + 1]), lnOmega[i])

# Calculate partition function and free energy per spin 
beta = np.arange(0.1, 10.01, 0.1)
FE = []
grids2 = np.array([0.5 * (grids[i] + grids[i + 1]) * (lattice_size ** 2) for i in range(n_grids)])
for i in range(len(beta)):
    exponent = lnOmega - beta[i] * grids2
    offset = np.max(exponent)
    Q = np.sum(np.exp(exponent - offset))
    FE.append(-1 * (offset + np.log(Q)) / beta[i] / (lattice_size ** 2))
print("\n")
for a, b in zip(beta, FE):
    print(a, b)

Problem 16.9

外部磁場$${h}$$が存在する場合の２次元イジング模型の分布関数をサンプリングするための単純な$${\rm M(RT)^2}$$モンテカルロプログラムを書き下せ。$${h}$$の大きさを変えた時にスピン―スピン相関関数$${\lang\sigma_i\sigma_j\rang-\lang\sigma_i\rang\lang\sigma_j\rang}$$の振る舞いがどう変わるかを$${T > T_c}$$と$${T < T_c}$$の両方に対して観察せよ。

スピン結合定数$${J}$$とし，$${K=\beta J}$$，$${h'=h/J\rightarrow h}$$とする。
pythonによるコード例を以下に示す。スピン―スピン相関関数は$${\lang\sigma_i\sigma_j\rang-\lang\sigma_i\rang\lang\sigma_j\rang}$$ではなく$${\lang\sigma_i\sigma_j\rang}$$で評価している。

import itertools
import numpy as np
import sys

np.random.seed(2023)

# Input parameters
K            = 0.3       # J * beta
h            = 0.5       # h / J
lattice_size = 20        # size of 2D spin lattice 
n_steps_eq   = 200       # Iteration steps for equilibration
n_steps_samp = 200       # Sampling steps
    
def calculateTotalEnergy(spin_lattice, K, h):
    E = 0
    for i in range(lattice_size):
        ii = 0 if i == lattice_size - 1 else i + 1
        E -= K * np.sum(spin_lattice[i, :] * spin_lattice[ii, :])
        E -= K * np.sum(spin_lattice[:, i] * spin_lattice[:, ii])
    E -= K * h * np.sum(spin_lattice)
    
    return E

def judgeAcceptance(E1, E2):
    exponent = E1 - E2

    if exponent > 0:
        return True
    elif np.exp(exponent) > np.random.rand():
        return True
    else:
        return False
        
def countSpinSpinCorrelation(spin_lattice, ss_corr):
    for i, j in itertools.product(range(lattice_size), repeat = 2):
        for k, l in itertools.product(range(i, lattice_size), range(j, lattice_size)):
            length2 = (i - k) ** 2 + (j - l) ** 2
            if not length2 in ss_corr.keys():
                ss_corr[length2] = [spin_lattice[i,j] * spin_lattice[k,l]]
            else:
                ss_corr[length2].append(spin_lattice[i,j] * spin_lattice[k,l])
    
    return ss_corr

# generate initial values of spin lattice
spin_lattice = np.random.rand(lattice_size ** 2).reshape((lattice_size,) * 2)
spin_lattice = np.where(spin_lattice >= 0.5, 1, -1)
#spin_lattice = np.ones(lattice_size ** 2).reshape((lattice_size,) * 2)

E1 = calculateTotalEnergy(spin_lattice, K, h) 
for _ in range(n_steps_eq):
    for i, j in itertools.product(range(lattice_size), repeat = 2):
        spin_lattice_candidate = spin_lattice.copy()
        spin_lattice_candidate[i,j] *= -1
        E2 = calculateTotalEnergy(spin_lattice_candidate, K, h)
        if judgeAcceptance(E1, E2):
            spin_lattice = spin_lattice_candidate
            E1 = E2

total_spin = []
ss_corr = {}
spin_lattice_sampling = []
for _ in range(n_steps_samp):
    for i, j in itertools.product(range(lattice_size), repeat = 2):
        spin_lattice_candidate = spin_lattice.copy()
        spin_lattice_candidate[i,j] *= -1
        E2 = calculateTotalEnergy(spin_lattice_candidate, K, h)
        if judgeAcceptance(E1, E2):
            spin_lattice = spin_lattice_candidate
            E1 = E2
    total_spin.append(np.sum(spin_lattice))
    ss_corr = countSpinSpinCorrelation(spin_lattice, ss_corr)

ave_spin = np.mean(total_spin) / (lattice_size ** 2)
ss_corr_keys = list(ss_corr.keys())
ss_corr_keys.sort()
for ss_corr_key in ss_corr_keys:
    mean_ss_corr = np.mean(ss_corr[ss_corr_key])# - ave_spin ** 2
    print(np.sqrt(ss_corr_key), mean_ss_corr)

$${K=0.3\ (T > T_c)}$$の場合のスピン―スピン相関関数の$${h}$$依存性を図9.1に示す。

$${K=1.0\ (T < T_c)}$$の場合のスピン―スピン相関関数の$${h}$$依存性を図9.2に示す。

参考文献

Mark E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation (Oxford Graduate Texts)
Michael Cross, Caltech: Physics 127b (Statistical Physics, Caltech: Physics 127b)
Masatsugu Sei Suzuki, Renormalization group for one dimensional Ising model
物理とはずがたり
David Chandler, Introduction to Modern Statistical Mechanices (Oxford University Press)