% % Created by Jaron Kent-Dobias on Thu Apr 20 12:50:56 EDT 2017. % Copyright (c) 2017 Jaron Kent-Dobias. All rights reserved. % \documentclass[aps,prl,reprint]{revtex4-1} \usepackage[utf8]{inputenc} \usepackage{amsmath,amssymb,latexsym,mathtools,xifthen} \usepackage{algorithm} \usepackage{algorithmic} % uncomment to label only equations that are referenced in the text %\mathtoolsset{showonlyrefs=true} % I want labels but don't want to type out ``equation'' \def\[{\begin{equation}} \def\]{\end{equation}} % math not built-in \def\arcsinh{\mathop{\mathrm{arcsinh}}\nolimits} \def\arccosh{\mathop{\mathrm{arccosh}}\nolimits} \def\ei{\mathop{\mathrm{Ei}}\nolimits} % exponential integral Ei \def\re{\mathop{\mathrm{Re}}\nolimits} \def\im{\mathop{\mathrm{Im}}\nolimits} \def\sgn{\mathop{\mathrm{sgn}}\nolimits} \def\dd{d} % differential \def\O{O} % big O \def\o{o} % little O \def\Z{\mathbb Z} \def\R{\mathbb R} \def\G{S} % subscript for ``critical'' values, e.g., T_\c \def\c{\mathrm c} % scaling functions \def\fM{\mathcal M} % magnetization \def\fX{\mathcal Y} % susceptibility \def\fF{\mathcal F} % free energy \def\fiF{\mathcal H} % imaginary free energy \def\fS{\mathcal S} % surface tension \def\fG{\mathcal G} % exponential factor \def\H{\mathcal H} % lattice types \def\sq{\mathrm{sq}} \def\tri{\mathrm{tri}} \def\hex{\mathrm{hex}} % dimensions \def\dim{d} \def\twodee{\textsc{2d} } \def\threedee{\textsc{3d} } \def\fourdee{\textsc{4d} } % fancy partial derivative \newcommand\pd[3][]{ \ifthenelse{\isempty{#1}} {\def\tmp{}} {\def\tmp{^#1}} \frac{\partial\tmp#2}{\partial#3\tmp} } \newcommand\nn[2]{\langle#1#2\rangle} \newcommand\avg[1]{\langle#1\rangle} \newcommand\eavg[1]{\langle#1\rangle_{\mathrm e}} \newcommand\mavg[1]{\langle#1\rangle_{\mathrm m}} \def\e{\mathrm e} \def\m{\mathrm m} \newcommand\new[1]{\tilde{#1}} \renewcommand\vec[1]{\mathbf{#1}} \newcommand\unit[1]{\hat{\mathbf{#1}}} % used to reformat display math to fit in two-column ``reprint'' mode \makeatletter \newif\ifreprint \@ifclasswith{revtex4-1}{reprint}{\reprinttrue}{\reprintfalse} \makeatother \begin{document} \algsetup{indent=1em} \title{A cluster algorithm for the Ising model in an external field} \author{Jaron Kent-Dobias} \author{James P.~Sethna} \affiliation{Laboratory of Atomic and Solid State Physics, Cornell University, Ithaca, NY, USA} \date\today \begin{abstract} An abstract. \end{abstract} \maketitle \section{Introduction} The Ising model is a simple model of a magnet. It is comprised of $n$ sites arranged in a regular lattice. At every site $i$ is a magnetic spin that can be oriented up or down, corresponding to a local magnetization $s_i$ with value $1$ or $-1$, respectively. Configurations of the spins are uniquely described by the local magnetization of every site, so the model's configuration space is $\G^n=\G\times\cdots\times\G$ for $\G=\{-1,1\}$ and states of the model are given by $s=(s_1,\ldots,s_n)\in\G^n$. The Ising Hamiltonian $\H:\G^n\to\R$ gives the energy of a particular state $s\in\G^n$ by \[ \H(s)=-\frac12\sum_{i=1}^n\sum_{j=1}^nJ_{ij}s_is_j-HM(s) \label{eq:ham} \] where $H$ is the external magnetic field, $M:\G^n\to\R$ is the magnetization and is defined for $s\in\G^n$ by \[ M(s)=\sum_{i=1}^ns_i, \] and the spin--spin coupling $J$ is defined by \[ J_{ij}=\begin{cases} 1 & \text{$i$, $j$ are neighboring sites}\\ 0 & \text{otherwise.} \end{cases} \] The quantities of interest in statistical physics are averaged out over timescales large compared to the dynamics of the system. If $A:\G^n\to\R$ is a physical quantity that depends on the model's configuration and $s(t)$ gives the configuration as a function of time $t$, then \[ \bar A=\lim_{T\to\infty}\frac1T\int_0^TA(s(t))\,\dd t \] Typically the average value of physical quantities is not actually computed by taking a time average. Instead, the ergodic hypothesis is applied, which stipulates that the system's dynamics access all states $s\in S^n$ with relative duration $e^{-\beta\H(s)}$, so that the time average can be written as the ensemble average, or the average over configuration space, \[ \avg A=\frac{\sum_se^{-\beta\H(s)}A(s)}{\sum_se^{-\beta\H(s)}} =\frac1Z\sum_se^{-\beta\H(s)}A(s), \label{eq:avg.obs} \] where $\beta=1/kT$ is the inverse temperature in energy units and \[ Z=\sum_se^{-\beta\H(s)} \] is the partition function. Monte Carlo methods attempt to estimate quantities of the form \eqref{eq:avg.obs} without doing the often intractable sum explicitly by sampling the set of possible states such that a state $s$ is sampled with probably $e^{-\beta\H(s)}$, appropriately weighting each state's contribution to the sum. A common way of doing this is the Markov-chain Monte Carlo approach, which steps through configuration space using a Markov process whose equilibrium distribution is Boltzmann. One of the most popular of these is the Metropolis--Hastings algorithm, shown in Algorithm \ref{alg:met-hast}. One spin in the current state $s$ is randomly chosen and perturbed. The new state is accepted if the resulting change in energy is negative, or with probability $e^{-\beta\Delta\H}$ if the change is positive. Thermodynamic averages of functions of the state can be made by periodically sampling from this loop for times much longer than the autocorrelation time of the Markov process. \begin{figure} \begin{algorithm}[H] \begin{algorithmic} \REQUIRE $s\in\G^n$ \LOOP \STATE $s'\gets s$ \STATE $i \gets\text{random integer between $1$ and $n$}$ \STATE $s'_i\gets-s_i$ \IF{$e^{-\beta(\H(s')-\H(s))}\geq\text{random number between $0$ and $1$}$} \STATE $s\gets s'$ \ENDIF \ENDLOOP \end{algorithmic} \caption{Metropolis--Hastings} \label{alg:met-hast} \end{algorithm} \end{figure} This algorithm is very generic and works well for many models. However, near the critical point of the Ising model it suffers from critical slowing down, where the autocorrelation time of the algorithm diverges with proximity to the critical point. The Wolff algorithm is such a method that is very efficient when $H=0$, i.e., in the absence of an external field. It is shown in Algorithm \ref{alg:wolff}. Briefly, whole clusters of spins are flipped in a way that correctly samples states from the Boltzmann distribution while never rejecting a flip attempt. A random spin is selected, then flipped. Each of the spin's neighbors that was initially oriented the same direction as the spin is then flipped with probability $1-e^{-2\beta}$, and the same is true recursively for each of those flipped spins' neighbors. The Wolff algorithm solves many of the problems of the Metropolis--Hastings algorithm, but only in the absence of a field, since it depends on the natural symmetry $\H(s)=\H(-s)$ which is true only when $H=0$. The Wolff algorithm can be extended to work in a field by combining it with aspects of the Metropolis--Hastings algorithm in what we will refer to as the hybrid Wolff/Metropolis--Hastings. In each step of this algorithm, a Wolff-style cluster is constructed exactly as above. Instead of moving on, however, the system is reverted to its previous state with probability $\max\{0,1-e^{\beta\Delta M H}\}$, where $\Delta M$ is the change in magnetization due to the cluster flip. This is precisely a Metropolis--Hastings step where the perturbation of $s$ to $s'$ is the Wolff cluster and, since the detailed balance with respect to the spin--spin energy is performed by the cluster creation, only the change in energy due to the external field is considered. For $H=0$ this reverts to the Wolff algorithm exactly. A strictly faster version of this algorithm is accomplished by evaluating the probability that the cluster-building is reverted at every step in its formation, so that large unlikely clusters abandoned before very much time is wasted in constructing them. This algorithm is given as follows. \begin{figure} \begin{algorithm}[H] \begin{algorithmic} \REQUIRE $s\in\G^n$ \STATE \textbf{let} $q$ be an empty stack \STATE \textbf{let} $i_0$ be the index of a random site \STATE $q.\mathrm{push}(i_0)$ \STATE $\sigma\gets s_{i_0}$ \WHILE {$q$ is not empty} \STATE $i\gets q.\mathrm{pop}$ \IF {$s_i=\sigma$} \FOR {$j=1$ \TO $n$} \IF {$1-e^{-2\beta J_{ij}s_is_j}\geq$ random number between 0 and 1} \STATE $q.\mathrm{push}(j)$ \ENDIF \ENDFOR \STATE $s_i\gets-s_i$ \ENDIF \ENDWHILE \end{algorithmic} \caption{Wolff (Zero-Field)} \label{alg:wolff} \end{algorithm} \end{figure} \begin{figure} \begin{algorithm}[H] \begin{algorithmic} \REQUIRE $s\in\G^n$ \STATE \textbf{let} $q_1$ and $q_2$ be empty stacks \STATE \textbf{let} $i_0$ be the index of a random site \STATE $q_1.\mathrm{push}(i_0)$ \STATE $\sigma\gets s_{i_0}$ \STATE stopped $\gets$ \FALSE \WHILE {$q_1$ is not empty} \STATE $i\gets q_1.\mathrm{pop}$ \IF {$s_i=\sigma$} \FORALL {$j>0$ such that $J_{ij}\neq0$} \IF {$s_is_j=1$ \AND $1-e^{-2\beta J_{ij}}\geq$ random number between 0 and 1} \STATE $q_1.\mathrm{push}(j)$ \ENDIF \ENDFOR \STATE $s_i\gets-s_i$ \STATE $q_2.\mathrm{push}(i)$ \IF {$1-e^{-2\beta\sigma H}>$ random number between 0 and 1} \STATE stopped $\gets$ \TRUE \STATE \textbf{break} \ENDIF \ENDIF \ENDWHILE \IF {stopped} \WHILE {$q_2$ is not empty} \STATE $i\gets q_2.\mathrm{pop}$ \STATE $s_i\gets-s_i$ \ENDWHILE \ENDIF \end{algorithmic} \caption{Hybrid Wolff/Metropolis--Hastings} \label{alg:h.wolff} \end{algorithm} \end{figure} \section{Cluster-flip in a Field} Consider a system whose state $s=\{s_0,s_1,\ldots,s_n\}$ $s\in\Z_2^{n+1}$ is specified by the value of $n+1$ boolean variables with the Hamiltonian \[ \new\H(s)=-\sum_{i=0}^n\sum_{j=0}^n\new J_{ij}s_is_j, \label{eq:new.ham} \] where \[ \new J_{ij}=\begin{cases} 1 & \text{$i,j>0$, $i$, $j$ are neighboring sites}\\ H & \text{$i=0$ or $j=0$}\\ 0 & \text{otherwise} \end{cases} \] gives the coupling. Since \eqref{eq:new.ham} is a discrete-spin Hamiltonian invariant under the transformation $s\to-s$, it directly admits analysis using the Wolff algorithm specified above. By rewriting the sums, however, we also see that \[ \new\H(s)=-\frac12\sum_{i=1}^n\sum_{j=1}^nJ_{ij}s_is_j-H\new M(s) \label{eq:mod_hamiltonian} \] where \[ \new M(s)=s_0\sum_{i=1}^ns_i \] This is precisely the same Hamiltonian as \eqref{eq:ham} but with a redefined magnetization whose sign depends on the state of the zeroth spin. Our algorithm for sampling the Ising model in a field, then, is simply the Wolff algorithm described above with $\tilde J_{ij}$ replacing $J_{ij}$. These are formally equivalent, since every state $s\in S$ corresponds to two states $\tilde s\in\tilde S$, $\tilde s=\pm(1,s)$, such that $\H(s)=\tilde\H(\tilde s)$. We therefore have \begin{align} \avg A_{\tilde S} &=\frac{\sum_{\tilde s}e^{-\beta\tilde\H(\tilde s)}A(s_0s)}{\sum_{\tilde s}e^{-\beta\tilde\H(\tilde s)}}\\ &=\frac{\sum_se^{-\beta\tilde\H(s_0,s)}A(s)+\sum_se^{-\beta\tilde\H(-s_0,-s)}A(s)}{\sum_se^{-\beta\tilde\H(1,s)}+\sum_se^{-\beta\tilde\H(-1,-s)}}\\ &=\avg A_S \end{align} $\tilde A(s_0,s)=A(s_0s)$. This is true for the mangetization $\tilde M$ as defined above, for instance. \begin{figure} \begin{algorithm}[H] \begin{algorithmic} \REQUIRE $s\in\G^n$ \STATE \textbf{let} $q$ be an empty stack \STATE \textbf{let} $i_0$ be the index of a random site \STATE $q.\mathrm{push}(i_0)$ \STATE $\sigma\gets s_{i_0}$ \WHILE {$q$ is not empty} \STATE $i\gets q.\mathrm{pop}$ \IF {$s_i=\sigma$} \FOR {$j=0$ \TO $n$} \IF {$1-\exp(-2\beta \tilde J_{ij}s_is_j)\geq$ random number between 0 and 1} \STATE $q.\mathrm{push}(j)$ \ENDIF \ENDFOR \STATE $s_i\gets-s_i$ \ENDIF \ENDWHILE \end{algorithmic} \caption{Wolff (Non-zero Field)} \label{alg:wolff-field} \end{algorithm} \end{figure} At sufficiently small $H$ both our algorithm and the hybrid Wolff--Metropolis reduce to the ordinary Wolff algorithm. At very large $H$, the hybrid Wolff--Metropolis behaves exactly like Metropolis, where almost only one spin is ever flipped at a time with probability $\sim e^{-2\beta H}$, since the energy is dominated by the field. We measured the correlation time $\tau$ of these three algorithms at various fields and temperatures. This was done using a batch mean. Time was measured as ``the number of spins that have either been attempted to flip or attempted to add to a cluster.'' $E_s=A_s+XB_s$ \[ \pd{\avg B}X=\beta\big(\avg{B^2}-\avg B^2\big) \] \[ \chi=\pd{\avg M}H=\beta\big(\avg{M^2}-\avg M^2\big) \] Since \eqref{eq:mod_hamiltonian} has the form of a simple (ferromagnetic) spin-spin Hamiltonian, its ensemble average can be taken using the Wolff algorithm while never rejecting a step. \section{Magnetization Estimator} For any finite system, the average magnetization at zero field is identically zero at all temperatures. However, people often want to use finite simulations to estimate the average magnetization of the Ising model in the thermodynamic limit, where due to a superselection principle there is one of two average magnetizations at zero field below the critical temperature. The ergodic principle that defends the esemble average replacing a time average is no longer valid, since once the system has taken a certain magnetization it has zero probability of ever changing direction. This is typically accomplished by using $\avg{|M|}$. But what is the best way to estimate the magnetization of a finite system in the presence of a field? $\tilde M$ is the true magnetization of the finite-size system. We can estimate the magnetization of system in the thermodynamic limit by enforcing the superselection principle: \[ S_\e=\{s\mid M(s)>0\} \] \[ \eavg{A}=\frac{\sum_{s\in S_\e}e^{-\beta\H(s)}A(s)}{\sum_{s\in S_\e}e^{-\beta\H(s)}}=\frac1{Z_\e}\sum_{s\in S_\e}e^{-\beta\H(s)}A(s) \] \[ S_\m=\{s\mid M(s)<0\} \] \[ \mavg{A}=\frac{\sum_{s\in S_\m}e^{-\beta\H(s)}A(s)}{\sum_{s\in S_\m}e^{-\beta\H(s)}}=\frac1{Z_\m}\sum_{s\in S_\m}e^{-\beta\H(s)}A(s) \] For $H=0$, $\eavg M=\avg{|M|}$, defending the canonical measure of the ferromagnetic phase magnetization at zero field. This can be seen first by identifying the bijection $s\to-s$ that maps $S_\e$ to $S_\m$. Then, since $\H(s)=\H(-s)$, $M(s)=-M(-s)$ when $H=0$, $Z_\e=Z_\m$, $Z=2Z_\e+Z_0$, and \begin{align} \avg{|M|} &=\frac1{2Z_\e+Z_0}\sum_{s\in S}e^{-\beta\H(s)}M(s)\\ &=\frac1{Z_\e+\frac12Z_0}\sum_{s\in S_\e}e^{-\beta\H(s)}M(s)\\ &=\frac1{1+\frac{Z_0}{2Z_\e}}\eavg M \end{align} At infinite temperature, $Z_0/Z_\e\simeq N^{-1/2}\sim L^{-1}$ for large $L$, $N$. At any finite temperature, especially in the ferromagnetic phase, $Z_0\ll Z_\e$ in a much more extreme way. If the ensemble average over only positive magnetizations can be said to converge to the equilibrium state of the Ising model in the thermodynamic limit, what of the average over only \emph{negative} magnetizations, ordinarily unphysical? Collapse of the magnetizations $\avg M$ and $\eavg M$ at low and high temperature \begin{acknowledgments} Thanks! \end{acknowledgments} %\bibliography{monte-carlo} \begin{itemize} \item metastable and stable state mixing at finite size? \item tracking the wolff algorithm cluster size versus correlation length \end{itemize} \end{document}