SK Model in Statistical Physics with Replica Method

Posted on 2023-08-02 In Statistical Physics Views:

Solving the analytical expression of the SK Model using the replica method.(in Chinese)

Sherrington–Kirkpatrick Model

Sherrington–Kirkpatrick模型(以下简称SK模型)是相互作用范围为无限的Edwards-Anderson模型，是描述自旋玻璃(Spin-Glass)的平均场模型。该模型在形式上基本等价于玻尔兹曼机(Boltzmann Machine)。SK模型源于Sherrington和Kirkpatrick为了使平均场理论有效而对哈密顿量的修改。

这里讨论SK模型是为了阐述使用副本方法(Replica Method)的基本计算方法。因为尽管SK模型很简单，但它的理论分析非常丰富，涉及连续的现象、序参量、状态的超度量空间和随机分支过程等。

理论模型

我们先从哈密顿量开始，SK 模型的哈密顿量为：
\begin{equation}
H = - \sum_{i<j}J_{ij}S_{i}S_{j} - h \sum_{i}S_{i}. \label{equ:skhamilton}
\end{equation}
其中，$S_{i} \in \{-1,+1\}$，在 Ising 模型中表示微观磁矩的指向是上方还是下方；第一项求和遍及$N$格点格子中的所有自旋对，但并不双重计数；$h$ 是均匀磁场；耦合常数$J_{ij}\sim\mathcal{N}(\frac{J_{0}}{N}, \frac{J^{2}}{N})$为高斯随机变量：

\begin{equation}
P(J_{ij}) = \frac{1}{J} \sqrt{\frac{N}{2\pi}} e^{-\frac{N}{2J^{2}} (J_{ij}-\frac{J_{0}}{N})^{2}}.
\end{equation}
宽度$J$与位置$i,j$的距离无关。注意到加入非零的相互作用均值$J_0$不会带来其他问题，只会让公式更复杂。这里对上式进行简化：

\begin{equation}\label{equ:JProb}
P(J_{ij}) = \frac{1}{J} \sqrt{\frac{N}{2\pi}} e^{-\frac{NJ_{ij}^2}{2J^{2}}}. \tag{1}
\end{equation}

对于无限范围内平移不变的相互作用，平均场理论是严格的，因此我们期望有一种平均场类型的理论能够对这个模型进行有效求解。

由于这里的无序是淬火¹的，因此我们必须计算自由能并按照概率分布函数式\eqref{equ:JProb}对其求期望，即：
\begin{equation}
f=-\frac{k_BT}{N}\prod_{i<j}^N\mathbb{E}_{J_{ij}}\left[\ln Z(\{J_{ij}\},h,T)\right]=-\frac{k_BT}{N}\prod_{i<j}^N\int {\mathrm d} J_{ij}P(J_{ij})\ln Z(\{J_{ij}\},h,T).
\end{equation}
其中$k_B$为玻尔兹曼常数，且有：
\begin{equation}\label{equ:skpartial}
\begin{aligned}
Z & = \underbrace{\sum_{S_1=\pm 1}\cdots\sum_{S_N=\pm 1}}_\text{N重求和} \exp\left\{ \beta \left[ \sum_{i<j}^N J_{ij}S_{i}S_{j} + h \sum_{i=1}^N S_{i}\right] \right\}\\
& = \sum_{S_1=\pm 1}\cdots\sum_{S_N=\pm 1}\exp\left\{ \beta \left[ \sum_{i<j}^N J_{ij}S_{i}S_{j} + \frac{h}{2(N-1)} \sum_{i<j}^N(S_i+S_j)\right] \right\}\\
& = \sum_{S_1=\pm 1}\cdots\sum_{S_N=\pm 1}\prod_{i<j}^N \exp\left\{ \beta \left[ J_{ij}S_{i}S_{j} + \frac{h}{2(N-1)}(S_i+S_j)\right] \right\}.
\end{aligned}\tag{2}
\end{equation}

¹. 如果我们从冷却出发，系统如果冷却得很快称为淬火，如果冷却得慢称为退火。如果我们从无序出发，如果一个系统依赖于不随时间演化的随机变量，那么它称为淬火无序；如果它依赖于随时间演化的随机变量，那么它称为退火无序。淬火无序的一个具体例子即是当温度变化时，磁性原子分布不会变化。 ↩

基于Replica Method的求解框架

这里我们可以使用狄拉克符号，引入矩阵$\boldsymbol{P}$，其矩阵元:
\begin{equation}
\langle S_i|\boldsymbol{P}|S_j\rangle = \exp\left\{ \beta \left[ J_{ij}S_{i}S_{j} + \frac{h}{2(N-1)}(S_i+S_j)\right] \right\}.
\end{equation}

注意到$S_i,S_j$只能取值$\pm 1$，因此矩阵$\boldsymbol{P}$为如下的$2\times 2$的矩阵：
\begin{equation}
\begin{aligned}
\boldsymbol{P} &=
\begin{pmatrix}\langle 1|\boldsymbol{P}|1\rangle & \langle 1|\boldsymbol{P}|-1\rangle\\
\langle -1|\boldsymbol{P}|1\rangle & \langle -1|\boldsymbol{P}|-1\rangle
\end{pmatrix}\\
&= \begin{pmatrix} e^{\beta J_{ij}+h\beta/(N-1)} & e^{-\beta J_{ij}}\\
e^{-\beta J_{ij}} & e^{\beta J_{ij}-h\beta/(N-1)}
\end{pmatrix}.
\end{aligned}
\end{equation}

那么配分函数\eqref{equ:skpartial}式，可以改写为²：
\begin{equation}
\begin{aligned}
Z & = \sum_{S_1=\pm 1}\cdots\sum_{S_N=\pm 1}\prod_{i<j}^N \langle S_i|\boldsymbol{P}|S_j\rangle\\
& =\text{Tr}_{\{S_i\}}\exp[-\beta H(\{J_{ij}\},\{S_{i}\})].
\end{aligned}
\end{equation}
事实上，计算至此，由于随机变量$J_{ij}$的存在，已经无法继续计算。为了交换对于耦合常数$J_{ij}$求期望和对所有自旋位形$\{S_{i}\}$求迹的顺序，我们才使用replica方法。³
因此，配分函数的$n$次幂为：
\begin{equation}
Z^n = \prod_{\alpha=1}^n\text{Tr}_{\{S_i^\alpha\}}\exp[-\beta H(\{J_{ij}\},\{S_i^\alpha\})] = \text{Tr}_{\{S_i^\alpha\}}\exp\left\{ \beta\sum_{\alpha=1}^n\left[ \sum_{i<j}J_{ij}S_i^\alpha S_j^\alpha + h \sum_{i}S_i^\alpha\right] \right\}.
\end{equation}

注意到⁴：
\begin{equation}
\mathbb{E}_{J_{ij}}\left[\exp(\beta J_{ij}\sum_{\alpha=1}^nS_i^\alpha S_j^\alpha)\right]=\exp\left[ \beta^2J^2/2N\sum_{\alpha,\gamma=1}^n S_i^\alpha S_j^\alpha S_i^\gamma S_j^\gamma\right].
\end{equation}

那么显然：
\begin{equation}
\mathbb{E}_{}\left[Z^n\right] = \text{Tr}_{\{S_i^\alpha\}}\exp\left[ \frac{\beta^2J^2}{2N}\sum_{i<j}^N \sum_{\alpha,\gamma=1}^n S_i^\alpha S_j^\alpha S_i^\gamma S_j^\gamma + \beta h \sum_{i=1}^N\sum_{\alpha=1}^n S_i^\alpha\right].
\end{equation}

我们得到了不同副本之间的自旋的有效相互作用为：
\begin{equation}\label{equ:freeengrep}
\begin{aligned}
-\beta f &= \lim_{n\to 0}\frac1n[\mathbb{E}_{}[Z^n]-1]\\
&=\lim_{n\to 0, N\to\infty}\frac{1}{nN}\left\{\text{Tr}_{\{S_i^\alpha\}}\exp\left[ \frac{\beta^2J^2}{2N}\sum_{i<j}^N \sum_{\alpha,\gamma=1}^n S_i^\alpha S_j^\alpha S_i^\gamma S_j^\gamma + \beta h \sum_{i=1}^N\sum_{\alpha=1}^n S_i^\alpha\right]-1\right\}\\
&\approx \lim_{n\to 0, N\to\infty}\frac{1}{nN}\left[\text{Tr}_{\{S_i^\alpha\}}\exp(\mathcal{H})-1\right].
\end{aligned}
\tag{3}
\end{equation}
事实上，至此，通用的计算就已经结束。下面至文末的计算是基于某一个具体模型的各种分析技巧，除了正交解耦技巧外，其他不具备普适性。

². 如果$J_{ij}$回退Ising模型，视为常数而非随机变量，那么： ↩

\begin{equation}
\begin{aligned}
Z & = \sum_{S_1=\pm 1}\cdots\sum_{S_N=\pm 1}\prod_{i<j}^N \langle S_i|\boldsymbol{P}|S_j\rangle\\
& =\sum_{S_1}\langle S_1|\boldsymbol{P}^{N(N-1)/2}|S_1\rangle\\
& = \text{Tr}\boldsymbol{P}^{N(N-1)/2}.
\end{aligned}
\end{equation}
这里求矩阵的迹Tr是严格的，下文中不严格的复用了这个表示，指代对所有可能求和。

³. 这便是自平均性的体现，当$N\gg \beta J\gg 1$时，以$1/N$为基展开，则有： ↩

\begin{equation}
\mathbb{E}_{}[Z^n] = \mathbb{E}_{}[Z]^n.
\end{equation}

⁴. 这里使用恒等式： ↩

\begin{equation}\label{equ:guassianintreq}
\sqrt{\frac{a}{2\pi}}\int_{\mathbb{R}}e^{-\frac12(ax^2-\lambda\sqrt{2}x)}{\mathrm d} x = e^{\lambda^2/4a}. \tag{4}
\end{equation}

继续求解

再次注意到：
\begin{equation}
\begin{aligned}
\sum_{i<j}^N \sum_{\alpha,\gamma=1}^n S_i^\alpha S_j^\alpha S_i^\gamma S_j^\gamma &= \frac12\sum_{i<j}^N \sum_{\alpha\neq\gamma}^n S_i^\alpha S_j^\alpha S_i^\gamma S_j^\gamma + \frac{nN(N-1)}{2} \\
& = \frac12\sum_{\alpha\neq\gamma}^n\left(\sum_{i=1}^N S_i^\alpha S_i^\gamma\right)^2+\frac{nN(N-1)}{2}.
\end{aligned}
\end{equation}

因此有：
\begin{equation}
\mathcal{H} = \frac{\beta^2J^2nN}{4} + \frac{\beta^2J^2}{2N}\sum_{\alpha,\gamma=1}^n\left(\sum_{i=1}^N S_i^\alpha S_i^\gamma\right)^2 + \beta h\sum_{i=1}^N\sum_{\alpha=1}^n S_i^\alpha.
\end{equation}
其中，副本的相互总用仅在第二项中出现。由于$\mathcal{H}$中不同位置的自旋总是耦合在一起的，因此，式\eqref{equ:freeengrep}中的求迹运算极端复杂(矩阵不是对角阵)。我们可以使用对角化技术，也即引入一组虚拟变量来消除这些耦合。对式\eqref{equ:guassianintreq}中的$a=N,\lambda = \sqrt{2}\beta J\sum_{i=1}^n S_i^\alpha S_i^\gamma$，有恒等式：
\begin{equation}
\sqrt{\frac{N}{2\pi}}\int_{\mathbb{R}}e^{-\frac12Nx_{\alpha\gamma}^2-\beta J\left(\sum_{i=1}^N S_i^\alpha S_i^\gamma\right) x_{\alpha\gamma} }{\mathrm d} x_{\alpha\gamma} = e^{\beta^2J^2\left(\sum_{i=1}^N S_i^\alpha S_i^\gamma\right)^2/2N}.
\end{equation}

这直接导致：
\begin{equation}
\begin{aligned}
e^{\mathcal{H}} & = \exp\left[\frac{\beta^2J^2nN}{4} + \beta h\sum_{i=1}^N\sum_{\alpha=1}^n S_i^\alpha\right]\prod_{\alpha,\gamma=1}^n \sqrt{\frac{N}{2\pi}}\int_{\mathbb{R}}e^{-\frac12Nx_{\alpha\gamma}^2-\beta J\left(\sum_{i=1}^N S_i^\alpha S_i^\gamma\right) x_{\alpha\gamma} }{\mathrm d} x_{\alpha\gamma}
\end{aligned}
\end{equation}
中不同位置的自旋不再耦合在一起，并且上式中求迹就化为对一个位置上(例如$i$)的自旋求迹，并有$N$项乘积。因此，忽略位置标记$i$，研究不同副本之间的自旋耦合上。

这里使用最速下降法(Method of Steepest Descent，细节可以参考上一篇博文)，处理上述积分，可以得到：
\begin{equation}
\text{Tr}e^{\mathcal{H}}=\exp\left[\frac{\beta^2J^2N}{4}[n-n(n-1)q^2]\right] \left[\text{Tr}\exp\left\{ \beta h\sum_{\alpha=1}^nS^\alpha + \beta^2J^2q\sum_{\alpha,\gamma=1}^nS^\alpha S^\gamma\right\}\right]^N
\end{equation}
其中$q$为该模型的序参量⁵。

再次使用式\eqref{equ:guassianintreq}技巧解耦自旋：
\begin{equation}
\exp\left[ \beta^2J^2q\sum_{\alpha,\gamma=1}^nS^\alpha S^\gamma \right] = \exp\left[-\frac{n\beta^2J^2q}{2}\right]\int_{\mathbb{R}}e^{-z^2/2}\exp\left[z\beta \sqrt{q}\sum_{\alpha=1}^nS^\alpha\right]\frac{\operatorname{d} z}{\sqrt{2\pi}},
\end{equation}
得到：
\begin{equation}
\begin{aligned}
\text{Tr}e^{\mathcal{H}}&=e^P\left[\text{Tr}\int_{\mathbb{R}}e^{-z^2/2}\exp\left[(z\beta J\sqrt{q} + \beta h)\sum_{\alpha=1}^nS^\alpha\right]\frac{\operatorname{d} z}{\sqrt{2\pi}}\right]^N\\
&=e^P\left[\int_{\mathbb{R}}e^{-z^2/2}(2\cosh\mathcal{Z})^n\frac{\operatorname{d} z}{\sqrt{2\pi}}\right]^N\\
& = e^P\exp\left[N\ln\int_{\mathbb{R}}e^{-z^2/2}(2\cosh\mathcal{Z})^n\frac{\operatorname{d} z}{\sqrt{2\pi}}\right]\\
&\approx e^P\exp\left[N\ln\int_{\mathbb{R}}e^{-z^2/2}(1+n\ln (2\cosh\mathcal{Z}))\frac{\operatorname{d} z}{\sqrt{2\pi}}\right]\\
&\approx e^P\exp\left[nN\ln\int_{\mathbb{R}}e^{-z^2/2}\ln(2\cosh\mathcal{Z})\frac{\operatorname{d} z}{\sqrt{2\pi}}\right].
\end{aligned}
\end{equation}
其中$P=N\beta^2J^2(n-n(n-1)q^2)/4-nN\beta^2J^2q/2$，$\mathcal{Z} = z\beta J\sqrt{q}+\beta h$。将上述结果带入\eqref{equ:freeengrep}式，并取极限：
\begin{equation}
\beta f = -\frac{\beta^2J^2}{4}(1-q)^2-\int_{\mathbb{R}}e^{-z^2/2}\ln(2\cosh\mathcal{Z})\frac{\operatorname{d} z}{\sqrt{2\pi}}.
\end{equation}
至此，理论计算求解结束。

⁵. 该序参量由下式给出： ↩

\begin{equation}
\begin{aligned}
q = \langle\langle S^\alpha S^\gamma\rangle\rangle &= \left. \frac{\text{Tr}S^\alpha S^\gamma e^{\mathcal{H}}}{\text{Tr}e^{\mathcal{H}}}\right|_{n\to 0}\\
& = \left.\frac{\int_{\mathbb{R}}e^{-z^2/2}(2\sinh\mathcal{Z})^2(2\cosh\mathcal{Z})^{n-2}\frac{\operatorname{d} z}{\sqrt{2\pi}}}{\int_{\mathbb{R}}e^{-z^2/2}(2\cosh\mathcal{Z})^n\frac{\operatorname{d} z}{\sqrt{2\pi}}}\right|_{n\to 0}\\
& = \int_{\mathbb{R}}e^{-z^2/2}\tanh^2\mathcal{Z}\frac{\operatorname{d} z}{\sqrt{2\pi}}.
\end{aligned}
\end{equation}