% % 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,upgreek} \usepackage{algcompatible} \usepackage{algorithm} % 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\tr{\mathop{\mathrm{Tr}}\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} \def\X{\mathbb X} \def\inf{\mathrm{inf}} % 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}}} \def\mcmc{\textsc{McMC}} \def\vsigma{\pmb{\upsigma}} % used to reformat display math to fit in two-column ``reprint'' mode \makeatletter \newif\ifreprint \@ifclasswith{revtex4-1}{reprint}{\reprinttrue}{\reprintfalse} \makeatother \begin{document} \title{An efficient cluster algorithm for spin systems in a symmetry-breaking 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} \end{abstract} \maketitle Let $G=(V,E)$ be a graph, where the set of vertices $V=\{1,\ldots,N\}$ enumerates the sites of a lattice and the set of edges $E$ contains pairs of neighboring sites. Let $R$ be a group and $X$ an $R$-set, with the action of group elements $r\in R$ on elements $s\in X$ denoted $r\cdot s$. $X$ is the set of states accessible by a spin, and $R$ is the \emph{symmetry group} of $X$. The set $X$ must admit a measure $\mu$ that is invariant under the action of $R$, e.g., for any $A\subseteq X$ and $r\in R$, $\mu(r\cdot A)=\mu(A)$. This trait is shared by the counting measure on any discrete set, or by any group acting by isometries on a Riemannian manifold, such as $O(n)$ on $S^{n-1}$ in the $n$-component model. Finally, the subset of elements in $R$ of order two must act transitively on $X$. This property, while apparently obscure, is shared by any symmetric space \cite{loos1969symmetric} or by any transitive, finitely generated isometry group. In fact, all the examples listed here have spins spaces with natural metrics whose symmetry group is the set of isometries of the spin spaces. We put one spin at each site of the lattice described by $G$, so that the state of the entire spin system is described by elements $\vec s\in X\times\cdots\times X=X^N$. The Hamiltonian of this system is a function $\H:X^N\to\R$ defined by \[ \H(\vec s)=-\!\!\!\!\sum_{\{i,j\}\in E}\!\!\!\!Z(s_i,s_j)-\sum_{i\in V}B(s_i), \] where $Z:X\times X\to\R$ couples adjacent spins and $B:X\to\R$ is an external field. $Z$ must be symmetric in its arguments and invariant under the action of any element of $R$ applied to the entire lattice, that is, for any $r\in R$ and $s,t\in X$, $Z(r\cdot s,r\cdot t)=Z(s,t)$. One may also allow $Z$ to also be a function of the edge---for modelling random-bond, long-range, or anisotropic interactions---or allow $B$ to be a function of site---for modelling random fields. All the formal results of this paper hold equally well for these cases, but we will drop the additional index notation for clarity. \begin{table*}[htpb] \begin{tabular}{l||ccccc} & Spin states ($X$) & Symmetry ($R$) & Action ($g\cdot s$) & Coupling ($Z(s,t)$) & Field ($B(s)$) \\ \hline\hline Ising & $\{-1,1\}$ & $\Z/2\Z$ & $0\cdot s\mapsto s$, $1\cdot s\mapsto -s$ & $st$ & $Hs$ \\ $n$-component & $S^{n-1}$ & $\mathrm O(n)$ & $M\cdot s\mapsto Ms$ & $s^{\mathrm T}t$ & $H^{\mathrm T}s$\\ Potts & $\mathbb Z/q\mathbb Z$ & $D_n$ & $r_m\cdot s=m+s$, $s_m\cdot s=m-s$ & $\delta(s,t)$ & $\sum_mH_m\delta(m,s)$\\ Clock & $\mathbb Z/q\mathbb Z$ & $D_n$ & $r_m\cdot s=m+s$, $s_m\cdot s=-m-s$ & $\cos(2\pi\frac{s-t}q)$ & $\sum_mH_m\cos(2\pi\frac{s-m}q)$\\ Roughening & $\mathbb Z$ & $D_\inf$ & $r_m\cdot s=m+s$, $s_m\cdot s=-m-s$ & $(s-t)^2$ & $Hs^2$\\ \end{tabular} \caption{Some models.} \label{table:models} \end{table*} The goal of statistical mechanics as applied to these systems is to compute expectation values of observables $A:X^N\to\R$. Assuming the ergodic hypothesis holds (for systems with broken-symmetry states, it does not), the expected value $\avg A$ of an observable $A$ is its average over every state $\vec s$ in the configuration space $X^N$ weighted by the probability $p(\vec s)$ of that state appearing, or \[ \avg A =\frac{\int_{X^N}A(\vec s)e^{-\beta\H(\vec s)}\,\dd\mu(\vec s)} {\int_{X^N}e^{-\beta\H(\vec s)}\,\dd\mu(\vec s)} \] We will first describe a generalized version of the celebrated Wolff algorithm in the standard case where $B(s)=0$. After reflecting on the technical requirements of that algorithm, we will introduce a transformation to our system and Hamiltonian that allows the same algorithm to be applied with nonzero, in fact \emph{arbitrary}, external fields. The Wolff algorithm proceeds in the following way. \begin{enumerate} \item Pick a random site and a random rotation $r\in R$ of order two, and add the site to a stack. \item While the stack isn't empty, \begin{enumerate} \item pop site $m$ from the stack. \item If site $m$ isn't marked, \begin{enumerate} \item mark the site. \item For every $j$ such that $\{m,j\}\in E$, add site $j$ to the stack with probability \[ p_r(s_m,s_j)=\min\{0,1-e^{\beta(Z(r\cdot s_m,s_j)-Z(s_m,s_j))}\}. \] \item Take $s_m\mapsto r\cdot s_m$. \end{enumerate} \end{enumerate} \end{enumerate} When the stack is exhausted, a cluster of connected spins will have been rotated by the action of $r$. In order for this algorithm to be useful, it must satisfy ergodicity and detailed balance. The probability $P(\vec s\to\vec s')$ that the configuration $\vec s$ is brought to $\vec s'$ by the flipping of a cluster $C\subseteq E$ is related to the probability of the reverse process $P(\vec s'\to\vec s)$ by \begin{widetext} \[ \begin{aligned} \frac{P(\vec s\to\vec s')}{P(\vec s'\to\vec s)} &=\prod_{\{i,j\}\in C}\frac{p_r(s_i,s_j)}{p_{r^{-1}}(s_i',s_j')}\prod_{\{i,j\}\in\partial C}\frac{1-p_r(s_i,s_j)}{1-p_{r^{-1}}(s'_i,s'_j)} =\prod_{\{i,j\}\in C}\frac{p_r(s_i,s_j)}{p_{r}(r\cdot s_i,r\cdot s_j)}\prod_{\{i,j\}\in\partial C}\frac{1-p_r(s_i,s_j)}{1-p_{r^{-1}}(r\cdot s_i,s_j)} \\ &=\prod_{\{i,j\}\in C}\frac{p_r(s_i,s_j)}{p_{r}(s_i,s_j)}\prod_{\{i,j\}\in\partial C}e^{\beta(Z(r\cdot s_i,s_j)-Z(s_i,s_j))} =\frac{e^{-\beta\H(\vec s)}}{e^{-\beta\H(\vec s')}} \end{aligned} \] \end{widetext} which satisfies detailed balance so long as $p_r(s,t)=p_{r^{-1}}(s,t)$. Therefore, for the algorithm to function one must choose rotations at random only from those group elements in $R$ that act by idempotents. Usually, as is the case in all the examples here, the action is natural enough that this is the same as only choosing group elements of order two to act on your spins. Ergodicity is achieved if the subset of symmetry rotations that obey this property allow any spin state to access any other under composition. This is true for all the models in Table \ref{table:models}. The function of the algorithm described above depends on the fact that the coupling $Z$ depends only on the relative orientation of the spins---global reorientations by acting by some rotation do not affect the Hamiltonian. The external field $B$ breaks this symmetry. However, this can be resolved. Define a new graph $\tilde G=(\tilde V,\tilde E)$, where $\tilde V=\{0,1,\ldots,N\}$ and \[ \tilde E=E\cup\big\{\{0,i\}\mid i\in V\big\}. \] We have introduced a new site to the lattice that neighbors every other site. Instead of assigning this site a spin whose value comes from the set $X$, we will assign it values $s_0\in R$, symmetry group elements, so that the new configuration space of the model is $R\times X^N$. We introduce a Hamiltonian $\tilde\H:R\times X^N\to\R$ defined by \[ \begin{aligned} \tilde\H(s_0,\vec s) &=-\!\!\!\!\sum_{\{i,j\}\in E}\!\!\!\!Z(s_i,s_j) -\sum_{i\in V}B(s_0^{-1}\cdot s_i)\\ &=-\!\!\!\!\sum_{\{i,j\}\in\tilde E}\!\!\!\!\tilde Z(s_i,s_j) \end{aligned} \] where the new coupling $\tilde Z:(R\cup X)\times(R\cup X)\to\R$ is defined for $s,t\in R\cup X$ by \[ \tilde Z(s,t) = \begin{cases} Z(s,t) & \text{if $s,t\in X$} \\ B(s^{-1}\cdot t) & \text{if $s\in R$} \\ B(t^{-1}\cdot s) & \text{if $t\in R$} \end{cases} \label{eq:new.z} \] Note that this modified coupling is invariant under the action of group elements: for any $r,s_0\in R$ and $s\in X$, \[ \begin{aligned} \tilde Z(rs_0,r\cdot s) &=B((rs_0)^{-1}\cdot (r\cdot s))\\ &=B((s_0^{-1}r^{-1})\cdot(r\cdot s))\\ &=B((s_0^{-1}r^{-1}r)\cdot s)\\ &=B(s_0^{-1}\cdot s) =\tilde Z(s_0,s) \end{aligned} \] The invariance $\tilde Z$ to rotations given other arguments follows from the invariance properties of $Z$. We have produced a system that incorporates the field function $B$ whose Hamiltonian is invariant to global rotations, but how does it relate to our previous system, whose properties we actually want to measure? If $A:X^N\to\R$ is an observable of the original system, one can construct an observable $\tilde A:R\times X^N\to\R$ of the new system defined by \[ \tilde A(s_0,\vec s)=A(s_0^{-1}\cdot\vec s) \] whose expectation value in the new system equals that of the original observable in the old system. First, note that $\tilde\H(1,\vec s)=\H(\vec s)$. Since the Hamiltonian is invarient under global rotations, it follows that for any $g\in R$, $\tilde\H(g,g\cdot\vec s)=\tilde\H(g^{-1}g,g^{-1}g\cdot\vec s)=\tilde\H(1,\vec s)=\H(\vec s)$. Using the invariance properties of the measure on $X$ and introducing a measure $\rho$ on $R$, it follows that \[ \begin{aligned} \avg{\tilde A} &=\frac{ \int_R\int_{X^N}\tilde A(s_0,\vec s)e^{-\beta\tilde\H(s_0,\vec s)}\,\dd\mu(\vec s)\,\dd\rho(s_0) } { \int_R\int_{X^N}e^{-\beta\tilde\H(s_0,\vec s)}\,\dd\mu(\vec s)\,\dd\rho(s_0) }\\ &=\frac{ \int_R\int_{X^N}A(s_0^{-1}\cdot\vec s)e^{-\beta\tilde\H(s_0,\vec s)}\,\dd\mu(\vec s)\,\dd\rho(s_0) } { \int_R\int_{X^N}e^{-\beta\tilde\H(s_0,\vec s)}\,\dd\mu(\vec s)\,\dd\rho(s_0) }\\ &=\frac{ \int_R\int_{X^N}A(\vec s')e^{-\beta\tilde\H(s_0,s_0\cdot\vec s')}\dd\mu(s_0\cdot\vec s')\,\dd\rho(s_0) } { \int_R\int_{X^N}e^{-\beta\tilde\H(s_0,s_0\cdot\vec s')}\dd\mu(s_0\cdot\vec s')\,\dd\rho(s_0) }\\ &=\frac{ \int_R\dd\rho(s_0)}{ \int_R\dd\rho(s_0)}\frac{\int_{X^N}A(\vec s')e^{-\beta\H(\vec s')}\dd\mu(\vec s') }{\int_{X^N}e^{-\beta\H(\vec s')}\dd\mu(\vec s') } =\avg A \end{aligned} \] To summarize, spin systems in a field may be treated in the following way. \begin{enumerate} \item Add a site to your lattice adjacent to every other site. \item Initialize a ``spin'' at that site that is a representation of a member of the symmetry group of your ordinary spins. \item Carry out the ordinary Wolff cluster-flip procedure on this new lattice, substituting $\tilde Z$ as defined in \eqref{eq:new.z} for $Z$. \end{enumerate} Ensemble averages of observables $A$ can then be estimated by sampling the value of $\tilde A$ on the new system. \section{Examples} \subsection{The Ising Model} In the Ising model, spins are drawn from the set $\{1,-1\}$. The symmetry group of this model is $C_2$, the cyclic group on two elements, which can be conveniently represented by the multiplicative group with elements $\{1,-1\}$, exactly the same as the spins themselves. The only nontrivial element is of order two. Because the symmetry group and the spins are described by the same elements, performing the algorithm on the Ising model in a field is very accurately described by simply adding an extra spin coupled to all others and running the ordinary algorithm. \subsection{The $n$-component Model} In the $n$-component model, spins are described by vectors on the $(n-1)$-sphere, so that $X=S^{n-1}$. The symmetry group of this model is $O(n)$, $n\times n$ orthogonal matrices. The symmetry group acts on the spins by matrix multiplication. The elements of $O(n)$ that are order two are reflections about some hyperplane through the origin and $\pi$ rotations about any axis through the origin. Since the former generate the entire group, the set of reflections alone suffices to provide ergodicity. Computation of the coupling of ordinary spins with the external field and expectation values requires a matrix inversion, but since the matrices in question are orthogonal this is quickly accomplished by a transpose. \subsection{The Potts \& Clock Models} In both the $q$-state Potts and clock models, spins are described by $\Z/q\Z$, the set of integers modulo $q$. The symmetry group of this model is the dihedral group $D_q=\{r_0,\ldots,r_{q-1},s_0,\ldots,s_{q-1}\}$, the group of symmetries of a regular $q$-gon. The element $r_n$ represents a rotation of the polygon by $2\pi n/q$, and the element $s_n$ represents a reflection composed with a rotation $r_n$. The group acts on the spins by permutation: $r_n\cdot m={n+m}\pmod q$ and $s_n\cdot m={-(n+m)}\pmod q$. Intuitively, this can be thought of as the natural action of the group on the vertices of a regular polygon that have been numbered $0$ through $q-1$. The elements of $D_q$ that are of order 2 are all reflections and $r_{q/2}$ if $q$ is even, though the former can generate the latter. While the reflections do not necessarily generate the entire group, for any $n,m\in\Z/q\Z$ there exists a reflection that takes $n\to m$, ensuring ergodicity. The elements of the dihedral group can be stored simply as an integer and a boolean that represents whether the element is a pure rotation or a reflection. The principle difference between the Potts and clock models is that, in the latter case, the form of the coupling $Z$ allows a geometric interpretation as being two-dimensional vectors fixed with even spacing along the unit circle. \subsection{Roughening} Though not often thought of as a spin model, simple roughening of surfaces can be described in this framework. The set of states is the integers $\Z$ and its symmetry group is the infinite dihedral group $D_\infty=\{r_i,s_i\mid i\in\Z\}$, where the action of the symmetry on the spins $j\in\Z$ is given by $r_i\cdot j=i+j$ and $s_i\cdot j=i-j$. These are shifts by $i$ and reflection about the integer $i$, respectively. The elements of order two are the reflections $s_i$, which suffice to provide ergodicity as any integer can be taken to any other in one step of this kind. The coupling is usually taken to be $Z(i,j)=|i-j|$, though it may also be any function of the absolute difference. Because randomly choices of integer will almost always result in energy changes so big that the whole system is always flipped, it is better to select random reflections about integers close to the average state of the system. Continuous roughening models---where the spin states are described by real numbers and the symmetry group is $\mathrm E(1)$, the Euclidean group for one-dimensional space---are equally well described. %\begin{figure} % \centering % \input{fig_correlation} % \caption{The autocorrelation time $\tau$ of the internal energy $\H$ for a $n=32\times32$ square-lattice Ising % model simulated by the three nonzero field methods detailed in this paper. % The top plot is for simulations at high temperature, $T=5.80225$, the middle %plot is for simulations at the critical temperature, $T=2.26919$, and the %bottom plot is for simulations at } % \label{fig:correlation} %\end{figure} \begin{figure} \centering \input{fig_correlation_collapse-hL} \caption{The correlation time $\tau$ as a function of the renormalization invarient $hL^{-\beta\delta/\nu}$ for the $N=L\times L$ square lattice Ising model for $L=8$, $16$, $32$, $64$, $128$, and $256$. $z=0.3$ } \end{figure} \begin{figure} \centering \input{fig_correlation-temp} \caption{The correlation time $\tau$ as a function of the renormalization invarient $ht^{-\beta\delta}$ for the $N=128\times 128$ square lattice Ising model. $z=0.3$ } \end{figure} \begin{acknowledgments} \end{acknowledgments} \bibliography{monte-carlo} \end{document}