path: root/monte-carlo.tex

%
%  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}