%% ================================================================
\chapter{Historical context}
\label{chap:history}

%% ----------------------------------------------------------------
\section{Theory}

In 1924, the physicist Satyendra Nath Bose\index{Bose,~S.N.} examined the
quantum statistics of \plainidx{photons}.  In 1925, collaborating with Bose,
Albert Einstein\index{Einstein,~A.} realized that the same could be done with
non-interacting massive particles.  He also discovered the
\plainidx{condensation} phenomenon:  a macroscopic occupation of the
(single-particle) \plainidx{ground state} of the external potential
\cite{bib:lssy}.  Moreover, Einstein predicted a \plainidx{critical
temperature} for the phenomenon.  This temperature was so low --- at the
nanokelvin scale --- that \plainidx{Bose-Einstein condensation} attracted
little interest in the physics community.

Feynman\index{Feynman,~R.P.} in 1953, along with Penrose and Onsager in 1956
\cite{bib:feynman,bib:penons}, developed the theoretical notion of
\plainidx{long permutation cycles} in the \plainidx{Feynman-Kac representation}
of the Bose gas.  Feynman claimed that long cycles correspond to Bose-Einstein
condensation.

Andr\'as S\"ut\H{o} referred to the existence of long permutation cycles as
\emphidx{cycle percolation}.  He proved in 1993 that BEC implies cycle
percolation in the non-ideal (interacting) gas\index{non-ideal gas}
\cite{bib:suto1}, and proved the converse in 2002 for the ideal
(non-interacting) gas\index{ideal gas}.  S\"ut\H{o} moreover proved in the 2002
paper that there are infinitely many macroscopic cycles in the condensation of
the non-ideal Bose gas.

For the ideal Bose gas, BEC is defined as the macroscopic occupation of the
single-particle ground state of the external potential.  For an interacting
Bose gas, Hamiltonian eigenfunctions do not factor and thus there are no
single-particle ground states.  BEC is carefully defined for interacting
systems \cite{bib:lssy} in terms of the largest eigenvalue of a
\plainidx{density-matrix operator}.  The 1983 work of Buffet and Pul\`e
\cite{bib:buffetpule} examines the macroscopic occupation of the \plainidx{zero
Fourier mode}.

%% ----------------------------------------------------------------
\section{Experiments}

Liquid helium was produced in the laboratory by Kammerlingh Onnes
\index{Onnes,~H.K.} in 1908; Fritz London\index{London,~F.} in 1938
\cite{bib:london} connected superfluidity of \plainidx{liquid helium} with
Bose-Einstein condensation.  Here, however, atoms of liquid helium are strongly
interacting --- they attract only weakly, due to helium being a noble gas, but
there are strong repulsive effects due to the high density of the liquid.
Thus, Einstein's non-interacting theory could not explain the phenomenon.

Several groups attempted during the 1990s to produce BECs in vapors of
spin-polarized hydrogen, but were not able to achieve low enough temperatures.
The group of Cornell and Wieman \cite{bib:boulder}, using hybrid cooling
methods, successfully brought \plainidx{rubidium} atoms to well below the
critical temperature and made numerous measurements on the resulting
condensates.  (The group received the 2001 Nobel prize in physics for this
work.)

Interested in BECs was sparked by this experimental success:  thousands of
papers, both theoretical and experimental, have been published on BECs in the
years since.  The work of Cornell and Wieman was of interest for several
reasons:

\begin{itemize}

\item Condensates were directly imaged.  Measurements were taken of temperature,
density, position, velocity, particle number, and the fraction of the
condensate occupying the ground state of the 3D harmonic trapping potential.

\item The method was able to vary temperature and density through wide ranges;
the \plainidx{condensate fraction} was varied from zero to 100 percent.

\item The gaseous rubidium condensate was \emphidx{weakly interacting} ---
permitting a perturbative analysis which liquid helium, with its strong
interactions, did not allow.  (Note in particular that the work of Ueltschi et
al.  is a weak-interaction theory; it is valid to first order in the
interparticle interaction strength.)

\end{itemize}

%% ----------------------------------------------------------------
\section{Critical temperature}
\label{sec:critical_temperature}

Recall that Einstein predicted a critical temperature $T_c$ for the ideal Bose
gas.  It is a long-standing question to discover the effects of interbosonic
\plainidx{interaction strength} \plainidx{$a$} on the critical temperature.
Moreover, one may fix the density $\rho_c{(a)}$ and obtain a critical
temperature $T_c{(a)}=1/\beta_c{(a)}$ or vice versa; also, both of these
critical parameters depend on the interaction strength $a$.  One expects the
critical combination of parameters to be a manifold in $(\rho,\beta,a)$ space.
(See figure \ref{fig:critical_manifold}.)

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics{figures/critical_manifold.eps}
\caption{Critical manifold in $(\rho,\beta,a)$ for small $a$.
	\label{fig:critical_manifold}}
\end{center}
\end{figure}

Much is known about the $a=0$ line of this critical manifold; off $a=0$, even
the crude shape has been under debate.  The following findings are described by
\cite{bib:baymetal}:  The superfluid transition temperature of liquid helium is
lower than that of an ideal gas of the same density.  Thus, assuming that
helium superfluidity is a strongly interacting BEC, one would expects
interactions to decrease the critical temperature for the strongly interacting
case.  Various theoretical work (tabulated below) suggested either an increase
or a decrease in critical temperature; path-integral simulations for low
density (i.e. weak interactions) showed an increase in critical temperature.
The emerging consensus is that
$$
	\Delta T_c(a) = \frac{T_c{(a)}-T_c{(0)}}{T_c{(0)}}
$$
is linear in $a$ for small $a$.

Ueltschi makes the following summary of the theoretical work on this question.
(See also \cite{bib:arnoldmoore,bib:kps,bib:nl04} for a review on the
widely varying analytical and simulational results on $T_c{(a)}$; see
\cite{bib:baymetal} for a thorough listing of the progress up to 2001 and
\cite{bib:su09} for an up-to-date survey.)

\begin{itemize}

\item 1964:
\emph{Huang}:
$\Delta T_c(a) \sim (a \rho^{1/3})^{3/2}$, increases

\item 1971:
\emph{Fetter \& Walecka}: % Analysis
$\Delta T_c(a)$ decreases

\item 1982:
\emph{Toyoda}: % Analysis
$\Delta T_c(a)$ decreases

\item 1992:
\emph{Stoof}: % Analysis
$
\Delta T_c(a) = c \, a \rho^{1/3} + o(a \rho^{1/3}), \quad c>0
$

\item 1996:
\emph{Bijlsma \& Stoof}: % Analysis
$c=4.66$

\item 1997:
\emph{Gr\"uter, Ceperley, Lalo\"e}: % PIMC
$c=0.34$

\item 1999:
\emph{Holzmann, Gr\"uter, Lalo\"e}: % Analysis
$c=0.7$;
{\emph{ Holzmann, Krauth}}: $c=2.3$;

\item 1999:
\emph{Baym et.\ al.}: % Analysis
$c=2.9$

\item 2000:
\emph{Reppy et.\ al.}: % Analysis
$c=5.1$

\item 2001:
\emph{Kashurnikov, Prokof'ev, Svistunov}: % PIMC
$c=1.29$

\item 2001:
\emph{Arnold, Moore}: % phi^4 field theory
$c=1.32$

\item 2004:
\emph{Kastening}: % Analysis
$c=1.27$

\item 2004:
\emph{Nho, Landau}: % PIMC
$c=1.32$

\end{itemize}

%% ----------------------------------------------------------------
\section{Context of dissertation}
\label{sec:context}

The work of Ueltschi and Betz \cite{bib:bu07,bib:ducycles,bib:qmath} extends
the permutation point of view originated by Feynman, Penrose, and Onsager,
drawing on the work of S\"ut\H{o}, Buffet, and Pul\`e
\cite{bib:feynman,bib:penons,bib:suto1,bib:suto2,bib:buffetpule} to develop a
model of random spatial permutations.  (See chapter \ref{chap:model_defns} for
precise definition of the model.) In the permutation representation, this
condensate transition manifests itself as the onset of long permutation cycles.
The central point of this approach is that the system energy has been recast in
terms of permutations, which are amenable to analysis and simulation.  This
permits a new perspective on an old question:
%how does the critical temperature of Bose-Einstein condensation depend on
%inter-particle interaction strength?
the main goal of Ueltschi and Betz's long-term project is to quantify the shift
in critical temperature, as a function of interaction strength, for non-ideal
Bose gases in the small-scattering-length regime.

The interaction terms for the permutation representation of the Bose gas are
difficult to compute.  Moreover, it is interesting to consider the model of
random spatial permutations (which we sometimes refer to as the random-cycle
model) for its own sake.  Thus, Ueltschi, Betz, Gandolfo, Ruiz, and the author
take various approaches with varying degrees of fidelity to the physical
Bose-gas model.  In section \ref{sec:prb_model}, we will see that the
random-cycle model for $N$ particles is parameterized by $N$ \emph{cycle
weights} $\{\alpha_\ell\} = \alpha_1$, $\ldots$, $\alpha_N$.  The remainder of
this section involves results that may be obtained, analytically or
simulationally, when various constraints are placed on the cycle weights.

In the papers \cite{bib:qmath,bib:bu07}, Betz and Ueltschi examine the Bose-gas
permutation weights with point positions allowed to vary on the continuum; an
exact expression for the critical temperature is stated and proved for a
simplified interaction model in which only two-cycles interact.  The
cycle-weight parameter $\alpha_2$ is expressed in terms of the bosonic
hard-core scattering length $a$; all other cycle weights are set to zero.  In
\cite{bib:bu08}, this approach is extended to a model in which all the cycle
weights $\alpha_\ell$ may vary, but with the hypothesis that $\alpha_\ell$ goes
to zero faster than $1/\log(\ell)$.  Here, an expression for the shift in
critical temperature is found, as a function of all $N$ cycle weights.  It is
key to note that these $\alpha_\ell$'s are not computed directly from the
physical scattering length $a$; rather, the result obtained is true for
\emph{any} cycle weights $\{\alpha_\ell\}$ satisfying the small-cycle-weight
hypothesis.  In \cite{bib:bu09a}, Ueltschi and Betz estimate, to first order,
cycle weights $\{\alpha_\ell\}$ for the Bose gas.

Betz, Ueltschi, and Velenik \cite{bib:buv09} examine cycle weights
$\{\alpha_\ell\}$ with various hypotheses, including the Ewens \cite{bib:ewens}
case in which cycle weights are constant for all cycle lengths $\ell$.  These
random permutations are non-spatial, i.e. $T=0$ in the vocabulary of chapter
\ref{chap:model_defns}.

As is often the case in statistical mechanics, the study of this interacting
system necessitates the use of computational methods --- specifically,
Markov-chain Monte Carlo.  In \cite{bib:gru}, a simulational approach is taken
for points held fixed on the cubic unit lattice in the non-interacting case
($\alpha_\ell \equiv 0$ in the language of chapter \ref{chap:model_defns}).

This dissertation, the only known simulation approach to the interacting model,
applies MCMC methods to the case where point positions are held fixed on the
cubic unit lattice, with Ewens cycle weights.  This extends from\cite{bib:gru}
as well as \cite{bib:buv09}.  We find that even though lattice positions are
used, and even though the small-cycle-weight hypothesis is invalidated, one
nonetheless recovers the shift in critical temperature as predicted in the
small-cycle-weight model of \cite{bib:bu08}.

%% ----------------------------------------------------------------
\section{Literature review}

In addition to the many references made in previous sections of this chapter,
we point out the following.

The papers
\cite{bib:gcl97},
\cite{bib:kps},
\cite{bib:ceperley}, and
\cite{bib:nl04}
are among path-integral Monte Carlo simulational approaches to Bose-Einstein
condensation --- perhaps the closest relatives to the numerical work done in
\cite{bib:gru} and in this dissertation.

The doctoral dissertations of Peter Gr\"uter and Markus Holzmann
\cite{bib:gruter1996,bib:holzmann2000} are paradigmatic examples of clear
dissertation writing; the latter also provided insight into finite-size scaling.
%\cite{bib:feynman}
% Feynman,~R.P.
% \newblock {\em Atomic Theory of the $\lambda$ Transition in Helium}.
% \newblock The Physical Review, vol. 91, no. 6 (1953).
%
% \cite{bib:penons}
% Penrose,~O. and Onsager,~L.
% \newblock {\em Bose-Einstein Condensation and Liquid Helium}.
% \newblock The Physical Review, vol. 104, no. 3 (1956).
%
% \cite{bib:london}
% London,~F.
% \newblock Phys. Rev. vol. 54, no. 947 (1938).
%
%\cite{bib:boulder}
% Anderson,~M.H., Ensher,~J.R., Matthews,~M.R., Wieman,~C.E., and Cornell,~E.A.
% \newblock {\em Observation of Bose-Einstein Condensation in a Dilute Atomic
%   Vapor}.
% \newblock \emph{Science}, vol. 269, issue 5221 (Jul. 14 1995), 198-201.
%
% \cite{bib:baymetal}
% Baym,~G., Blaizot,~J.-P., Holzmann,~M., Lalo\"e,~F., and Vautherin,~D.
% \newblock {\em Bose-Einstein transition in a dilute interacting gas}.
% \newblock \texttt{arXiv:cond-mat/0107129v2}.
%
% \cite{bib:lssy}
% Lieb,~E.H., Seiringer,~R., Solovej,~J.P., and Yngvason,~J.
% \newblock {\em The Mathematics of the Bose Gas and its Condensation}.
% \newblock Oberwolfach Seminars, vol. 34.
% \newblock Birkh\"auser, 2005.
%
% \cite{bib:su09},
% Seiringer,~R. and Ueltschi,~D.
% \newblock {\em Rigorous upper bound on the critical temperature of dilute Bose
% gases}.
% %\newblock \texttt{arXiv.org:0904.0050}.
% \newblock Phys. Rev. B 80, 014502 (2009).
%
% \cite{bib:arnoldmoore}.
% Arnold,~X. and Moore,~X.
% \newblock {\em BEC transition temperature of a dilute homogeneous imperfect
%   Bose gas}.
% \newblock Phys. Rev. Lett. vol. 87, no. 12, Sept. 2001.

%% ----------------------------------------------------------------
%The model of random spatial permutations is developed by Ueltschi and Betz
%in
%\cite{bib:ducycles},
% Ueltschi,~D.
% \newblock {\em Feynman cycles in the Bose gas.}
% \newblock \texttt{arXiv:math-ph/0605002v2}.
%\cite{bib:bu07},
% Betz,~V. and Ueltschi,~D.
% \newblock{\em Spatial random permutations and infinite cycles}.
% %\newblock \texttt{arXiv:0711.1188}.
% \newblock Commun. Math. Phys. 285, 469-501 (2009).
%
%\cite{bib:qmath}.
% Ueltschi,~D.
% \newblock {\em The model of interacting spatial permutations and its relation
% to the Bose gas.}
% %\newblock \texttt{arXiv:0712.2443v3}.
% \newblock Mathematical Results in Quantum Mechanics, pp. 225-272,
% World Scientific (2008).
%and \cite{bib:bu08}.
% Betz,~V. and Ueltschi,~D.
% \newblock {\em Spatial random permutations with small cycle weights.}
% %\newblock \texttt{arXiv:0812.0569v1}.
%Their work uses the Fourier method of Buffet and Pul\'e \cite{bib:buffetpule}.
% \cite{bib:buffetpule}
% Buffet,~E. and Pul\`e,~J.V.
% \newblock {\em Fluctuation properties of the imperfect Bose gas}.
% \newblock J. Math. Phys. vol. 24, no. 6, June 1983.
%
% \cite{bib:suto1} and \cite{bib:suto2}
% S\"ut\H{o},~A.
% \newblock {\em Percolation transition in the Bose gas.}
% \newblock J. Phys. A: Math. Gen. \textbf{26} (1993) 4689-4710.
%
% S\"ut\H{o},~A.
% \newblock {\em Percolation transition in the Bose gas II.}
% \newblock J. Phys. A: Math. Gen. \textbf{35} (2002) 6995-7002.
% \newblock Probabl. Th. Rel. Fields (2010).
%The only published work on simulation of the random-cycle model, which treats
%the non-interacting model, is \cite{bib:gru}.
% Gandolfo,~D., Ruiz,~J., and Ueltschi,~D.
% \newblock {\em On a model of random cycles}.
% %\newblock \texttt{arXiv:cond-mat/0703315}.
% \newblock Statist. Phys. 129, 663-676 (2007).
%This dissertation is the first published work on simulations of the
%interacting model.

% %% ----------------------------------------------------------------
Mean longest cycle for uniformly distributed (i.e. non-interacting,
non-spatial) permutations ($T=0$ and $\alpha=0$ in the language of chapter
\ref{chap:model_defns}) was discussed by \cite{bib:shepplloyd}, following a
question posed by Golomb on the basis of experimental data \cite{bib:golomb}.
See also section \ref{sec:conj}.
% Shepp,~L.A. and Lloyd,~S.P.
% \newblock{\em Ordered Cycle Length in a Random Permutation}.
% \newblock Trans. Amer. Math. Soc. 121, (1966), 340-357.
%
% Golomb,~S.W.
% \newblock {\em Random permutations}.
% \newblock Bull. Amer. Math. Soc. 70 (1964), 747.
%
Non-uniformly distributed non-spatial permutations with constant cycle weights
($T=0$ and $\alpha\ne 0$ in the language  of chapter \ref{chap:model_defns})
arose in mathematical biology \cite{bib:ewens}.
% Ewens,~W.J.
% \newblock {\em The sampling theory of selectively neutral alleles}.
% \newblock Theor. Popul. Biol. 3, 87-112 (1972).

% %% ----------------------------------------------------------------
Quantum mechanics and statistical mechanics, as background for the work
of Ueltschi and Betz, may be found in
%
\cite{bib:griffiths},
% Griffiths,~D.J.
% \newblock {\em Introduction to Quantum Mechanics} (2nd ed.).
% \newblock Pearson Prentice Hall, 2005.
%
\cite{bib:huang},
% Huang,~K.
% \newblock {\em Introduction to Statistical Physics}.
% \newblock CRC Press, 2001.
%
and \cite{bib:sakurai}.
% Sakurai,~J.J.
% \newblock {\em Modern Quantum Mechanics} (2nd ed.).
% \newblock Addison Wesley, 1993.

% %% ----------------------------------------------------------------
Worm algorithms for path-integral Monte Carlo, which inspired my
random-spatial-permutation worm algorithm of chapter \ref{chap:worm}, are used
throughout simulational physics.  See in particular
%
\cite{bib:bps06} and
% Boninsegni,~M., Prokof'ev,~N.V., and Svistunov,~B.V.
% \newblock {\em Worm algorithm and diagrammatic Monte Carlo:  A new approach
%   to continuous-space path integral Monte Carlo simulations.}
% \newblock Physical Review E \textbf{74}, 036701 (2006).
%
\cite{bib:pst98}.
% Prokof'ev,~N.V., Svistunov,~B.V., and Tupitsyn,~I.S.
% \newblock {\em Exact, complete, and universal continuous-time worldline Monte
%   Carlo approach to the statistics of discrete quantum systems.}
% \newblock Journal of Experimental and Theoretical Physics, vol. 87, no. 2
%   (1998).

% %% ----------------------------------------------------------------
Finite-size scaling techniques are employed for path-integral Monte Carlo
simulations in
%
\cite{bib:gcl97},
% Gr\"uter,~P. and Ceperley,~D.
% \newblock {\em Critical temperature of Bose-Einstein condensation of
%   hard-sphere gases.}
% \newblock Physical Review Letters, vol. 79, no. 19 (1997).
%
\cite{bib:hk99},
% Holzmann,M. and Krauth,~W.
% \newblock {\em Transition temperature of the homogeneous, weakly interacting
%   Bose gas.}
% \newblock Physical Review Letters, vol. 83, no. 14 (1999).
%
\cite{bib:kps},
% Kashurnikov,~V.A., Prokof'ev,~N.V., and Svistunov,~B.V.
% \newblock {\em Critical temperature shift in weakly interacting Bose gas}.
% \newblock Physical Review Letters, vol. 87, no. 12 (2001).
%
\cite{bib:nl04},
% Nho,~K. and Landau,~D.P.
% \newblock {\em Bose-Einstein condensation temperature of a homogeneous
%   weakly interacting Bose gas:  Path integral Monte Carlo study.}
% \newblock Physical Review A \textbf{70}, 053614 (2004)
%
\cite{bib:pc87},
% Pollock,~E.L. and Ceperley,~D.M.
% \newblock {\em Path-integral computation of superfluid densities.}
% \newblock Physical Review B, vol. 36, no. 16 (1987).
%
\cite{bib:pgp08},
% Pilati,~S., Giorgini,~S., and Prokof'ev,~N.V.
% \newblock {\em Critical temperature of interacting Bose gases in two and
%   three dimensions.}
% \newblock Physical Review Letters \textbf{100}, 140405 (2008).
%
and \cite{bib:pr92}.
% Pollock,~E.L. and Runge,~K.J.
% \newblock {\em Finite-size-scaling analysis of a simulation of the
%   ${}^4$He superfluid transition}.
% \newblock Physical Review B, vol. 46, no. 6 (1992).
%
Citation trails in the above-cited works lead back to \cite{bib:barber}.
%
Some background information is found in \cite{bib:landaubinder};
% Landau,~D.P. and Binder,~K.
% \newblock {\em A Guide to Monte Carlo Simulations in Statistical Physics}
%   (2nd ed.).
% \newblock Cambridge University Press, 2005.
%
an excellent survey, encompassing and explicating all the above methods ---
truly a blessing for the aspiring learner --- is \cite{bib:cggp}.
% Caracciolo,~S., Gambassi,~A., Gubinelli,~M., and Pelisetto,~A.
% \newblock {\em Finite-Size Scaling in the Driven Lattice Gas.}
% \newblock Journal of Statistical Physics, vol. 115, Nos. 1/2, April 2004.
% \newblock \texttt{arxiv:cond-mat/0312175}.

% %% ----------------------------------------------------------------
Markov chain Monte Carlo methods are discussed in \cite{bib:landaubinder};
% Landau,~D.P. and Binder,~K.
% \newblock {\em A Guide to Monte Carlo Simulations in Statistical Physics}
%   (2nd ed.).
% \newblock Cambridge University Press, 2005.
%
this dissertation has been influenced most heavily by \cite{bib:berg}.  Indeed,
my appendix \ref{appx:eta} is an elaboration on Berg's discussion of integrated
autocorrelation time.
% Berg,~B.
% \newblock {\em Markov Chain Monte Carlo Simulations and Their Statistical
%   Analysis}.
% \newblock World Scientific Publishing, 2004.
%
The probability background necessary for either Landau and Binder or Berg
may be found, with increasing levels of sophistication, in
\cite{bib:lawler}, \cite{bib:gs}, and \cite{bib:oksendal}.
% Grimmett,~G. and Stirzaker,~D.
% \newblock {\em Probability and Random Processes}, 3rd ed.
% \newblock Oxford, 2001.
%
% Lawler,~G.
% \newblock {\em Introduction to Stochastic Processes} (2nd ed.).
% \newblock Chapman \& Hall/CRC, 2006.
%
% {\O}ksendal,~B.
% \newblock {\em Stochastic Differential Equations} (6th ed.).
% \newblock Springer, 2007.
%
The standard reference for statistical analysis, including confidence
intervals, is \cite{bib:cb}.
% Casella,~G. and Berger,~R.L.
% \newblock {\em Statistical Inference} (2nd ed.).
% \newblock Duxbury Press, 2001.
More practical aspects of the statistical reduction of experimental data
are found in \cite{bib:young}.
% Young,~H.D.
% \newblock {\em Statistical Treatment of Experimental Data}.
% \newblock McGraw-Hill, 1962.

% %% ----------------------------------------------------------------
The Mersenne Twister \cite{bib:mersenne} is the pseudo-random-number generator
used in the this dissertation's computational work.  Another good generator is
pseudo-DES \cite{bib:nr}.  Moreover, any numerical dissertation without a
reference to \emph{Numerical Recipes} is incomplete; its inclusion here is as
good a point as any to end the literature review.
% Matsumoto,~M. and Nishimura,~T.
% \newblock{Mersenne twister: a 623-dimensionally equidistributed uniform
% pseudo-random number generator}.
% \newblock ACM Transactions on Modeling and Computer Simulation,
% Volume 8, Issue 1 (Jan. 1998), pp 3-30.
%
% Press,~W. et al.
% \newblock {\em Numerical Recipes} (2nd ed.).
% \newblock Cambridge, 1992.

Last, we delineate the originality of work presented in this dissertation.
Chapters \ref{chap:history}-\ref{chap:RV_theory} of this dissertation are a
rephrasing and an elaboration on \cite{bib:bu07,bib:bu08,bib:qmath}.  Chapter
\ref{chap:MCMC} is quite standard; the contribution made here is to present
familiar general results in the specific context of random permutations.  The
essential SO algorithm of chapter \ref{chap:SO_SAR}, with a small modification,
was presented in \cite{bib:gru}; likewise for the SO $\Delta H$ computations in
chapter \ref{chap:Delta_H}.  The treatment here is the first correctness proof
of the SO algorithm.  The SAR algorithm of chapter \ref{chap:SO_SAR} is due to
Daniel Ueltschi.  The band-update algorithm (chapter \ref{chap:band_updates})
is due to the author.  The worm algorithm and its correctness proof (chapter
\ref{chap:worm}) are due to the author, along with the remaining $\Delta H$
computations of chapter \ref{chap:Delta_H}.  The remaining chapters,
\ref{chap:sw}-\ref{chap:future_work}, are also original work.  Appendix
\ref{appx:eta} is a new take on an old question; see also \cite{bib:berg}.  The
correlated-uniform Markov process is original, as is the explicit comparison of
batched and non-batched means for exponentially correlated stationary Markov
processes.

%% ================================================================
\chapter{The model of random spatial permutations}
\label{chap:model_defns}

Here we review concepts from \cite{bib:bu07,bib:bu08}, fixing notation and
intuition to be used in the rest of the paper.

%% ----------------------------------------------------------------
\section{The probability model}
\label{sec:prb_model}

The state space is $$ \Omega_{\Lambda,N} = \Lambda^N \times \SN $$ where
$\Lambda=[0,L]^3$ with periodic boundary conditions and $\SN$ is the group of
\plainidx{permutations} of $N$ points\footnote{One may of course consider
$\Lambda=[0,L]^d$ for $d=1,2$, but for this paper, $d=3$ only.}. Point
positions are $\vecX = (\vecx_1, \ldots, \vecx_N)$ for $\vecx_1, \ldots,
\vecx_N \in \Lambda$.  These are called spatial permutations in that they
involve the permutation $\pi$ as well as the $N$ point positions $\vecx_i$.
See figure \ref{fig:pmt1}.

The \plainidx{probability measure} on this state space
is constructed via Gibbs measure on a \index{$H$}\plainidx{Hamiltonian}.  The
background probability measure is discrete (uniform) in $\pi$ and continuous
(Lebesgue) in $\vecX$.  The Hamiltonian takes
one of two forms.  In the first, relevant to the Bose gas, we have
\begin{align}
\label{eqn:RCM_H_gen_V}
	H(\vecX,\pi) &= \frac{T}{4} \sum_{i=1}^N \|\vecx_i-\vecxpii\|^2
		+ \sum_{1 \le i<j \le N} V(\vecx_{i}, \vecx_{\pi(i)}, \vecx_{j}, \vecx_{\pi(j)})
\end{align}
where $T=1/\beta$ and the $V$ terms are interactions between permutation jumps.
The operator $H$ is unbounded\index{unbounded operator}, but it is
symmetric\index{symmetric operator} so we consider its \plainidx{self-adjoint
extension}.  We take its domain to be $f$ in $C^2(\Lambda^N)$ with periodic
boundary conditions.  The permutation jump $\vecxi\mapsto\vecxpii$ interacts
with the permutation jump $\vecxj\mapsto\vecxpij$.  The temperature scale
factor $T/4$, not $\beta/4$, is atypical in statistical mechanics.  For
purposes of the current work, this may be considered an \plainidx{ansatz}; in
\cite{bib:bu07}, this choice of scale factor is be shown to be appropriate for
permutation representation of the Bose gas.  In the second form of the
Hamiltonian, considered in this paper, we use interactions which are dependent
solely on cycle length:
\begin{align}
\label{eqn:RCM_H_rell_V}
	H(\vecX,\pi) &= \frac{T}{4} \sum_{i=1}^N \|\vecx_i-\vecxpii\|^2
		+ \sum_{\ell=1}^N \alpha_\ell r_\ell(\pi),
\end{align}
where $r_\ell(\pi)$ is the number of $\ell$-cycles in $\pi$ and the
$\alpha_\ell$'s are free parameters, called \emph{cycle weights}.  One
ultimately hopes to choose the $\alpha_\ell$'s appropriately for the Bose gas;
even if not, the model is well-defined and of its own mathematical interest.

\begin{figure}
\begin{center}
\psfragscanon
\includegraphics[scale=1.0]{figures/pmt1.eps}
\caption[A spatial permutation on $N=26$ points.]
	{A spatial permutation on $N=26$ points.  There are 11 one-cycles, three
	two-cycles, one four-cycle, and one five-cycle.  We say $r_1(\pi)=11$,
	$r_2(\pi)=3$, $r_4(\pi)=1$, $r_5(\pi)=1$, and $r_\ell(\pi) = 0$ for all
	other $\ell$.
	\label{fig:pmt1}}
\end{center}
\end{figure}

The first contribution to the energy\footnote{The papers \cite{bib:bu07} and
\cite{bib:qmath} generalize from $\|\vecx_i-\vecxpii\|$ to
$\xi(\vecx_i,\vecxpii)$ where $\xi$ is a spherically symmetric
non-negative-valued function on $\R^d$ having integrable $e^{-\xi}$.  This
generalization is not of interest in this dissertation.} is the sum of squares
of permutation jump lengths.  This makes permutations with long jumps
disfavored; permutations with many short jumps will be less strongly
disfavored.  The second contribution to the energy involves cycle weights.  We
consider only small cycle weights, which perturb the critical temperature but
which do not qualitatively modify the effects of the distance-related terms.
More intuition for the model will be presented in section
\ref{sec:qual_long_cycles}.

We immediately restrict our attention to interactions which are
\plainidx{conjugacy class functions}\index{class functions}, i.e. which are
constant on permutation cycles.  That is, for such a random variable $X$, we
have $X(\sigma^{-1}\pi\sigma) = X(\pi)$ for all $\sigma,\pi\in\SN$.  Also, we
have $X(\pi) = X(r_1(\pi), \ldots, r_N(\pi))$ where $r_\ell(\pi)$ counts the
number of $\ell$-cycles in $\pi$ for $\ell$ from $1$ to $N$.

Choices of point positions $\vecX$ yield two cases: (1) In the
\emphidx{annealed model}, point positions are variable on the
\plainidx{continuum} (here, uniform random) and are averaged over.  This model
is examined analytically in \cite{bib:bu07,bib:qmath,bib:bu08}.  (2) In the
\emphidx{quenched model}, point positions are held fixed (here, on the cubic
unit \plainidx{lattice}).  This model is examined simulationally in
\cite{bib:gru}, and in this dissertation.  We often write $H(\pi)$ in place of
$H(\vecX,\pi)$ since we either work on a lattice where the $\vecx_i$'s are held
fixed, or on the continuum where the $\vecx_i$'s are integrated out.  Thus, the
system energy $H$ (as well as all other random variables we consider) is a
function of the permutation $\pi$.

%\textbf{Partition functions:}
We consider two \index{partition function}partition functions, for a fixed
point configuration $\vecX$ (quenched model) and for an average over point
configurations (annealed model), respectively:
$$
	Y(\Lambda,\vecX) = \sum_{\sigma\in\SN} e^{-H(\vecX,\sigma)}
	\qand
	Z(\Lambda,N) = \frac{1}{N!} \int_{\Lambda^N} Y(\Lambda,\vecX)\,d\vecX.
$$
%\textbf{Probability measures:}
Fixing (quenched) point positions $\vecX$, we have a \plainidx{discrete
distribution} on $\pi$:
\begin{align}
\label{eqn:RCM_Gibbs_P}
	Y(\Lambda,\vecX) &= {\sum_{\sigma\in\SN} e^{-H(\vecX,\sigma)}}, &
	P(\pi) = P_{\Lambda,\vecX}(\pi)
	&= \frac{e^{-H(\vecX,\pi)}} {Y(\Lambda,\vecX)}.
\end{align}
For varying point positions (e.g. for considerations of the Bose gas), we have
a \plainidx{joint distribution} which is continuous in $\vecX$ and discrete in
$\pi$:
\begin{align*}
	P_{\Lambda,N}(\vecX,\pi)\,d\vecX &= \frac{e^{-H(\vecX,\pi)} d\vecX}
		{Z(\Lambda,N)}.
\end{align*}
From this we obtain two \plainidx{marginal distributions}.  If we integrate
over point configurations $\vecX$, then we obtain a discrete distribution on
$\SN$:
\begin{align*}
	P_{\Lambda,N}(\pi) &=
		\frac{1}{N!} \int_{\Lambda^N} d\vecX\, P_{\Lambda,N}(\vecX,\pi)
	= \frac{
			\frac{1}{N!} \int_{\Lambda^N} d\vecX\, e^{-H(\vecX,\pi)}
		}{
			\frac{1}{N!} \int_{\Lambda^N} d\vecX\,
				\sum_{\sigma\in\SN} e^{-H(\vecX,\sigma)}
		}
	= \frac{
			\int_{\Lambda^N} d\vecX\, e^{-H(\vecX,\pi)}
		}{
			Z(\Lambda,N)\,N!
		}.
\end{align*}
If, on the other hand, we sum over permutations, then we obtain a
\plainidx{continuous distribution} for point configurations:
\begin{align}
\label{eqn:true_dist}
	P_{\Lambda,N}(\vecX)\,d\vecX
	&= \sum_{\pi\in\SN} P_{\Lambda,N}(\vecX,\pi)\,d\vecX
	= \frac{Y(\Lambda,\vecX)\,d\vecX}{Z(\Lambda,N)}.
\end{align}
This continuous distribution is certainly of interest:  it is the point
distribution for the Bose gas, when the Hamiltonian is appropriately chosen.
However, it is very difficult to compute \cite{bib:lls}.  From here on, we
consider the two discrete distributions on $\SN$, namely,
\begin{align}
\label{eqn:RCM_P}
	P_{\Lambda,\vecX}(\pi)
	&= \frac{e^{-H(\vecX,\pi)}}{Y(\Lambda,\vecX)}
	\qand
	P_{\Lambda,N}(\pi)
	= \frac{\int_{\Lambda^N}d\vecX e^{-H(\vecX,\pi)}}{Z(\Lambda,N)N!}.
\end{align}

%\textbf{Expectations:}
For\index{expectation} a random variable $X(\pi)$, we have
$$
	\bbE_{\Lambda,\vecX}[X] =
	\frac{
		\sum_{\pi\in\SN} X(\pi) e^{-H(\vecX,\pi)}
	}{Y(\Lambda,\vecX)}
	\qand
	\bbE_{\Lambda,N}[X] =
	\frac{
		\int_{\Lambda^N} d\vecX \sum_{\pi\in\SN} X(\pi) e^{-H(\vecX,\pi)}
	}{Z(\Lambda,N) N!}.
$$
In either case, we also write the probability $P(\pi)$ and the expectation as
\begin{align}
\label{eqn:E_pi}
	\bbE[X] &= \sum_{\pi\in\SN} P(\pi) X(\pi).
\end{align}

%% ----------------------------------------------------------------
\section{Model variants by choice of cycle weights}

The model of random spatial permutations is in fact a family of models.  As
described in the previous section, point positions may be annealed or quenched,
the latter being the case for this dissertation.  Likewise, various constraints
may be places on the cycle weights.  Recall from equation
\ref{eqn:RCM_H_rell_V} that
\begin{align}
	H(\vecX,\pi) &=
	\frac{1}{4\beta} \sum_{i=1}^N \|\vecx_i-\vecxpii\|^2
	+ \sum_{\ell=1}^N \alpha_\ell r_\ell(\pi).
\end{align}
There are $N$ free parameters $\alpha_\ell$, and thus many models of spatial
permutations.

If $\alpha_\ell=0$ for all $\ell$, one obtains the \emphidx{non-interacting
case}.
%
When $\alpha_2=\alpha$ and $\alpha_\ell=0$ for $\ell \ne 2$, we have the
\emphidx{two-cycle model} \cite{bib:bu07,bib:qmath}.  As $\alpha$ is increased,
two-cycles are discouraged.
%
Otherwise, we have the \emphidx{general-cycle model}.  This splits into (at
least) three submodels:
%
Betz and Ueltschi, in \cite{bib:bu08}, consider the case where $\alpha_\ell$
tend toward zero faster than $1/\log(\ell)$.  (Note that this includes the
$r_2$ model as a special case.)
%
Ideally, one would have $\alpha_\ell$'s which match the Brownian-bridge
interactions for the Bose gas.  Work in this direction has recently been done
by Betz and Ueltschi \cite{bib:bu09a}, but is beyond the scope of this
dissertation.
%
For this dissertation, we consider the case where $\alpha_\ell=\alpha$ is
constant in $\ell$.  This is the \emphidx{Ewens distribution} \cite{bib:ewens}.
In summary, the model considered in this dissertation uses point positions
quenched on the cubic unit lattice, with small non-negative  Ewens cycle
weights.

%% ----------------------------------------------------------------
\section{Qualitative characterization of long cycles}
\label{sec:qual_long_cycles}

Having seen definitions for the probability model, one next asks what a typical
random spatial permutation looks like.  In this section we develop intuition;
in section \ref{sec:ordps}, we construct quantitative descriptions of the ideas
presented here.

\begin{figure}[!htb]
\psfragscanon
\begin{center}
\includegraphics[scale=0.6]{figures/dis_dot_plot_high_T.eps}
\qquad
\includegraphics[scale=0.6]{figures/dis_dot_plot_medium_T.eps}
\qquad
\includegraphics[scale=0.6]{figures/dis_dot_plot_low_T.eps}
\caption
	[Some typical permutations for varying high $T$.]
	{Points and permutation jumps, for a typical permutation at
	high $T$ (short cycles of short jumps), medium but subcritical $T$ (all
	jump lengths are short, with occasional long cycles thereof), and low $T$
	(jump lengths are arbitrary).
	\label{fig:dot_plots}}
\end{center}
\end{figure}

As $T \to \infty$, the probability measure becomes supported only on the
identity permutation:  the distance-dependent terms in equation
\ref{eqn:RCM_H_rell_V} are large whenever any jump has non-zero length.  For
large but finite $T$, the length-dependent terms penalize permutation jumps
from a site to any site other than itself.  Thus, for large $T$ we expect the
identity permutation to be the most likely, with occasional 2-cycles, 3-cycles,
etc. involving nearby points.  On the other extreme, as $T \to 0$,
length-dependent terms go to zero and the probability measure approaches the
uniform distribution on $\SN$:  the distance-dependent terms all go to zero.
For intermediate $T$, one observes that the length $\|\pi(\vecx)-\vecx\|$ of
each permutation jump remains small, increasing smoothly as $T$ drops.  See
figure \ref{fig:mean_max_jump_len_LTa} for more precise information.

The intermediate temperature regime is the one of interest.  As is found in
theoretical and simulational work, as detailed through the rest of this
dissertation, this regime has the following properties.  For $T$ below a
critical temperature $T_c$, there is a \plainidx{phase transition}: while
individual jump lengths remain short, \emph{arbitrarily long cycles form}.  See
figure \ref{fig:dot_plots} for depictions of typical permutations at high $T$,
subcritical $T$, and low $T$.  In the non-interacting case, $T_c$ is
approximately 6.86; interactions in the form of positive $\alpha$ terms
interactions increase $T_c$.  Quantifying that dependence, i.e.
\begin{align}
\label{eqn:Delta_Tc}
	\Delta T_c(\alpha) &= \frac{T_c(\alpha)-T_c(0)}{T_c(0)},
\end{align}
as a function of $\alpha$ for small positive $\alpha$, is the central goal of
this dissertation.

From figures such as \ref{fig:dot_plots}, one can see long cycles visually.
How do we measure them numerically?  Let $\ellmax(\pi)$ be the length of the
longest cycle in $\pi$, with $\bbE[\ellmax]$ its mean over all permutations.
We observe that for $T$ above $T_c$, $\bbE[\ellmax]$ is constant as
$L\to\infty$: the mean length of the longest cycle does not increase as $L$
increases.  That is to say, $\bbE[\ellmax]/N$ goes to zero as $L\to\infty$.
For $T$ below $T_c$, the length of the longest cycle does increase as $L$
increases.  We find that $\bbE[\ellmax]$ scales with $N$, i.e.
$\bbE[\ellmax]/N$ approaches a temperature-dependent constant as $L\to\infty$:
there are arbitrarily long cycles, or infinite cycles, in the infinite-volume
limit.  See figure \ref{fig:E_ellmax_N} for plots of $\bbE[\ellmax]/N$ as a
function of $T$ for various system sizes with $N=L^3$.  See also figure
\ref{fig:r_small_and_lmax} of section \ref{sec:comput_lmax_GRU_quot}.  Precise
information about $\bbE[\ellmax]/N$ and other quantities is found in section
\ref{sec:ordps}.

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.6]{../tests/fss3a/resultsfigs/graph_rell_alpha0_0.000_E_ellmax_N.eps}
\includegraphics[scale=0.6]{../tests/fss3a/resultsfigs/graph_rell_alpha0_0.001_E_ellmax_N.eps}
\caption
	[Order parameter $f_M = \bbE\ellmax/N$.]
	{Order parameter $f_M = \bbE[\ellmax]/N$ for finite systems, with
	$\alpha=0.000,0.001$.  Interactions increase the critical temperature.
	\label{fig:E_ellmax_N}}
\end{center}
\end{figure}

%% ----------------------------------------------------------------
\section{Known results and conjectures}
\label{sec:conj}

In this dissertation we study chiefly the $\alpha$-dependent shift in critical
temperature for Ewens cycle weights and cubic-lattice point positions.  We also
consider an $\alpha$-dependent macroscopic-cycle quotient.

Known results for point locations averaged over the continuum (as discussed in
section \ref{sec:prb_model}) are obtained largely using Fourier methods
\cite{bib:bu08} which are unavailable for point positions held fixed on the
lattice.  Betz and Ueltschi have determined $\Delta T_c(\alpha)$, to first
order in $\alpha$, for two-cycle interactions \cite{bib:bu07} and small cycle
weights \cite{bib:bu08}.  The critical $(\rho,T,\alpha)$ manifold relates
$\rho_c$ to $T_c$.  Specifically, Betz and Ueltschi obtain the following
\cite{bib:bu08} with the small-cycle-weight constraint on the cycle weights
$\{\alpha_\ell\}$:
\begin{align*}
	\rho_c(\alpha_1, \alpha_2, \ldots)
	&\approx \frac{1}{(4\pi\beta)^{3/2}}
		\sum_{\ell \ge 1} e^{-\alpha_\ell} \ell^{-3/2}.
\end{align*}
In the non-interacting special case, all cycle weights are zero, and we have
\begin{align*}
	\rho_c(0)
	&\approx \frac{1}{(4\pi\beta)^{3/2}}
		\sum_{\ell \ge 1} \ell^{-3/2}
	= \frac{\zeta(3/2)}{(4\pi\beta)^{3/2}}.
\end{align*}
For the two-cycle special case \cite{bib:bu07}, $\alpha_2$ is non-zero, and the
other cycle weights are zero.  We have
\begin{align*}
	\rho_c(\alpha)
	&= \frac{1}{(4\pi\beta)^{3/2}}
		\left(
		e^{-\alpha} 2^{-3/2} +
		\sum_{\ell \ne 2} e^{-\alpha_\ell} \ell^{-3/2} \right)
	= \rho_c(0) - \frac{(1-e^{-\alpha})}{(8\pi\beta)^{3/2}}.
\end{align*}
The Ewens case does not satisfy the small-cycle-weight constraint; all the
$\alpha_\ell$'s are the same.  Nonetheless, using Betz and Ueltschi's formula
we obtain
\begin{align*}
	\rho_c(\alpha)
	&\approx \frac{1}{(4\pi\beta)^{3/2}}
		\sum_{\ell \ge 1} e^{-\alpha} \ell^{-3/2}
	= \frac{e^{-\alpha} \zeta(3/2)}{(4\pi\beta)^{3/2}}.
\end{align*}

These expressions give $\rho_c$ as a function of $\beta$.  For the cubic
lattice, we fix $\rho \equiv 1$ and thus obtain $\beta_c$, or equivalently
$T_c$:
\begin{align*}
	T_c(\alpha_1, \alpha_2, \ldots) &\approx
	\frac{4 \pi}
		{\left(\sum_{\ell \ge 1} e^{-\alpha_\ell}\ell^{-3/2}\right)^{2/3}}.
\end{align*}
For the non-interacting special case, this is
\begin{align}
\label{eqn:theo_Tc}
	T_c(0) &\approx \frac{4 \pi}{\zeta(3/2)^{2/3}} \approx 6.625.
\end{align}
For the two-cycle special case, one has
\begin{align*}
	T_c(\alpha) &\approx \frac{4 \pi}
		{\left[\zeta(3/2) - 2^{3/2}(1-e^{-\alpha}) \right]^{2/3}}.
\end{align*}
For the Ewens case, again we violate the small-cycle-weight hypothesis, but
we may nonetheless write down the expression
\begin{align*}
	T_c(\alpha) &\approx \frac{4 \pi}
		{\left[e^{-\alpha} \zeta(3/2)\right]^{2/3}}
		\approx 6.625 e^{2\alpha/3}.
\end{align*}
Taylor-expanding in the small parameter $\alpha$, the shift in critical
temperature is then
\begin{align}
\label{eqn:rho_c}
	\Delta T_c(\alpha)
	&= \frac{T_c(\alpha)-T_c(0)}{T_c(0)}
	\approx e^{2\alpha/3} - 1 \approx \frac{2\alpha}{3} &
	c &\approx 0.667.
\end{align}
(Note that this is not in conflict with the constant $c$ in section
\ref{sec:context}:  there, one examines $\Delta T_c(a)$ where $a$ is the
scattering length of the interacting Bose gas; here, one has $\Delta
T_c(\alpha)$ for free parameter $\alpha$.) We inquire whether this result,
obtained for small cycle weights with point positions varying on the continuum,
holds for Ewens weights with point positions held fixed on the lattice.  We
suspect that the fine details of point positions are unimportant for the shift
in critical temperature.  For Ewens interactions, though, $\Delta T_c(\alpha)$
is theoretically unknown for Ewens interactions with points either on the
continuum or on the lattice.  The simulational treatment in this paper is the
only known attack on this question.

In addition to our main result on the shift in critical temperature, we also
examine the fraction of sites in macroscopic cycles.  Namely, for $\alpha_\ell
\equiv 0$ (the non-interacting model), we observe empirically that the
\emph{macroscopic-cycle quotient} $\bbE[\ellmax] / N f_I$ is constant for $T$
below but near $T_c$.  (That is, the two order parameters $f_I$ and
$\bbE[\ellmax]/N$ have the same critical exponent.) For uniform-random
(non-spatial) permutations, Shepp and Lloyd 1966 \cite{bib:shepplloyd} solved
Golomb's 1964 question \cite{bib:golomb}: $\bbE[\ellmax]/N \approx 0.6243$.
Unpublished work of Betz and Ueltschi has found $\bbE[\ellmax]/N f_I$ holds
that same value for random spatial permutations in the non-interacting case
$\alpha_\ell \equiv 0$.  The intuition is that long cycles are in some sense
uniformly distributed within the zero Fourier mode.  We conjecture that
$\bbE[\ellmax] / N f_I$ is $\alpha$-dependent but constant in $T$ (for $T$
below but near $T_c$) for all interaction models.

%% ================================================================
\chapter{Random variables}
\label{chap:RV_theory}

Specific random variables under consideration include the following.  For each,
the relevant theory sections (in this chapter) and experiment sections (in
chapter \ref{chap:sw}) are indicated.

\begin{itemize}

\item System energy $H(\pi)$, with non-interacting terms $D(\pi)$ and
interacting terms $V(\pi)$.  Energy per site $h(\pi)=H(\pi)/N$,
$d(\pi)=D(\pi)/N$, $v(\pi)=V(\pi)/N$.  These are as discussed in section
\ref{sec:prb_model}; explicit computation is discussed in section
\ref{sec:comput_HDV}.

\item $r_\ell(\pi)$:  the number of $\ell$-cycles in the permutation for
$\ell=1, 2, \ldots, N$.  The definitions are familiar from elementary algebra;
computation is discussed in section \ref{sec:comput_rell}.

\item Cycle length $\ellxpi$, spatial cycle length $s_\vecx(\pi)$; correlation
lengths $\ellbar$ and $\xi:=\essbar$: sections \ref{sec:theory_cyclen_corrlen}
and \ref{sec:comput_cyclen_corrlen}.

\item Mean jump length and maximum jump length: sections \ref{sec:theory_mjl}
and \ref{sec:comput_mjl}.  These are used to confirm the hypothesis of short
jump lengths as mentioned in section \ref{sec:theory_wno}.

% \ref{sec:inf_cyc_existence}
\item $\rho_{mn}$ and $f_I$:  section \ref{sec:theory_fI} and (theory) and
\ref{sec:comput_fI} (experiment).  These are the density of sites in specified
lengths and the fraction of sites in infinite cycles, respectively.

\item Longest cycle length $\ellmax(\pi)$ and macroscopic-cycle \cite{bib:gru}
quotient $E[\ellmax] / N f_I$:  sections \ref{sec:theory_lmax_GRU_quot} and
\ref{sec:comput_lmax_GRU_quot}.

\item Winding numbers $W_x(\pi)$, $W_y(\pi)$, $W_z(\pi)$, scaling winding
number $f_S$, and fraction $f_W$ of sites in winding cycles:  sections
\ref{sec:theory_wno} and \ref{sec:comput_wno}.

\end{itemize}

Given a random variable $X$, let $Q = \bbE[X]$.  For each quantity $Q$, one
should distinguish between the finite-volume value $Q_L(T)$ and the
infinite-volume limit $Q_\infty(T) = \lim_{L\to\infty} Q_L(T)$.  For this
chapter, and omitted subscript indicates the infinite-volume value.  The
difference becomes clear in chapter \ref{chap:sw}; at that point, I will
carefully distinguish between $Q_L(T)$ and $Q_\infty(T)$.

%% ----------------------------------------------------------------
\section{Cycle lengths and correlation length $\xi$}
\label{sec:theory_cyclen_corrlen}

\begin{defn}
\label{def:ell}
Fixing\index{$\ell$}\index{$s$} $\pi\in\SN$ and $\vecx \in \Lambda$, the
\plainidx{cycle length} $\ellxpi$ is the smallest positive number $a$ such that
$\pi^a(\vecx) = \vecx$.  This is nothing more than the length of the cycle
containing $\vecx$: if $\vecx$ is in a 3-cycle, then $\ellxpi = 3$.  We define
the \plainidx{spatial cycle length}
$$
	\essx(\pi) = \sum_{j=1}^{\ellxpi} \|\pi^j(\vecx)-\pi^{j-1}(\vecx)\|.
$$
This is simply the sum of Euclidean jump distances for all permutation jumps in
the cycle containing $\vecx$.  We may at times instead write, respectively,
$$
	\ell_i(\pi) = \ell_{\vecx_i}(\pi)
	\qqor
	s_i(\pi) = s_{\vecx_i}(\pi).
$$
\end{defn}

For example, for the point configuration $\vecX$ and the permutation $\pi$
in figure \ref{fig:pi8}, we have
$$
	\ell_1(\pi) =
	\ell_2(\pi) =
	\ell_3(\pi) =
	\ell_4(\pi) = 4,
	\quad
	\ell_5(\pi) =
	\ell_6(\pi) =
	\ell_7(\pi) = 3,
	\qand
	\ell_8(\pi) = 1.
$$

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics{figures/pi8.eps}
\caption{A configuration of $\vecX$ and $\pi$ with $N=8$.
	\label{fig:pi8}}
\end{center}
\end{figure}

\begin{defn}
\label{def:ellxy}
If $\vecx$ and $\vecy$ are in a common cycle of $\pi$, we say that $\vecx$
is \plainidx{connected to} $\vecy$; otherwise we say that $\vecx$ is not
connected to $\vecy$.  These are written\index{$\ctdto$}\index{$\vecx\ctdto\vecy$}
$$
	\vecx \ctdto \vecy \qqand
	\vecx \notctdto \vecy.
$$
In the former case, we write
$$
	\ellxypi
$$
for the smallest positive integer $a$ such that $\pi^a(\vecx) = \vecy$.  This
is the number of permutation jumps from $\vecx$ to $\vecy$.
\end{defn}

\begin{defn}
We may average cycle lengths and spatial cycle lengths (definition
\ref{def:ell}) over all points $\vecx$:
\begin{align*}
	\ellbar(\pi) &=
		\frac{1}{N} \sum_{\vecx \in \Lambda} \ellx(\pi) &
	\essbar(\pi) &=
		\frac{1}{N} \sum_{\vecx \in \Lambda} \essx(\pi).
\end{align*}
\end{defn}

Now, $\bbE[\ellx] = \bbE[\elly]$ for all $\vecx,\vecy\in\Lambda$, and both are
the same as $\bbE[\ellbar]$ as shown in the lemma below.  Experimentally,
however, as discussed at the beginning of chapter \ref{chap:MCMC}, we cannot
compute expectations over all $\pi\in\SN$.  We must content ourselves with
sample means over some sequence of $M$ permutations $\pi_1$, $\ldots$, $\pi_M$,
obtained in an MCMC simulation.  We may increase the sample size by a factor
$N$ (thereby decreasing the sample variance by a factor of $N$, by the central
limit theorem) if we average over all points $\vecx$.  For each permutation, we
may compute $\ellbar(\pi)$ by averaging $\ellx(\pi)$ over all $N$ point
positions, and likewise for $\essx(\pi)$ and $\essbar$.

\begin{lem}
For all $\vecx \in \Lambda$, $\bbE[\ellbar] = \bbE[\ellx]$.
\begin{proof}
The left-hand side is
\begin{align*}
	\bbE[\ellbar] &= \sum_{\pi\in\SN} P(\pi) \ellbar(\pi).
\end{align*}
Since both sums are finite, we have
\begin{equation}
\label{eqn:ellbar_sums_commute}
\begin{aligned}
	\bbE[\ellbar]
	&= \sum_{\pi\in\SN} P(\pi) \ellbar(\pi)
	= \sum_{\pi\in\SN}
		P(\pi)
		\frac{1}{N}
		\sum_{\vecx \in \Lambda}
		\ellx(\pi) \\
	&=
		\frac{1}{N}
		\sum_{\vecx \in \Lambda}
		\sum_{\pi\in\SN}
		P(\pi)
		\ellx(\pi)
	= \overline{\bbE[\ellx]}.
\end{aligned}
\end{equation}
The equality $\overline{\bbE[\ellx]} = \bbE[\ellx]$ follows from translation
invariance.
\end{proof}
\end{lem}

In summary, we have
\begin{align*}
	\ellxpi &= \min\{a > 0: \pi^a(\vecx) = \vecx\} &
	\essxpi &= \sum_{j=1}^{\ellxpi} \|\pi^j(\vecx)-\pi^{j-1}(\vecx)\| \\
%
	\bbE[\ellx] &= \sum_{\pi\in\SN} P(\pi) \ellx(\pi) &
	\bbE[\essx] &= \sum_{\pi\in\SN} P(\pi) \essx(\pi) \\
%
	\ellbar(\pi) &= \frac{1}{N} \sum_{\vecx\in\Lambda} \ellx(\pi) &
	\essbar(\pi) &= \frac{1}{N} \sum_{\vecx\in\Lambda} \essx(\pi) \\
%
	\bbE[\ellbar] &= \sum_{\pi\in\SN} P(\pi) \ellbar(\pi) &
	\bbE[\essbar] &= \sum_{\pi\in\SN} P(\pi) \essbar(\pi).
\end{align*}
The quantities on the second line are convenient for theoretical use; the
quantities on the fourth line (due to the larger sample size) are preferable
for experimental use.

We define\index{$\xi$} a \plainidx{correlation length} $\xi$ to be
\begin{align}
	\xi &= \essbar.
\end{align}

Computational details are in section \ref{sec:comput_cyclen_corrlen}.

%% ----------------------------------------------------------------
\section{Mean and maximum jump length}
\label{sec:theory_mjl}

\begin{defn}
Let\index{$j$}\index{$j_\vecx$}
\begin{align*}
	j_\vecx(\pi) &= \|\pi(\vecx)-\vecx\|
\end{align*}
be the length of the permutation jump starting at site $\vecx$; let
\begin{align*}
	j(\pi) &= \frac{1}{N} \sum_{i=1}^N j_{\vecx_i}(\pi).
\end{align*}
\end{defn}
The \plainidx{mean jump length} at site $\vecx$ is simply
\begin{align*}
	\bbE[j_\vecx] &= \sum_{\pi\in\SN} P(\pi) \|\pi(\vecx)-\vecx\|.
\end{align*}
By linearity of expectation, this is the same as the average over
all sites:
\begin{align*}
	\bbE[j(\pi)] &= \bbE[j_{\vecx_1}(\pi)].
\end{align*}
As was the case for $\ellbar$ and $\essbar$ in section
\ref{sec:theory_cyclen_corrlen}, we approximate the uncomputably large sum over
all $N!$ permutations by a random sequence of $M$ permutations, and the sample
mean is random.  By the central-limit theorem argument in section
\ref{sec:theory_cyclen_corrlen}, the variance of the sample mean of $j$ is a
factor of $N$ smaller than the variance of the sample mean of $j_\vecx$, since
the sample size is $MN$ instead of $M$.  Additional computational details
are discussed in section \ref{sec:comput_mjl}.

%% ----------------------------------------------------------------
\section{Fraction of sites in infinite cycles $f_I$}
\label{sec:theory_fI}

If\index{$f_I$} one wants to quantify the temperature-dependent onset of long
cycles (section \ref{sec:qual_long_cycles}), then one can define a random
variable which counts the fraction of sites in long cycles.  In the infinite
limit, one looks for infinite cycles.  One could define
\begin{align}
\label{eqn:f_I_infty}
	f_I(\infty) &= 1 - \sum_{k \ge 1} P(\ell_0 = k)
\end{align}
where $\ell_0(\pi)$ is as defined in section \ref{sec:theory_cyclen_corrlen}.
For finite volume, where simulations are done, the right-hand side of
\ref{eqn:f_I_infty} is always 0:  every site is in a cycle of some finite
length.  One might then define
\begin{align}
\label{eqn:f_I_N}
	f_I(N) &= 1 - \sum_{k = 1}^{\veps(N)-1} P(\ell_0 = k)
\end{align}
where $\veps(N)$ is a small fraction of $N$, say, $\sqrt{N}$.  (This is chosen
so that as $N\to\infty$, one obtains $f_I(\infty)$.)

For practical computation, \cite{bib:gru} begin by defining
$$
	\rho = \frac{N}{V},
$$
i.e. $\rho$ is the \plainidx{particle density}.  For $1 \le m \le n \le N$,
define
$$
	\varrho_{m,n}(\pi)
	= \frac{1}{V} \;
	\# \left\{ i=1, \ldots, N: m \le \ell_i(\pi) \le n \right\}
$$
This random variable, taking values between 0 and $\rho$, is the
\plainidx{density of sites} in cycles of specified length.  One may also
consider the related random variable
$$
	f_{m,n}(\pi)
	= \frac{1}{N} \;
	\# \left\{ i=1, \ldots, N: m \le \ell_i(\pi) \le n \right\}
$$
which is $\varrho_{m,n}/\rho$.  (Of course, on the unit lattice where $\rho=1$,
the two random variables are identical.)  This runs from 0 to 1 and is the
\plainidx{fraction of sites} in cycles of specified length.  For figure
\ref{fig:pi8} on page \pageref{fig:pi8}, we have $f_{2,3}(\pi) = 3/8$.  Then
$\bbE[f_{\veps(N), N}]$ matches equation \ref{eqn:f_I_N} as follows:
\begin{align*}
	\bbE[f_{\veps(N),N}(\pi)] &= \frac{1}{N}
		\bbE\left[ \#\{i=1, \ldots, N : \ell_i(\pi) \ge \veps(N)\} \right] \\
	&= \frac{1}{N} \bbE\left[
		\sum_{i=1}^N 1_{\ell_i(\pi) \ge \veps(N)}(\pi)\} \right]
	= \frac{1}{N} P(\ell_i(\pi) \ge \veps(N)).
\end{align*}
This is the same as $P(\ell_0(\pi) \ge \veps(N))$ by translation invariance.
Then
\begin{align*}
	P(\ell_0(\pi) \ge \veps(N))
	&= \sum_{k=\veps(N)}^N P(\ell_0=k)
	= 1 - \sum_{k=1}^{\veps(N)-1} P(\ell_0=k) = f_I(N).
\end{align*}

In practice, a single cutoff of the form $\veps(N)$ is not used; one estimates
the infinite-limit behavior in a different way.  To see how to do this, we next
invoke results of \cite{bib:suto2} and \cite{bib:buv09} regarding the behavior
of $\bbE[f_{1,k}]$ as a function of $k/N$ in the infinite limit.  The former,
\cite{bib:suto2}, applies in the non-interacting case; the latter,
\cite{bib:buv09}, applies in the Ewens case but with non-spatial permutations,
which are equivalent to random spatial permutations with $T=0$.  For the
non-interacting case, $\bbE[f_{1,k}]$ is a straight line of slope 1 as shown in
the upper left of figure \ref{fig:f_I_why}.  At $T=0$, it fills the full
diagonal; cycle lengths have uniform distribution.  At $T>T_c$, the diagonal
vanishes into the upper-left corner; there are no long cycles.  The transition
to criticality occurs at $T$ such that the diagonal becomes visible.

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics{figures/f_I_why.eps}
\caption[Expectation of $f_{1,k}$ as a function of $k/N$.]
	{Qualitative behavior of $\lim_{N\to\infty} \bbE[f_{1,k}]$ as a function of
	$k/N$ in the non-interacting (upper left) and non-spatial large-$\alpha$
	Ewens-interacting (upper right) cases.  Lower left and lower right show the
	behavior for finite $N$.
	\label{fig:f_I_why}}
\end{center}
\end{figure}

In finite volume, $\bbE[f_{1,k}]$ is rounded as shown in the bottom left of
figure \ref{fig:f_I_why}.  Thus, one wishes to draw a tangent-line
approximation for the infinite-volume behavior, and take $f_I$ to be one minus
the vertical intercept.  This avoids use of a specific, arbitrary cutoff
$\veps(N)$, replacing it instead with a graphical estimator which makes use of
all available data.  See also section \ref{sec:comput_fI} for computational
details, including the handling of sampling variability.

For Ewens interactions, the diagonals are curved as shown in the upper right of
figure \ref{fig:f_I_why}.  For the small $\alpha$ values considered in this
dissertation, however, the Betz-Ueltschi-Velenik deviation of the curve from
pure linearity is less than the sampling variability in the data itself (see
the plots in section \ref{sec:comput_fI}).  Thus, the sag is not relevant to
our discussion.  It should also be remarked that the author has performed
larger-$\alpha$ simulations for which the $\bbE[f_{1,k}]$ sag is indeed
observed simulationally.

%% ----------------------------------------------------------------
\section{Macroscopic-cycle quotient $f_M/f_I$}
\label{sec:theory_lmax_GRU_quot}

\begin{defn}
For\index{$\ellmax$} a permutation $\pi$ in $\SN$, define
\begin{align}
	\ellmax(\pi) &= \max_{1 \le i \le N} \ell_i(\pi).
\end{align}
For a spatial permutation, this is precisely the same as
\begin{align}
	\ellmax(\pi) &= \max_{\vecx \in \Lambda} \ell_\vecx(\pi).
\end{align}
We write
\begin{align}
	f_M &= \bbE[\ellmax]/N.
\end{align}
\end{defn}

\begin{defn}
The \emphidx{macroscopic cycle quotient}, written $f_M/f_I$ for brevity,
is given by
\begin{align}
	\textrm{macroscopic cycle quotient}
	&= \begin{cases}
		\frac{\bbE[\ellmax]}{N f_I}, &\quad f_I \ne 0 \\
		0,                           &\quad f_I  =  0.
	\end{cases}
\end{align}
\end{defn}
Intuition was discussed in section \ref{sec:conj}; computational details are
discussed in section \ref{sec:comput_lmax_GRU_quot}.

% %% ----------------------------------------------------------------
% \section{Existence of infinite cycles}
% \label{sec:inf_cyc_existence}
%
% Here\index{infinite cycles} we describe theorem 1 of \cite{bib:qmath}, which
% is proved in section 1 of \cite{bib:bu07}.  This applies to the $U\equiv 0$
% case.
%
% %Here we describe theorem 1 of \cite{bib:qmath}, deferring [yyy?] a sketch
% %of the proof to section [yyy].
%
% \textbf{Thermodynamic limit:}  We\index{thermodynamic limit} inquire about
% the fraction of sites participating in short and long cycles (as quantified
% below) in the \plainidx{infinite-volume limit}.  Namely, we let $V,N \to
% \infty$ with fixed ratio $\rho = N/V$, and we ask about the\index{density}
% cycle-length distribution as a function of $\rho$.
%
% One does not need to construct an infinite-volume model, although this is
% done in section 3, for pure interest:  We will examine limits of expectations
% of random variables, where the limit is taken as the number of points $N$ of
% the model goes to infinity.  The limits then are in $\R$.
%
% \textbf{Critical density:}
% We define the \plainidx{critical density}, $\rhocz$, by the following
% formula.  Here, this is simply an \plainidx{ansatz}.
% \begin{equation}
% \label{eqn:ideal_rho_c}
% 	\rhocz = \int_{\R^d} \frac{d\veck}{e^{4\beta\pi^2\|\veck\|^2}-1}.
% \end{equation}
%
% \begin{rem}
% One may compute
% \begin{align}
% \label{eqn:ideal_rho_c_zeta}
% 	\rhocz &= \int_{\R^d} \frac{d\veck}{e^{4\beta\pi^2\|\veck\|^2}-1}
% 	= \frac{\zeta(3/2)}{(4\pi\beta)^{3/2}}
% \end{align}
% where $\zeta(z)$ is the \plainidx{Riemann zeta function}.
% \end{rem}
%
% % yyy Note that $\exp(-4\beta\pi^2|\veck|^2)$ is the \plainidx{Fourier %
% transform} of $\exp(\beta|\vecx|^2/4)$.
%
% % yyy wik bec article:
% % N = V int d3k/(2pi)^3 p(k)/(1-p(k))
% %   = V int d3k/(2pi)^3 1 / (exp(k^2/2mT) - 1)
% % w/ p(k) = exp(-k^2/2mT).
% % "maximum number of excited particles p/(1-p)".
%
% \textbf{Statement of the theorem:}
% For any $0<a<b<1$ (nominally, $a$ is just above 0 and $b$ is just below 1)
% and any $s \ge 0$,
% \begin{equation}
% \label{eqn:thm_1_rho}
% \begin{aligned}
% %
% 	\lim_{V \to \infty} \bbE_{\Lambda,N}(\varrho_{1,N^a}) &=
% 		\begin{cases}
% 		\rho,   & \rho \le \rhocz \\
% 		\rhocz, & \rhocz \le \rho
% 		\end{cases} \\
% %
% 	\lim_{V \to \infty} \bbE_{\Lambda,N}(\varrho_{N^a,N^b}) &=
% 		0 \\
% %
% 	\lim_{V \to \infty} \bbE_{\Lambda,N}(\varrho_{N^b,sN}) &=
% 		\begin{cases}
% 		0,           & \rho \le \rhocz \\
% 		\rho-\rhocz, & \rhocz \le \rho \le s+\rhocz \\
% 		s,           & s+\rhocz \le \rho.
% 		\end{cases}
% \end{aligned}
% \end{equation}
% (See figure \ref{fig:limErho}.)
%
% \begin{figure}[!htb]
% \begin{center}
% \psfragscanon
% \includegraphics{figures/limErho.eps}
% \caption[Distribution of cycle lengths as function of density.]
% 	{Distribution of cycle lengths as function of density in the
% 	infinite-volume limit.
% 	\label{fig:limErho}}
% \end{center}
% \end{figure}
%
% This statement, in terms of the density of sites in cycles of specified
% lengths, may be more intuitively rephrased in terms of the fraction of sites
% in
% cycles of specified length:
% \begin{align*}
% %
% 	\lim_{V \to \infty} \bbE_{\Lambda,N}(f_{1,N^a}) &=
% 		\begin{cases}
% 		1,   & \rho \le \rhocz\\
% 		\rhocz/\rho, & \rhocz \le \rho
% 		\end{cases} \\
% %
% 	\lim_{V \to \infty} \bbE_{\Lambda,N}(f_{N^a,N^b}) &=
% 		0 \\
% %
% 	\lim_{V \to \infty} \bbE_{\Lambda,N}(f_{N^b,sN}) &=
% 		\begin{cases}
% 		0,             & \rho \le \rhocz \\
% 		1-\rhocz/\rho, & \rhocz \le \rho \le s+\rhocz \\
% 		s/\rho,        & s+\rhocz \le \rho.
% 		\end{cases}
% \end{align*}
%
% The theorem has the following interpretation:
% \begin{itemize}
% \item At density below $\rhocz$, all sites are in \plainidx{short cycles} and
% no sites are in \plainidx{long cycles}.
% \item As density increases past $\rhocz$ and $\rhocz+s$, fewer sites are in
% short cycles and more sites are in long cycles.  In particular, a strictly
% positive fraction of sites are in long cycles.
% \item Asymptotically, all sites are in long cycles.
% \end{itemize}

%% ----------------------------------------------------------------
\section{Winding numbers, $f_S$, and $f_W$}
\label{sec:theory_wno}

The box $\Lambda=[0,L]^3$ with periodic boundary conditions is topologically
equivalent to the 3-torus.  Permutation cycles wind around the 3-torus some
number of times in the $x$, $y$, and/or $z$ directions.  The sub-$T_c$ onset of
long cycles corresponds to the appearance of cycles which wrap around the torus
in one or more of the three axes.  If a cycle goes around once in the clockwise
direction, we want to say it has sign $+1$; likewise, we want sign $-1$ for the
counterclockwise direction.  The following definition formalizes this
intuition.

\begin{defn}
\label{def:wno}
The \emphidx{winding number} (really a 3-tuple of numbers) of a permutation
$\pi$ is
\begin{align}
\label{eqn:def_wno}
	\vecW(\pi) &= (W_x(\pi), W_y(\pi), W_z(\pi))
	= \frac{1}{L} \sum_{i=1}^N \tilde{\vecd}(\pi(\vecx_i), \vecx_i).
\end{align}
This simply counts the integer number of wraps of $\pi$'s cycles around the
3-torus in each of the three directions.  We also write
\begin{align}
\label{eqn:def_W2}
	\vecW^2(\pi) &= \vecW(\pi) \cdot \vecW(\pi)
	= W_x(\pi)^2 + W_y(\pi)^2 + W_z(\pi)^2.
\end{align}
\end{defn}

\begin{rem}
\label{rem:differences_on_torus}
The formula for $\vecW(\pi)$ is more subtle than it first appears.  We work in
the regime where single permutation jumps are on the order of 1, 2, 3, or so
--- specifically, jumps are of length less than $L/2$.  The differences
$\tilde{\vecd}$ in \ref{eqn:def_wno} are to be interpreted in terms of lifts to
a separate Euclidean chart for each site $\vecx_i$.  For example, let $L=10$.
If $\vecx=(0,0,8)$ and $\pi(\vecx) = (0,0,9)$ then
$\tilde{\vecd}(\pi(\vecx),\vecx)$ is $(0,0,1)$.  However, if $\vecx=(0,0,0)$ and
$\pi(\vecx) = (0,0,9)$, which is equivalent mod $L$ to $\pi(\vecx) = (0,0,-1)$,
then $\tilde{\vecd}(\pi(\vecx),\vecx)$ is $(0,0,-1)$ rather than $(0,0,9)$.
\end{rem}

\begin{defn}
The\index{$f_S$} \emph{scaled winding number} is
\begin{align*}
	f_S &= \frac{\bbE[\vecW^2] L^2}{3 \beta N}
	= \frac{\bbE[\vecW^2]}{3 \beta L}.
\end{align*}
See \cite{bib:pc87} for the physical derivation.  For us, it simply needs to be
scaled by $L^2/N = 1/L$ in order to be an intensive parameter.
\end{defn}

\begin{defn}
Let $c_\vecx(\pi)$ consist of all sites in the same cycle as $\vecx$:
\begin{align*}
	c_\vecx(\pi) &= \{\vecy \in \Lambda:
		\vecy = \pi^a(\vecx), a = 0, 1, 2, \ldots, N-1\}.
\end{align*}
Let $\vecw_\vecx(\pi)$ be the \emphidx{winding vector} for $\vecx$ (which is
clearly the same for all sites in the same cycle as $\vecx$):
\begin{align*}
	\vecw_\vecx(\pi) &=
	\frac{1}{L} \sum_{\vecy \in c_\vecx(\pi)} \tilde\vecd(\pi(\vecy) - \vecy),
\end{align*}
where distances are again interpreted as in remark
\ref{rem:differences_on_torus}.  Note that all three slots of
$\vecw_\vecx(\pi)$ are necessarily integer-valued.  We say that $\pi$
\emphidx{winds through} $\vecx$ if $\vecx$ has a non-zero winding vector:
\begin{align*}
	t_\vecx(\pi) &= \begin{cases}
		1, &\quad \vecw_\vecx(\pi) \ne (0,0,0); \\
		0, &\quad \vecw_\vecx(\pi)  =  (0,0,0).
	\end{cases}
\end{align*}
We use these to define the \emphidx{fraction of sites in winding cycles}:
\begin{align*}
	f_W(\pi) &= \frac{1}{N} \sum_{\vecx \in \Lambda} t_\vecx(\pi).
\end{align*}
\end{defn}
Computational details are discussed in section \ref{sec:comput_wno}.

%% ----------------------------------------------------------------
\section{Order parameters: quantifying long cycles}
\label{sec:ordps}

Of the random variables presented in this chapter, the following, referred to
as \emph{order parameters}, may be used to locate the critical temperature
$T_c(\alpha)$ below which long cycles begin to appear.  The first four are
non-zero for $T<T_c$ and zero for $T>T_c$.  Since $\xi$ blows up below $T_c$,
$1/\xi$ goes to zero below $T_c$.  See also figure
\ref{fig:intuit_order_params}.
\begin{itemize}
\item Fraction of sites in infinite cycles $f_I$.
\item $f_M := \bbE[\ellmax]/N$.
\item Scaled winding number $f_S$.
\item Fraction $f_W$ of sites in cycles which wind.
\item Reciprocal correlation length $1/\xi$.
\end{itemize}

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.6]{../tests/fss3a/disfigs/rv_overview_alpha0_0.000_fI.eps}
\includegraphics[scale=0.6]{../tests/fss3a/disfigs/rv_overview_alpha0_0.000_fM.eps}

\includegraphics[scale=0.6]{../tests/fss3a/disfigs/rv_overview_alpha0_0.000_fS.eps}
\includegraphics[scale=0.6]{../tests/fss3a/disfigs/rv_overview_alpha0_0.000_fW.eps}

\includegraphics[scale=0.6]{../tests/fss3a/disfigs/rv_overview_alpha0_0.000_recip_xi.eps}
\caption[Behavior of order parameters as functions of $L$ and $T$.]
	{Behavior of order parameters as functions of $L$ and $T$, for the
	non-interacting model.  Each of the following occurs at a critical
	temperature $T_c$, in the limit $L \to \infty$: onset of long cycles
	($f_I$), onset of $N$-scaling of longest cycle length ($f_M$), onset of
	winding cycles ($f_S$ and $f_W$), and blow-up of correlation length
	(vanishing of $1/\xi$).  For finite $L$, the transitions are smooth.
	Interactions increase the critical temperature, shifting these graphs to
	the right.
	\label{fig:intuit_order_params}}
\end{center}
\end{figure}

%% ================================================================
\chapter{Markov chain Monte Carlo methods}
\label{chap:MCMC}

In this chapter we discuss the need for random-sampling methods, and justify
their use rigorously.  Given a \plainidx{random variable} $X(\pi)$, such as any
of those presented in chapter \ref{chap:RV_theory}, the
\plainidx{expectation} of $X$ is (equation \ref{eqn:E_pi})
\begin{align*}
	\bbE[X] &= \sum_{\pi\in\SN} P(\pi) X(\pi).
\end{align*}
This is a real number, with no uncertainty.  The problem is that the number of
permutations, $N!$, grows intractably in $N$:  even for $L=10$ (and we consider
$L$ up to $80$), $N=1000$ and $N!$ is a number with over 5,000 digits.  The
true expectation is effectively incomputable.  Expectations are instead
\emph{estimated}\index{estimator} by summing over some number $M$ (in the
current work, $10^5$ or $10^6$) of typical permutations.  The sample mean
\begin{align}
\label{eqn:pigmX}
	\pigm X &= \frac{1}{M} \sum_{i=1}^M X(\pi_i)
\end{align}
depends on the random sequence $\pi_1, \ldots, \pi_M$.  It is now a random
variable with its own variance.

To create such a sequence of system states (for us, permutations), the method
used throughout the computational physics community
\cite{bib:berg,bib:landaubinder} --- and this dissertation is no exception ---
is Markov chain Monte Carlo.  Namely, given a permutation $\pi_i$, one selects
a successor permutation $\pi_{i+1}$ in some random way.  This is the Monte
Carlo part.  Moreover, the transition probabilities from $\pi_i$ to each
candidate $\pi_{i+1}$ depend only on $\pi_i$, and not on any previous choices.
This is the Markov chain part.

In the next sections we show (1) we can construct Markov chains which sample
from the Gibbs distribution $P(\pi)$ (equation \ref{eqn:RCM_Gibbs_P}); (2)
other distributions (induced by the selection of initial state) converge to the
Gibbs distribution; (3) the sample mean $\pigm X$ converges almost surely to
the true expectation $\bbE[X]$; and (4) the variance of $\pigm X$ can be
estimated, allowing us to place error bars on our estimates of $\bbE[X]$.

Most of the material in this chapter is familiar: see
\cite{bib:berg,bib:landaubinder,bib:cb,bib:fg,bib:gs} to name only a few.
Results are restated here for self-containment of presentation, and they are
phrased in the language of random permutations for convenience.

%% ----------------------------------------------------------------
\section{Invariant and limiting distributions}
\label{sec:invt}

The only probability distribution for random permutations considered thus far
is the Gibbs distribution of equation \ref{eqn:RCM_Gibbs_P} on page
\pageref{eqn:RCM_Gibbs_P}.  However, others exist.  If one is asked for a
permutation $\pi$ and one always replies with the identity, then the
distribution of that answer is the singleton measure on the identity.  This is
not the Gibbs measure (unless $\alpha = 0$ and $\beta = 0$).  A third
distribution is the uniform distribution, where each permutation has
probability $1/N!$.  (This is the same as the Gibbs measure only for $\alpha =
0$ and $T=0$.) We will continue to write the Gibbs measure as $P(\pi)$; we will
use Greek letters such as $\lambda(\pi)$ or $\theta(\pi)$ for other
distributions.

\begin{defn}
Let $\Omega$ be a finite set with a sequence of probability measures
$\lambda_1$, $\lambda_2, \lambda_3, \ldots$.  (Since $\Omega$ is finite, the
$\sigma$-fields may be taken to be all of $2^\Omega$.)  Let $S_1, S_2, S_3,
\ldots$ be a random sequence of states where each $S_i \in \Omega$ is
distributed according to $\lambda_i$.  Let $\Lambda_k$ be the joint probability
distribution on $\Omega^k$.  That is, the $S_i$'s are not necessarily either
independent or identically distributed; the $\lambda_i$'s are the marginals of
$\Lambda_k$ on the $i$th slot.  We say the sequence $S_1, S_2, S_3, \ldots$ is
a \emphidx{Markov chain} if, for all $i=1,2,3,\ldots$ and all $\omega_i \in
\Omega$,
\begin{align*}
	\Lambda_{i+1}(S_{i+1}=\omega_{i+1} \mid S_1=\omega_1, \ldots, S_i=\omega_i)
	&=
	\Lambda_{i+1}(S_{i+1}=\omega_{i+1} \mid S_i=\omega_i).
\end{align*}
That is, in a Markov chain, the transition probabilities from the current state
to the next depend only on the current state and not on any deeper history.
Furthermore, if for all $i \ge 1$,
\begin{align*}
	\Lambda_{i+1}(S_{i+1} = \omega_{i+1} \mid S_{i} = \omega_i)
	&=
	\Lambda_{2}(S_2 = \omega_{i+1} \mid S_1 = \omega_i),
\end{align*}
then we say the Markov chain is \emphidx{homogeneous}.
(All Markov chains considered in this dissertation are homogeneous.  For us, of
course, $\Omega$ is $\SN$ and the states $S_i$ are permutations $\pi_i$.)
\end{defn}

\begin{defn}
The \emphidx{Markov matrix} associated with a homogeneous $K$-state Markov
chain is
\begin{align*}
	\vecM &= (M_{ij})_{1 \le i,j \le K}
	= \Lambda_2(S_2 = \omega_j \mid S_1 = \omega_i).
\end{align*}
We also write
\begin{align*}
	M(S_2 = \omega_j \mid S_1 = \omega_i)
	&= \Lambda_2(S_2 = \omega_j \mid S_1 = \omega_i).
\end{align*}
(For us, of course, $K=N!$.)
\end{defn}

Note that a Markov chain is determined by the probability distribution
$\lambda_1$ for its initial state, along with the Markov transition matrix
$\vecM$.  To see this, use the definition of  conditional probability to split
the joint distribution into a sequencing of states:
\begin{align*}
	\Lambda_k&(S_1 = \omega_1,S_2 = \omega_2, \ldots, S_k = \omega_k)
	= \Lambda_k(S_1=\omega_1)
	\cdot \Lambda_k(S_2 = \omega_2 \mid S_1 = \omega_1) \\
	&\cdots
	\Lambda_k(S_k = \omega_k \mid S_1 = \omega_1, \cdots, S_{k-1}=\omega_{k-1}) \\
	=& \Lambda_k(S_1=\omega_1)
	\cdot \Lambda_k(S_2 = \omega_2 \mid S_1 = \omega_1)
	\cdots \Lambda_k(S_k = \omega_k \mid S_{k-1}=\omega_{k-1}) \\
	=& \lambda_1(S_1=\omega_1)
	\cdot M(S_2 = \omega_2 \mid S_1 = \omega_1)
	\cdots M(S_k = \omega_k \mid S_{k-1}=\omega_{k-1}).
\end{align*}

\begin{defn}
A probability distribution $\theta(S)$ is \emphidx{invariant} for a Markov
transition matrix $\vecM$ if $\theta = \theta \vecM$, that is, if for all $i=1,
\ldots, K$,
\begin{align*}
	\theta(S_2 = \omega_j) &= \sum_{i=1}^K
		M(S_2 = \omega_j \mid S_1 = \omega_i)
		\theta(S_1 = \omega_i).
\end{align*}
In other words, $\theta$ is invariant for $\vecM$ if $\vecM$-transitions
preserve the distribution $\theta$.  If the distribution of $S_1$ is $\theta$,
then $S_2$ will also be distributed according to $\theta$, and so on.  Such a
sequence of states $S_1, S_2, \ldots$, while not in general independent, is
identically distributed.
\end{defn}

\begin{rem}
A homogeneous chain is not the same as a stationary sequence.  In the former
cases, the transition matrix is the same at each time step; in the latter case,
the probability distributions are the same at each time step.
\end{rem}

\beginex An illustrative example uses die-tipping.  Recall that an ordinary
six-sided die has pips on opposite faces summing to seven.  There are six
states, which we assume to be uniformly distributed if the die is rolled.
Given that the die has $n$ pips facing up, we may \emph{tip} the die by picking
one of the four sides at uniform random and putting that side up.  E.g. if 1 is
up, then after the tip, 2, 3, 4, or 5 will appear each with probability $1/4$;
1 or 6 will appear with probability $0$.

The transition matrix, with rows indexing the current state and columns
indicating the successor state, is
\begin{align*}
	\vecM = \prmatrix{
	0   & 1/4 & 1/4 & 1/4 & 1/4 & 0   \\
	1/4 & 0   & 1/4 & 1/4 & 0   & 1/4 \\
	1/4 & 1/4 & 0   & 0   & 1/4 & 1/4 \\
	1/4 & 1/4 & 0   & 0   & 1/4 & 1/4 \\
	1/4 & 0   & 1/4 & 1/4 & 0   & 1/4 \\
	0   & 1/4 & 1/4 & 1/4 & 1/4 & 0
	}.
\end{align*}
Suppose the die is initially set on table with $1$ up, i.e. $\lambda_1$
is the singleton distribution with $\lambda_1(1)=1$.  Then
\begin{align*}
	\lambda_2 &= \prmatrix{
		0   & 1/4 & 1/4 & 1/4 & 1/4 & 0
	}
\end{align*}
If, instead, the die is initially rolled, then $\lambda_1$ is uniform.  One
computes $\lambda_2$ to be uniform as well.  The 1-up distribution is not
invariant for the die-tipping transition rule, but the die-roll distribution
is.
\mendex

Given a Markov matrix $\vecM$, one may wish to find an invariant distribution
$\theta$.  Going the other way, given a distribution $\theta$, one may wish to
construct a Markov matrix $\vecM$ such that $\theta$ is invariant with respect
to $\vecM$.  The latter is our main goal here:  the distribution of interest is
the Gibbs distribution $P$.  The next theorem is key.  First, we define the
terminology necessary to state it.

\begin{defn}
A Markov chain $(\vecM, \lambda_1)$ on a state space $\Omega$ is
\emphidx{irreducible} if for all $S, S' \in \Omega$ there exists an $n$ such
that $M^n(S,S') > 0$.
\end{defn}

\begin{defn}
\label{defn:period}
In a Markov chain $(\vecM,\lambda_1)$ on a state space $\Omega$, a state
$\omega\in\Omega$ has \emphidx{period} $p$ if any return to $S$ must occur at
multiples of $p$ steps.  A state is \emphidx{aperiodic} if $p=1$.  Precisely,
\begin{align*}
	p = \gcd\{n : \lambda(S_n = \omega \mid S_0 = \omega) > 0\}.
\end{align*}
Furthermore, the chain is said to be aperiodic if all states are aperiodic.
\end{defn}

\begin{defn}
A Markov matrix $\vecM$ on a state space $\Omega$ (with $\#\Omega=K$) and a
distribution $\theta$ on $\Omega$ are \emphidx{reversible}, or satisfy
\emphidx{detailed balance}, if for all $1 \le i, j \le K$,
\begin{align*}
	M_{ij} \theta(S_i) &= M_{ji} \theta(S_j).
\end{align*}
\end{defn}

\begin{rem*}
An irreducible, aperiodic chain on a finite state space is sometimes called
\emphidx{ergodic}.
\end{rem*}

\begin{thm}[Invariant-distribution theorem]
\label{thm:irr_aper_DB_invt}
Fix a finite set $\Omega$.  Let $\vecM$ be a Markov transition matrix on
$\Omega$, as above, and let $\theta$ be a probability distribution on $\Omega$.
If the Markov chain $(\vecM, \theta)$ is \plainidx{irreducible},
\plainidx{aperiodic}, and satisfies \plainidx{detailed balance}, then $\theta$
is invariant for $\vecM$.  That is, the Markov chain $(\vecM,\theta)$ is
invariant.  Furthermore, if the Markov chain $(\vecM, \lambda_1)$ is
irreducible, aperiodic, and satisfies detailed balance, then for all $S \in
\Omega$, $\lambda_n(S) \to \theta(S)$ as $n \to \infty$.
\begin{rem*}
Some authors call this this \emphidx{ergodic theorem}; yet, others call our
theorem \ref{thm:sampling} the ergodic theorem.
\end{rem*}
\begin{proof}
See \cite{bib:lawler}, sections 1.2 and 1.3.
\end{proof}
\end{thm}

We say $\theta$ ($P$ in the context of random permutations) is the
\emphidx{stationary distribution} or \emphidx{invariant distribution} for
$\vecM$.  We say that it is also a \emphidx{limiting distribution} for any
initial distribution $\lambda_1$ satisfying the above hypotheses.  Importantly,
for simulations, we need not know what a typical permutation looks like.  We
may start always with the identity permutation, i.e. $\lambda_1$ is the
singleton distribution supported on the identity permutation.  We may then run
the Markov chain, producing a sequence of permutations.  As the number of
iterations goes to infinity, the distribution of permutations approaches $P$.

The specific number $M$ of iterations needed for convergence of $\lambda_M$ to
$\pi$ is another matter entirely.  The theory given by Lawler, as noted above,
guarantees that a Markov matrix $\vecM$ with $P$ as its invariant distribution
has no other invariant distribution:  $\vecM$ has a single eigenvalue 1 with
eigenvector $P$.  The rate of convergence of $\lambda_n$ to $P$ depends on the
second-largest eigenvalue for $\vecM$, which one in general does not know how
to compute.  In practice, this \emphidx{mixing time} or \emphidx{thermalization
time} is estimated, using techniques such as those in section
\ref{sec:thermalization}.

The theory above also does not tell us how to construct a Markov matrix $\vecM$
having a desired distribution $P$ as its invariant distribution --- it simply
tells us what we can do once we have constructed such a matrix.  A specific
technique is due to N. Metropolis \cite{bib:berg,bib:landaubinder,bib:cb}.  The
essence is that if the invariant probability distribution $P$ is defined as a
Gibbs measure via an energy function $H$ on $\Omega$, then proposed successor
states $S_{i+1}$ are accepted with probability $\min\{1, e^{-\Delta H}\}$ where
$\Delta H$ is the energy difference for the state change.  In this
dissertation, I directly prove that such methods result in detailed balance.
Thus, the reader is referred to propositions \ref{prop:SO_DB} and
\ref{prop:worm_DB} for details.

%% ----------------------------------------------------------------
\section{The sampling theorem}
\label{sec:sampling_thm}

Here we show that the sample mean $\pigm X$ (equation \ref{eqn:pigmX} on page
\pageref{eqn:pigmX}) converges to the true mean $\bbE[X]$ (equation
\ref{eqn:E_pi} on page \pageref{eqn:E_pi}).

\begin{thm}[Sampling theorem]
\label{thm:sampling}
Let $X$ be a random variable on the finite probability space $(\Omega,2^\Omega,
\theta)$.
If the stationary Markov chain $(\vecM,\theta)$ satisfies the hypotheses
of theorem \ref{thm:irr_aper_DB_invt}, then
\begin{align*}
	\frac{1}{M} \sum_{i=1}^M X(S_i) &\to \bbE[X]
	\quad\textrm{as}\quad
	M \to \infty.
\end{align*}
\begin{rem*}
Some authors call this this \emphidx{ergodic theorem}; yet, others call our
theorem \ref{thm:irr_aper_DB_invt} the ergodic theorem.
\end{rem*}
\begin{proof}
See \cite{bib:fg}, section 28.4.  Note that, trivially, if the distribution on
states $\omega$ is invariant, then so is the distribution on real-value random
variable $X(\omega)$.
\end{proof}
\end{thm}

% yyy maybe prove that?

%% ----------------------------------------------------------------
\section{Sample means and their variances}
\label{sec:body_tauint}

There is one caveat to replacing the true mean $\bbE[X]$ with the MCMC sample
mean $\pigm X$:  the naive computation of the standard deviation of the sample
mean, which is correct for independent identically distributed states, is a
wildly incorrect underestimate for the standard deviation of the sample mean in
the case of identically distributed but correlated states.  This issue is
so important that appendix \ref{appx:eta} is devoted to it.
The key result of the appendix is that the standard deviation (error bar) of
the sample mean $\pigm X$ is off by a factor of the \emphidx{integrated
autocorrelation time} $\sqrt{\hattauint(X)}$, which is computed as described in
section \ref{sec:estimate_meanvar}.

%% ----------------------------------------------------------------
\section{Simple example: 1D Ising}
\label{sec:1D_Ising_metro}

For an example which is more complex than die-tipping but less complex than
random spatial permutations, consider the 1D $N$-point Ising model.  Namely,
the configuration space is $\Omega=\{\pm 1\}^N$, i.e. $N$ particles which may
be in either an up (filled) or a down (hollow) state:
\begin{center}
\psfragscanon
\includegraphics[scale=1.0]{figures/ising1.eps}
\end{center}
A state is described by
$
	\vecomega = (\omega_1, \ldots, \omega_n).
$
The configuration space $\Omega$ has $2^N$ possible configurations.  The system
is endowed with an \plainidx{energy function}.  For the 1D Ising model, one has
$$
	H(\vecomega) = \sum_{i=1}^n \sum_{j=1}^n S_{ij} \omega_i \omega_j
		+ \sum_{i=1}^n h_i \omega_i.
$$
where the $S_{ij}$'s are interaction terms (non-interacting, nearest neighbor,
mean-field, etc.) and the $h_i$'s are magnetization terms.  Given a
\plainidx{temperature}-related parameter $\beta$, one sets the (Gibbs)
probability of each configuration to be proportional to $e^{-\beta H}$.

One picks an \plainidx{initial configuration}.  There are three obvious
choices: (1) Start with all spins down, i.e. $\vecomega = (-1, \ldots, -1)$.
(2) Start with all spins up, i.e. $\vecomega = (+1, \ldots, +1)$.  (3) Start
with $\vecomega$ selected from a uniform probability distribution on $\Omega$.
Then, one selects a \plainidx{site} $i$ and decides whether to flip $\omega_i$ to
$-\omega_i$.

\begin{center}
\psfragscanon
\includegraphics[scale=1.0]{figures/ising2.eps}
\end{center}

This decision is made using the \plainidx{Metropolis} prescription, namely:
\begin{itemize}
\item One computes the change in energy $\Delta H = H(\vecomega')
	- H(\vecomega)$ which would be obtained if $\vecomega$ were
	sent to $\vecomega'$ by flipping $\omega_i$.
\item One may compute $\Delta H$ by separately computing
	$H(\vecomega')$ and $H(\vecomega)$ and subtracting the two.
	However, since the only change is at the site $i$, one may do some
	ad-hoc algebra to derive an expression for $\Delta H$ which is less
	computationally expensive.
\item One accepts the change with probability $\min\{1, e^{-\Delta H}\}$.
\end{itemize}
This is called a \emphidx{Metropolis step}\index{step}.

Looping through all $n$ sites from $i=1$ to $i=n$, performing a Metropolis step
at each site $i$, is called a \emphidx{Metropolis sweep}\index{sweep}.  If one
realizes a random variable $X(\vecomega)$ at each of $M$ sweeps, averaging $X$
over the $M$ sweeps, one obtains an approximation $\pig{X}_M$ for the
expectation $\bbE[X]$.  One should first run $L$ Metropolis sweeps of the
system, discarding the realizations of the random variable $X$, before running
the $M$ sweeps in which data are accumulated.  The $L$ sweeps are called the
\plainidx{thermalization phase}; the $M$ sweeps are called the
\plainidx{accumulation phase}.

%% ----------------------------------------------------------------
\section{Recipe for MCMC algorithms}
\label{sec:MCMC_recipe}

The naive outline of an MCMC run is simple:
\begin{itemize}
\item Use a Markov chain, discussed at the end of this section,
to generate a sequence $\pi_1, \pi_1, \ldots, \pi_M$ of permutations.
\item For each permutation $\pi_i$, for each random variable $X$ of interest,
remember the value $X_i = X(\pi_i)$.
\item Compute the sample mean $\Xbar = \frac{1}{N}\sum_{i=1}^N X_i$.
Also compute the sample standard deviation, and any other desired
statistics.
\item Display the statistics.
\end{itemize}
Since the initial permutation $\pi_1$ is the identity, the initial distribution
is the singleton supported at the identity, which is not the Gibbs distribution
$P$.  Furthermore, as a very low-energy state, the identity highly non-typical
with respect to the Gibbs distribution (equation \ref{eqn:RCM_P}) for the model
of spatial permutations.  As discussed at the end of section \ref{sec:invt},
one runs the chain until it is thermalized.  (The number of steps $\tau$
required for this is random, but it turns out to fall within a narrow range.)
Renumbering $\pi_\tau$ to $\pi_0$, one then accumulates statistics over the $M$
permutations $\pi_1, \ldots, \pi_M$.  See section \ref{sec:thermalization} for
the thermalization-detection algorithm used in this dissertation.

Thus the computational recipe is as follows:

\begin{itemize}
\item Start with the initial permutation being the identity permutation.
\item Run the Markov chain, generating a sequence of permutations until
	thermalization has been detected.  At that point, call the current
	permutation $\pi_0$.  \item Continue generating a sequence $\pi_1, \ldots,
	\pi_M$ of permutations.
\item For each permutation $\pi_i$, for each random variable $X$ of interest,
remember the value $X_i = X(\pi_i)$.
\item After all $M$ permutations have been generated, compute the sample mean
	$\Xbar = \frac{1}{N}\sum_{i=1}^N X_i$.  Also compute the sample standard
	deviation and its error bar (using estimated integrated autocorrelation
	time), and any other desired statistics, such as histograms.
\item Display the statistics.
\end{itemize}

Given the framework established by previous sections of this chapter, the
recipe to prove correctness of this algorithm reduces to the following: define
a Markov chain, then prove irreducibility, aperiodicity, and detailed balance.
We devote the next chapters to present three Markov chains:  the swap-only
algorithm, the swap-and-reverse algorithms, and the worm algorithm.

%% ================================================================
\chapter{The swap-only and swap-and-reverse algorithms}
\label{chap:SO_SAR}

This chapter and chapter \ref{chap:worm} present Markov chains for MCMC
sampling of the model of random spatial permutations (chapter
\ref{chap:model_defns}) within the framework of chapter \ref{chap:MCMC}.  The
chains are proved correct, then compared and contrasted.  Computational results
using the swap-and-reverse algorithm are given in chapter \ref{chap:results}.

%% ----------------------------------------------------------------
\section{The swap-only algorithm}
\label{sec:SO_algorithm}


The \plainidx{swap-only algorithm} for transitioning from $\pi$ to $\pi'$,
within the context of the recipe in section \ref{sec:MCMC_recipe}, is as
follows:
\begin{itemize}
\item Call the current permutation $\pi$.
\item Sweep through sites $\vecx$ of the
lattice in either lexical or uniform-random order.  \item For each site
$\vecx$, do a \emphidx{Metropolis step}:
	\begin{itemize}
	\item Choose a site $\pi(\vecy)$ from among the six nearest
	neighbors of $\pi(\vecx)$.
	\item Propose to change $\pi$ to the permutation $\pi'$ which
		has $\pi'(\vecz)=\pi(\vecz)$ for all $\vecz \ne \vecx, \vecy$
		but $\pi'(\vecx) = \pi(\vecy)$ and $\pi'(\vecy) = \pi(\vecx)$.
		(See figure \ref{fig:soswap}.)
%		\begin{align*}
%		\pi: &
%		\pcmatrix{
%			\ldots & \vecx      & \ldots & \vecy      & \ldots \\
%			       & \downarrow &        & \downarrow &        \\
%			\ldots & \pi(\vecx) & \ldots & \pi(\vecy) & \ldots \\
%		} &
%		\pi': &
%		\pcmatrix{
%			\ldots & \vecx      & \ldots & \vecy      & \ldots \\
%			       & \downarrow &        & \downarrow &        \\
%			\ldots & \pi(\vecy) & \ldots & \pi(\vecx) & \ldots \\
%		} \\
%		\end{align*}

	\item With probability proportional to $\min\{1, \edh\}$
	where $\Delta H = H(\pi') - H(\pi)$, accept the change.
	(If the change is rejected, set $\pi'=\pi$.)
	\end{itemize}
\item After a sweep, obtain a value of random variable(s) for inclusion in
computation of their sample means.  Label the permutations at the start and
ends of the sweep $\pi_i$ and $\pi_{i+1}$, respectively.  All the intermediate
$\pi'$ are not retained.
\end{itemize}

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.7]{figures/gkswap.eps}
\caption{Metropolis moves for the swap-only algorithm.
	\label{fig:soswap}}
\end{center}
\end{figure}

\begin{defn}
\label{defn:SO_trivial}
A swap is \emphidx{trivial} if $\vecx = \vecy$.
\end{defn}

%% ----------------------------------------------------------------
\section{Explicit construction of the Markov matrix}
\label{sec:SO_Markov_matrix}

For section \ref{sec:SO_correctness} we will need an explicit construction of
the Markov matrix corresponding to the swap-only algorithm as described in
section \ref{sec:SO_algorithm}.

The Markov perspective on the algorithm is that the distribution
$\lambda_0(\pi)$ of the first permutation is either supported solely on the
identity, or uniform on all $N!$ permutations.  The distribution for subsequent
permutations is
\begin{align*}
	\lambda_{i+1}(\pi') &= \sum_{\pi \in \SN} \lambda_{i}(\pi) M(\pi, \pi')
\end{align*}
or, in matrix/vector notation,
\begin{align*}
	\lambda_{i+1} &= \lambda_{i} \vecM.
\end{align*}
In this section we precisely describe the matrix $\vecM$; in section
\ref{sec:SO_correctness} we show that $\lambda_i$ approaches the Gibbs
distribution (equation \ref{eqn:RCM_P}).

The matrix $\vecM$ is $N! \times N!$: rows are indexed by $\pi_1, \ldots,
\pi_{N!}$ and columns are indexed by $\pi'_1, \ldots, \pi'_{N!}$.  Most of the
entries of $\vecM$ are zero:  Metropolis steps change only two permutation
sites whereas most $\pi,\pi'$ differ at more than two sites.

\begin{defn}
For $\pi,\pi'\in\SN$, define
\begin{align*}
	d(\pi,\pi') &= \#\{i=1,2,\ldots,N : \pi(i) \ne \pi'(i)\}.
\end{align*}
\end{defn}

\begin{rem*}
Note that $d(\pi,\pi') \ne 1$ since if two permutations agree on $N-1$
sites, they must agree on the remaining site.
\end{rem*}

\begin{lem}
\label{lem:d_metric}
The function $d(\pi,\pi')$ is a \plainidx{metric} on $\SN$.
\begin{proof}
Symmetry is obvious, as is non-negativity.  For positive definiteness, note
that $d(\pi,\pi')=0$ iff $\pi=\pi'$.  For the triangle inequality, let
$\pi,\pi',\pi'' \in \SN$.  Partition the set $\{1, 2, \ldots, N\}$ into the
four disjoint sets
\begin{align*}
	A &= \{i = 1,2,\ldots,N : \pi(i) =   \pi'(i), \pi'(i) =   \pi''(i)\}, \\
	B &= \{i = 1,2,\ldots,N : \pi(i) =   \pi'(i), \pi'(i) \ne \pi''(i)\}, \\
	C &= \{i = 1,2,\ldots,N : \pi(i) \ne \pi'(i), \pi'(i) =   \pi''(i)\}, \\
	D &= \{i = 1,2,\ldots,N : \pi(i) \ne \pi'(i), \pi'(i) \ne \pi''(i)\}.
\end{align*}
Then $\pi=\pi''$ on all of $A$;
$\pi\ne \pi''$ on all of $B$ and $C$;
and $\pi,\pi''$ may or may not agree on various elements of $D$:
$$
\begin{tabular}{|c||c||c||c|}
\hline $A$ & $B$ & $C$ & $D$ \\
\hline $\pi  = \pi' $ & $\pi  = \pi' $ & $\pi  \ne \pi' $ & $\pi  \ne \pi'$ \\
\hline $\pi' = \pi''$ & $\pi' \ne \pi''$ & $\pi' = \pi''$ & $\pi' \ne \pi''$ \\
\hline $\pi  = \pi''$ & $\pi \ne \pi''$ & $\pi  \ne \pi''$ & Varies \\
\hline
\end{tabular}
$$
That is,
\begin{align*}
	d(\pi, \pi')  &=   \#C + \#D, \\
	d(\pi',\pi'') &=   \#B + \#D, \\
	\#B + \#C &\le d(\pi, \pi'') \le \#B + \#C + \#D.
\end{align*}
Then
\begin{align*}
	d(\pi,\pi'') &\le \#B + \#C + \#D \le \#B + \#C + 2\#D
	= d(\pi,\pi') + d(\pi',\pi'').
\end{align*}
\end{proof}
\end{lem}

\begin{defn}
Lattice sites $\vecx,\vecy$ are \emphidx{nearest-neighbor} if
$\|\vecx-\vecy\| = 1$.
\end{defn}

%\begin{defn}
%For $\pi,\pi' \in \SN$, let
%$$
%	\Delta H = H(\pi') - H(\pi).
%$$
%If we write $\pi',\pi \in \SN$, then we mean
%$$
%	\Delta H = H(\pi) - H(\pi').
%$$
%\end{defn}

\begin{defn}
\label{defn:R_pi_and_ctdto}
For $\pi\in\SN$, define
\begin{align*}
	R(\pi) &= \{\pi'\in\SN: d(\pi,\pi')=2 \textrm{ and }
		\|\pi(\vecx)-\pi(\vecy)\|=1\}
\end{align*}
where the $\vecx$ and $\vecy$ are taken to be the two points at which
$\pi,\pi'$ differ.  Then $R(\pi)$ is the set of permutations $\pi'$ reachable
from $\pi$ on a swap.
We also write\index{$\ctdto$}\index{connected to}\index{$\pi'\ctdto\pi$}
$$
	\pi' \ctdto \pi
$$
if $\pi' \in R(\pi)$ or $\pi'=\pi$.
\end{defn}

The Metropolis steps are then described as follows.  For each $\pi\in\SN$,
\begin{align*}
	M(\pi,\pi') &= \begin{cases}
		C\left(1 \wedge e^{-H(\pi')+H(\pi)}\right),
		& \pi' \in R(\pi), \\
		%
		\displaystyle
		1 - \sum_{\pi' \in R(\pi)}
		C\left(1 \wedge e^{-H(\pi')+H(\pi)}\right),
		& \pi=\pi'; \\
		%
		0,
		& \textrm{ otherwise}.
	\end{cases}
\end{align*}

To find out what $C$ is:
\begin{itemize}
\item There are $N$ choices of lattice points $\vecx$.
\item For each $\vecx$, there are 6 choices of $\piy$ which are
nearest neighbors to $\pix$.
\item This double-counts the $3N$ distinct choices of $\pi'$ reachable from
$\pi$ in a single Metropolis step, since choosing $\vecx$ and then $\vecy$
results in the same Metropolis step as choosing $\vecy$ and then $\vecx$.
\end{itemize}

% yyy To do: prove $3N$ distinct.

We require row sums to be 1 for stochastic matrices, so we may take
\begin{align*}
	C &= \frac{1}{3N}.
\end{align*}
This prefactor of $1/3N$ is the same for all non-zero entries of $\vecM$
for swap-only Metropolis.  Then
\begin{align}
\label{eqn:SO_M}
	M(\pi,\pi') &= \begin{cases}
		\frac{1}{3N}\left(1 \wedge e^{-H(\pi')+H(\pi)}\right),
		& \pi' \in R(\pi), \\
		%
		\displaystyle
		1 - \sum_{\pi' \in R(\pi)}
		C\left(1 \wedge e^{-H(\pi')+H(\pi)}\right),
		& \pi=\pi'; \\
		%
		0,
		& \textrm{ otherwise}.
	\end{cases}
\end{align}

%% ----------------------------------------------------------------
\section{Correctness of the swap-only algorithm}
\label{sec:SO_correctness}

It is clear that the swap-only algorithm produces a sequence of permutations,
but with what distribution?  From the Markov-chain theory in section
\ref{sec:invt}, we know that if the chain is irreducible, aperiodic,
and satisfies detailed balance, then the chain has the Gibbs distribution
(equation \ref{eqn:RCM_P}) as its unique invariant distribution.

Also note from Markov-chain theory \cite{bib:lawler} that all states in a
recurrence class have the same period.  Thus, if we can show that the chain is
irreducible (i.e. the entire state space is a single recurrence class), then
for aperiodicity of the chain it suffices to show that a single state (e.g. the
identity permutation) has period 1.

\begin{prop}[Irreducibility]
\label{prop:SO_irreducible}
For all $\pi,\pi'$, there is an $n$ such that $M^n(\pi,\pi') > 0$.
That is, any permutation is reachable from any other.
\begin{proof}
Transpositions generate $\SN$:  for all $\pi\in\SN$, there exist transpositions
$\sigma_1$, $\ldots$, $\sigma_m$ such that
%\begin{align*}
$	\pi = \prod_{j=1}^m \sigma_j$.
%\end{align*}
Thus, it suffices to show that given any permutation $\pi$ and any two points
$\vecx$ and $\vecz$, so $\pi: \vecx \mapsto \pi(\vecx)$ and $\pi: \vecz \mapsto
\pi(\vecz)$, we can construct a sequence of swaps sending $\pi$ to $\pi'$ so
that $\pi': \vecx \mapsto \pi(\vecz)$, $\pi': \vecz \mapsto \pi(\vecx)$, and
$\pi'(\vecy) = \pi(\vecy)$ for all $\vecy \ne \vecx,\vecz$.  (If $\pi(\vecx)$
and $\pi(\vecz)$ are nearest-neighbor lattice sites, of course, then a single
swap does the job.)

Define $G_{\veca,\vecb}:\SN\to\SN$ to be the swap operator for
nearest-neighbor lattice sites $\pi(\veca)$ and $\pi(\vecb)$.  Write
$\pi'=G_{\veca,\vecb} \pi$.  Given $\vecx$ and $\vecz$, there is a (non-unique)
sequence of lattice sites $\vecy_0, \vecy_1, \vecy_2, \ldots, \vecy_n$ such
that $\vecy_0=\vecx$, $\vecy_n=\vecz$, and
$\|\pi(\vecy_{i+1})-\pi(\vecy_i)\|=1$ for $i=0, 1, \ldots, n-1$.  (See figure
\ref{fig:irrprop}.)  We will construct a sequence of swaps along this
nearest-neighbor path whose end result is to swap the permutation arrows
starting at $\vecx$ and $\vecz$, leaving all other arrows unchanged.  We first
need a lemma about compositions of swaps.

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.6]{figures/irrprop.eps}
\caption[Sequence of local swaps resulting in a non-local swap.]
	{A sequence of (nearest-neighbor) swaps which results in a
	non-nearest-neighbor swap.
	\label{fig:irrprop}}
\end{center}
\end{figure}

\begin{notn}
Given $\vecx_1, \ldots, \vecx_N$ and a permutation $\pi$, we may write
$\pi$ as an \emph{image map} with the $\vecx_i$'s along the top row and
their images along the bottom row:
\begin{align*}
	\pcmatrix{\vecx_1 & \ldots & \vecx_N \\
	\pi(\vecx_1) & \ldots & \pi(\vecx_N)}
\end{align*}
\end{notn}

\begin{lem}
The composition $G_{\veca,\vecc} \circ G_{\veca,\vecb}$ behaves as follows:
\begin{align*}
	\pi =
	\pcmatrix{
		\ldots &     \veca  &     \vecb &      \vecc & \ldots\\
		\ldots & \pi(\veca) & \pi(\vecb) & \pi(\vecc)& \ldots
	}
	\mapsto
	\pi' =
	\pcmatrix{
		\ldots &     \veca  &     \vecb &      \vecc & \ldots\\
		\ldots & \pi(\vecc) & \pi(\veca) & \pi(\vecb)& \ldots
	}
\end{align*}
\begin{proof}
The first map, $G_{\veca,\vecb}$, does
\begin{align*}
	\pcmatrix{
		\ldots &     \veca  &     \vecb &      \vecc & \ldots\\
		\ldots & \pi(\veca) & \pi(\vecb) & \pi(\vecc)& \ldots
	}
	\mapsto
	\pcmatrix{
		\ldots &     \veca  &     \vecb &      \vecc & \ldots\\
		\ldots & \pi(\vecb) & \pi(\veca) & \pi(\vecc)& \ldots
	};
\end{align*}
$G_{\veca,\vecc}$ sends this to
\begin{align*}
	\pcmatrix{
		\ldots &     \veca  &     \vecb &      \vecc & \ldots\\
		\ldots & \pi(\vecc) & \pi(\veca) & \pi(\vecb)& \ldots
	}.
\end{align*}
\end{proof}
\end{lem}

\begin{cor}
The composition
$
	G_{\vecy_0,\vecy_n} \circ
	G_{\vecy_0,\vecy_{n-1}} \circ \ldots \circ
	G_{\vecy_0,\vecy_2} \circ
	G_{\vecy_0,\vecy_1}
$, sending $\pi\mapsto\pi'$, performs the right cyclic shift on images of
$\vecy_0, \ldots, \vecy_n$ given by
\begin{align*}
	%\pi =
	\pcmatrix{
	    \vecy_0 &    \vecy_1 &\ldots&    \vecy_{n-1}&     \vecy_n & \\
	\pi(\vecy_0)&\pi(\vecy_1)&\ldots&\pi(\vecy_{n-1})&\pi(\vecy_n)&
	}
	\mapsto
	%\pi' =
	\pcmatrix{
	    \vecy_0 &    \vecy_1 &\ldots&    \vecy_{n-1}&      \vecy_n \\
	\pi(\vecy_n)&\pi(\vecy_0)&\ldots&\pi(\vecy_{n-2})&\pi(\vecy_{n-1})
	}.
\end{align*}
Likewise,
$
	G_{\vecy_n,\vecy_1} \circ
	G_{\vecy_n,\vecy_2} \circ \ldots \circ
	G_{\vecy_n,\vecy_{n-2}} \circ
	G_{\vecy_n,\vecy_{n-1}}
$
leaves the image of $\vecy_0$ unchanged and performs the left cyclic shift on
images of $\vecy_1, \ldots, \vecy_n$ given by
\begin{align*}
	\pcmatrix{
	    \vecy_0 &    \vecy_1 &\ldots&    \vecy_{n-1} &    \vecy_n  \\
	\pi(\vecy_0)&\pi(\vecy_1)&\ldots&\pi(\vecy_{n-1})&\pi(\vecy_n)
	}
	\mapsto
	\pcmatrix{
	    \vecy_0 &    \vecy_1 &\ldots&    \vecy_{n-1} &    \vecy_n \\
	\pi(\vecy_0)&\pi(\vecy_2)&\ldots&\pi(\vecy_{n})  &\pi(\vecy_{1})
	}.
\end{align*}
\begin{proof}
These follow from the lemma by induction on $n$.
\end{proof}
\end{cor}

Composing these two maps, we find that
\begin{align*}
	(
	G_{\vecy_n,\vecy_1} \circ
	G_{\vecy_n,\vecy_2} \circ \ldots \circ
	G_{\vecy_n,\vecy_{n-2}} \circ
	G_{\vecy_n,\vecy_{n-1}}
	)\circ(
	G_{\vecy_0,\vecy_n} \circ
	G_{\vecy_0,\vecy_{n-1}} \circ \ldots \circ
	G_{\vecy_0,\vecy_2} \circ
	G_{\vecy_0,\vecy_1}
	)
\end{align*}
swaps the images of $\vecx=\vecy_0$ and $\vecz=\vecy_n$ while leaving
all other images unchanged, that is,
\begin{align*}
	\pcmatrix{
	    \vecy_0  &     \vecy_1  & \ldots&     \vecy_{n-1}  &     \vecy_n \\
	\pi(\vecy_0) & \pi(\vecy_1) & \ldots& \pi(\vecy_{n-1}) & \pi(\vecy_n)
	}
	\mapsto
	\pcmatrix{
	    \vecy_0  &     \vecy_1  & \ldots&     \vecy_{n-1}  &     \vecy_n \\
	\pi(\vecy_n) & \pi(\vecy_1) & \ldots& \pi(\vecy_{n-1}) & \pi(\vecy_0)
	}.
\end{align*}
\end{proof}
\end{prop}

\begin{rem*}
Below we will discuss winding cycles, and the empirical fact that the swap-only
algorithm reaches them only rarely.  The chain is irreducible but the non-zero
transition probability can still be very small.
\end{rem*}

\begin{prop}[Aperiodicity]
\label{prop:SO_aperiodic}
The swap-only algorithm's Markov chain is \plainidx{aperiodic}.
\begin{proof}
From Markov-chain theory \cite{bib:lawler}, we know that all states in a
recurrence class have the same period (see definition \ref{defn:period} on
page \pageref{defn:period}). Thus, it suffices to show that one state, e.g. the
identity, is aperiodic.  If there is non-zero probability of transitioning from
the identity to itself, the identity is an aperiodic state.  But this in fact
the case:  The identity has zero energy (equation \ref{eqn:RCM_H_rell_V}).  Any
state proposed as a transition from the identity is necessarily a two-cycle,
for which the distance terms in equation \ref{eqn:RCM_H_rell_V} are positive.
The cycle-weight terms or interaction terms are non-negative.  Thus $\Delta H >
0$, and the swap-only Metropolis transition probability is strictly less than
1.  Therefore there is a non-zero probability of staying at the identity
permutation.
\end{proof}
\end{prop}

\begin{prop}[Detailed balance]
\label{prop:SO_DB}
For all $\pi,\pi' \in \SN$,
\begin{align}
\label{eqn:SO_DB}
	P(\pi) M(\pi, \pi') &= P(\pi') M(\pi', \pi).
\end{align}
\begin{proof}
The detailed-balance statement in terms of the Gibbs distribution (equation
\ref{eqn:RCM_P}) and the swap-only Metropolis transition matrix (equation
\ref{eqn:SO_M}) is
\begin{align*}
	\frac{e^{-H(\pi)}}{Z}
	\left(1 \wedge e^{-H(\pi')} e^{H(\pi)}\right)
	&\maybeeq
	\frac{e^{-H(\pi')}}{Z}
	\left(1 \wedge e^{-H(\pi)} e^{H(\pi')}\right).
\end{align*}
The $Z$'s cancel.
The lemma below shows that $M(\pi,\pi') \ne 0$ iff $M(\pi',\pi) \ne 0$.
If $M(\pi,\pi') = 0$, then detailed balance holds.
If $M(\pi,\pi') \ne 0$, then there are two cases.
If $H(\pi') \le H(\pi)$, then
\begin{align*}
 	e^{-H(\pi)} \left(1\right)
	&=
	e^{-H(\pi')} \left(e^{-H(\pi)} e^{H(\pi')}\right).
\end{align*}
If $H(\pi') > H(\pi)$,
\begin{align*}
	e^{-H(\pi)} \left(e^{-H(\pi')} e^{H(\pi)}\right)
	&=
 	e^{-H(\pi')} \left(1\right).
\end{align*}
In all cases, detailed balance holds.
\end{proof}
\end{prop}

\begin{lem}
\label{lem:nzdb1}
For all $\pi,\pi' \in \SN$,
\begin{align*}
	M(\pi,\pi') \ne 0 \iff M(\pi', \pi) \ne 0.
\end{align*}
\begin{proof}
If $d(\pi,\pi') > 2$, then $M(\pi,\pi') = 0$ and $M(\pi',\pi) = 0$.

If $d(\pi,\pi') = 2$, then $\pi$ and $\pi'$ differ at two sites
$\vecx$ and $\vecy$.  Put
\begin{align*}
	\pi (\vecx) &= \vecu, &
	\pi'(\vecx) &= \vecv, \\
	\pi (\vecy) &= \vecv, &
	\pi'(\vecy) &= \vecu.
\end{align*}
If $\|\vecu-\vecv\| = 1$ then $M(\pi,\pi') \ne 0$ and $M(\pi',\pi) \ne 0$
since $\vecu$ and $\vecv$ are nearest-neighbor lattice sites.
If $\|\vecu-\vecv\| > 1$ then $M(\pi,\pi') = 0$ and $M(\pi',\pi) = 0$
since $\vecu$ and $\vecv$ are not nearest-neighbor lattice sites.

Last, if $\pi=\pi'$, then $M(\pi,\pi') = M(\pi',\pi)$.
\end{proof}
\end{lem}

This lemma completes the proof that the swap-only algorithm satisfies detailed
balance and thus has the Gibbs distribution as its invariant distribution.

The following proposition is not a correctness result, but rather a sanity
check.  It shows that cycles may grow or shrink upon swaps.

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.7]{figures/gkcycles.eps}
\caption[Swaps merge disjoint cycles and split single cycles.]
	{Swaps merge disjoint cycles and split single cycles.
	The left-hand permutation can be reached from the right-hand permutation
	via a swap, and vice versa.
	\label{fig:socycles}}
\end{center}
\end{figure}

\begin{prop}
\label{prop:SO_merge_or_split}
If $\vecx$ and $\vecy$ are\index{split}\index{merge}\index{cycle
split}\index{cycle merge} in disjoint cycles before a non-trivial swap then
they are in the same cycle afterward and vice versa (see figure
\ref{fig:socycles}).
\begin{proof}
First suppose that $\vecx$ and $\vecy$ are in disjoint cycles.  Let the
respective cycle lengths be $\ell(\vecx) = a$ and $\ell(\vecy)=b$.
Recall that $1 \le a,b \le N$.
Those cycles are
$$
	\vecx        \mapsto
	\pi(\vecx)   \mapsto
	\pi^2(\vecx) \mapsto
	\ldots       \mapsto
	\pi^{a-1}(\vecx) \mapsto
	\vecx
$$
and
$$
	\vecy        \mapsto
	\pi(\vecy)   \mapsto
	\pi^2(\vecy) \mapsto
	\ldots       \mapsto
	\pi^{b-1}(\vecy) \mapsto
	\vecy.
$$
Since these are disjoint cycles, all elements listed are distinct
lattice sites.  After the swap, we have
$$
	\vecy        \mapsto
	\pi(\vecx)   \mapsto
	\pi^2(\vecx) \mapsto
	\ldots       \mapsto
	\pi^{a-1}(\vecx) \mapsto
	\vecx
$$
and
$$
	\vecx        \mapsto
	\pi(\vecy)   \mapsto
	\pi^2(\vecy) \mapsto
	\ldots       \mapsto
	\pi^{b-1}(\vecy) \mapsto
	\vecy.
$$
This is a single cycle of length $a+b$, starting with $\vecy$, including
$\vecx$, and returning to $\vecy$.

Second,  suppose that $\vecx$ and $\vecy$ are in the same cycle.  Let
$a=\ellxypi$ and $b=\ellxypi$ (see definition
\ref{def:ellxy}).  (These numbers are both positive since the swap is
non-trivial, i.e.  $\vecx \ne \vecy$.) Then we have
$$
	\vecx        \mapsto
	\pi(\vecx)   \mapsto
	\pi^2(\vecx) \mapsto
	\pi^{a-1}(\vecx) \mapsto
	\vecy        \mapsto
	\pi(\vecy)   \mapsto
	\pi^2(\vecy) \mapsto
	\ldots       \mapsto
	\pi^{b-1}(\vecy) \mapsto
	\vecx.
$$
This is a single cycle of length $a+b$; all lattice sites listed are
distinct.  After the swap, we have
$$
	\vecy        \mapsto
	\pi(\vecx)   \mapsto
	\pi^2(\vecx) \mapsto
	\pi^{a-1}(\vecx) \mapsto
	\vecy
\qand
	\vecx        \mapsto
	\pi(\vecy)   \mapsto
	\pi^2(\vecy) \mapsto
	\ldots       \mapsto
	\pi^{b-1}(\vecy) \mapsto
	\vecx.
$$
These are disjoint cycles of length $a$ and $b$, respectively; the first
contains $\vecx$ and the second contains $\vecy$.
\end{proof}
\end{prop}

\begin{rem}
\label{rem:SO_merge_or_split_cyclens}
If $\vecx \ctdto \vecy$ (i.e. the swap splits\index{split} their common
cycle): the old cycle lengths $\ellxpi$ and $\ellypi$ are
equal, and the new cycle lengths after the swap are
\begin{align*}
	\ellxpiprime &= \ellyxpi \qqand
	\ellypiprime = \ellxypi.
\end{align*}
Otherwise (i.e. the swap merges\index{merge} the disjoint cycles):
\begin{align*}
	\ellxpiprime &= \ellypiprime = \ellxpi + \ellypi.
\end{align*}
\end{rem}

%% ----------------------------------------------------------------
\section{Winding numbers and the swap-and-reverse algorithm}
\label{sec:wno_SAR}

The propositions of section \ref{sec:SO_correctness} showed that the swap-only
algorithm is correct, asymptotically in the number of Metropolis sweeps --- in
particular, any permutation is reachable from any other with non-zero
probability.  However, in practice some of these transition probabilities can
be quite small.  In particular, we observe empirically that the swap-only
algorithm always generates permutations with zero winding number.
This is readily proved.

\begin{prop}
The swap step of the swap-only algorithm preserves winding number.
\begin{proof}
The permutations before and after the swap have winding numbers
\begin{align*}
	\vecW(\pi) &= \colvecthree{W_x(\pi)}{W_y(\pi)}{W_z(\pi)}
		= \sum_{i=1}^N \tilde\vecd(\pi(\vecx_i),\vecx_i),
	&
	\vecW(\pi') &= \colvecthree{W_x(\pi')}{W_y(\pi')}{W_z(\pi')}
		= \sum_{i=1}^N \tilde\vecd(\pi'(\vecx_i),\vecx_i)
\end{align*}
where the difference vector $\tilde\vecd$ is as in definition \ref{def:wno} and
remark \ref{rem:differences_on_torus}.  Since we work in the regime of short
jumps, and since $\pi'(\vecx_i)$ is a nearest neighbor of $\pi(\vecx_i)$, the
Euclidean charts overlap and the $\vecx_i$'s cancel when we subtract
$\vecW(\pi)$ from $\vecW(\pi')$.  Also, $\pi$ agrees with $\pi'$ except at the
two swap points $\vecx$ and $\vecy$; we have $\pi'(\vecx)=\pi(\vecy)$ and vice
versa.  Thus the change in winding number is computed entirely in the same
chart and we have
\begin{align*}
	\vecW'-\vecW
	&= \frac{1}{L} \sum_{i=1}^N \tilde\vecd(\pi'(\vecx_i), \pi(\vecx_i)) \\
	&= \frac{1}{L} \left[ \tilde\vecd(\pi'(\vecx), \pi(\vecx))
	+ \tilde\vecd(\pi'(\vecy), \pi(\vecy)) \right] \\
	&= \frac{1}{L} \left[ \tilde\vecd(\pi(\vecy), \pi(\vecx))
	+ \tilde\vecd(\pi(\vecx), \pi(\vecy)) \right]
	= 0.
\end{align*}
\end{proof}
\end{prop}

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics{figures/wno_zero_pf.eps}
\caption
	[Illustration of winding-number conservation.]
	{Example permutations $\pi$ (left) and $\pi'$ (right)
	illustrating winding-number conservation.
	\label{fig:wno_zero_pf}}
\end{center}
\end{figure}

An example is shown in figure \ref{fig:wno_zero_pf} in dimension $d=2$
with $N=8$ points on a lattice of width $L=4$.  The sites
affected by the permutation swap are $\vecx=\vecx_3$ and $\vecy=\vecx_7$.  The
change in winding numbers is
\begin{align*}
	\vecW(\pi') - \vecW(\pi)
	&= \frac{1}{4}\left[ \tilde\vecd(\pi'(\vecx_3),\vecx_3)
	+  \tilde\vecd(\pi'(\vecx_7),\vecx_7)
	-  \tilde\vecd(\pi (\vecx_3),\vecx_3)
	-  \tilde\vecd(\pi (\vecx_7),\vecx_7) \right] \\
	&= \frac{1}{4}\left[
			\colvectwo{1}{0} +
			\colvectwo{-1}{0}
		\right]
	- \frac{1}{4}\left[
			\colvectwo{0}{-1} +
			\colvectwo{0}{1}
		\right]
		= \colvectwo{0}{0} \\
	&= \frac{1}{4}\left[ \tilde\vecd(\pi'(\vecx),\pi(\vecx)) \right]
	+  \frac{1}{4}\left[ \tilde\vecd(\pi'(\vecy),\pi(\vecx)) \right] \\
	&= \frac{1}{4}\left[
			\colvectwo{1}{1} +
			\colvectwo{-1}{-1}
		\right].
\end{align*}

A partial solution is explained intuitively by figure \ref{fig:sowinding}.
Part 1 of the figure shows a permutation $\pi$ with a long cycle on the torus
which almost meets itself in the $x$ direction.  In part 2, after a Metropolis
step sending $\pi$ to $\pi'$, one cycle winds by $+1$ and the other by $-1$.
Metropolis steps create winding cycles only in opposite-direction pairs; the
total $W_x(\pi)$ is still zero.  Part 3 of the figure shows that if we reverse
one cycle (which is a zero-energy move), $W_x(\pi)$ is now 2.  In general,
winding numbers of even parity can be generated.

\begin{figure}
\begin{center}
\psfragscanon
\includegraphics[scale=0.8]{figures/gkwinding.eps}
\caption[Conservation of winding number in the SO algorithm.]
	{Conservation of winding number in the swap-only algorithm, and a
	partial solution provided by the swap-and-reverse algorithm.
	\label{fig:sowinding}}
\end{center}
\end{figure}

The \emphidx{swap-and-reverse} algorithm adds a second type of sweep to the
swap-only algorithm.  Namely:  (1) In a swap-only sweep, For each lattice site,
one does a Metropolis step as above.  (2) In a cycle-reversing sweep, for each
cycle in the permutation, one reverses the direction of the cycle with
probability $1/2$.  This permits winding numbers of even parity in each of the
three axes.  The correctness proof is unaffected, since cycle reversal is a
zero-energy change.  The time required to reach permutations with non-zero
winding numbers, which the asymptotic correctness proof does not address, is
reduced.  The additional penalty in terms of CPU time consumed by cycle
reversal is found to be negligible.

Other solutions exist to generate winding numbers of arbitrary parity: the
\plainidx{band-update} method of chapter \ref{chap:band_updates} and the worm
algorithm of chapter \ref{chap:worm}.  As will be shown, they suffer from
too-low acceptance rate and too-long stopping time, respectively.  Therefore,
the swap-and-reverse algorithm is our current best algorithm; it is used to
generate all the results discussed in chapter \ref{chap:results}.  We believe
that sampling permutations with even winding number does not harm our final
result.  The order parameters $f_S$ and $f_W$ depend on winding phenomena, but
the other three, $1/\xi$, $f_I$, and $f_M$, do not; furthermore, results
obtained in chapter \ref{chap:results} using each of the five order parameters
are compatible.

%% ================================================================
\chapter{Band updates} \label{chap:band_updates}

Winding cycles are a global phenomenon, while the SO swaps of section
\ref{sec:SO_algorithm} are local.  The Swendsen-Wang algorithm for the planar
Ising model addresses a global/local problem not unlike our winding-number
problem, which was presented in detail in section \ref{sec:wno_SAR}.  Motivated
by Swendsen-Wang, we attempt to define a Metropolis step which changes one of
the winding-number components $W_x$, $W_y$, or $W_z$ by $\pm 1$.  Our attempt
at non-local updates has an unreasonably low acceptance rate, namely, on the
order of $e^{-L}$ where $L$ is the box length.  Nonetheless, this concept may
provide fodder for better, future ideas.

%% ----------------------------------------------------------------
\section{The algorithm}

As discussed in section \ref{sec:SO_correctness}, a Metropolis transition from
$\pi$ to $\pi'$ may be thought of as composition with another permutation
$\tau$, i.e. $\pi' = \tau \pi$.  We replace the swap operator
$\tau=G_{\veca,\vecb}$ of section \ref{sec:SO_correctness}, which was a
two-cycle, with a permutation $\tau$ which has a winding $L$-cycle.  Without
loss of generality, it suffices to discuss a permutation $\tau$ which sends
every lattice point to itself except for one line of points with $y$ and $z$
coordinates equal to zero (figure \ref{fig:band_updates_tau}).  This $\tau$
will have winding number $(W_x,W_y,W_z)=(+1,0,0)$.  If we can do that, then by
reflection and rotation symmetries we can construct similar $\tau$'s with
$W_x$, $W_y$, or $W_z$ equal to $\pm 1$.

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.7]{figures/band_updates.eps}
\caption
	[Band update as composition with an $L$-cycle.]
	{Example of band update as composition with an $L$-cycle: $\pi' = \tau\pi$.
	The winding number $W_x$ increases by $+1$.
	\label{fig:band_updates_tau}}
\end{center}
\end{figure}

Figure \ref{fig:band_updates_tau} displays the idea.  The permutation $\tau$ is
an $L$-cycle; we put $\pi' = \tau\pi$.  It seems clear from the picture that
the winding number is modified in the desired manner, and moreover it seems
plausible to conjecture that winding numbers are additive with respect to
permutation composition, i.e. that $\vecW(\pi') = \vecW(\tau) + \vecW(\pi)$.
Proving either of these statements rigorously would be worthwhile if band
updates were worth pursuing.  However, as shown in the next section, they are
not.

%% ----------------------------------------------------------------
\section{Acceptance rate}

We examine the acceptance rate of the band-update algorithm using a
semi-empirical method.  Recall from section \ref{sec:theory_mjl} that we define
the random variable $j_\vecx(\pi)$ to be $\|\pi(\vecx)-\vecx\|$.  By
examination of histograms acquired over MCMC simulations, we see that this
parameter is non-critical; its values change smoothly with $T$ near $T_c$.  For
$T=6$ and $7$, respectively, we find the following distributions for $j_\vecx$:
$$
\begin{tabular}{|c|c|c|c|}
\hline $j$ & $P(J = j), T=6$ & $j$ & $P(J=j), T=7$ \\
\hline
\hline 0.0000000 & 0.4801 & 0.0000000 & 0.7745 \\
\hline 1.0000000 & 0.3361 & 1.0000000 & 0.1822 \\
\hline 1.4142136 & 0.1545 & 1.4142136 & 0.0378 \\
\hline 1.7320508 & 0.0216 & 1.7320508 & 0.0042 \\
\hline 2.0000000 & 0.0036 & 2.0000000 & 0.0012 \\
\hline 2.2360680 & 0.0032 & 2.2360680 & 0.0001 \\
\hline 2.4494897 & 0.0009 &           &        \\
\hline
\end{tabular}
$$
Thus for $T=6,7$, respectively, mean jump lengths are $0.608$ and $0.245$,
while mean squared jump lengths are $0.746$ and $0.276$.  Roughly, the latter
is of order $0.5$.  Recall in particular that it is the squared jump
length which provides the main contribution to the system energy (equation
\ref{eqn:RCM_H_rell_V}).

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.7]{figures/band_updates_2.eps}
\caption
	[Change in jump lengths on a band update.]
	{Change in jump lengths for a point affected by a band update.
	\label{fig:band_updates_delta}}
\end{center}
\end{figure}

The left-hand side of figure \ref{fig:band_updates_delta} shows (in two
dimensions only) various possibilities for a permutation arrow from a point
$\vecx$.  The right-hand side of the figure shows what happens to the jump
lengths at $\vecx$ on a right shift.  (All $L$ points in the $L$ cycle will be
affected similarly.) An up arrow of length 1 becomes a diagonal up-left arrow
of length $\sqrt 2$, a right arrow of length 1 becomes a right arrow of length
$2$, and so on.  Given an attempted band update $\pi \to \pi'$, the typical
value of the energy change $\Delta H$ at a single point (remembering that
cycle-weight corrections are low-order perturbations, since $\alpha$ is small)
is
\begin{align*}
	\frac{T}{4} \left( \bbE[j_\vecx(\pi')^2] - \bbE[j_\vecx(\pi)^2] \right).
\end{align*}
From the table just above, we can compute the expected values of
$j_\vecx(\pi)^2$ and $j_\vecx(\pi')^2$.  We find the energy difference per
point to be on the order of $+0.8 \cdot T/4 \approx +1.2$.  Since $L$ points
are involved in a band update, this means $\bbE[\Delta H] \approx 1.2 L \approx
L$.  Transitions are accepted with probability $\min\{1,e^{-\Delta H}\}$
(section \ref{sec:MCMC_recipe}) which is approximately $e^{-L}$.  We need to
consider $L$ from approximately $30$ and upward to get past the most severe
finite-size effects; $e^{-30} \approx 10^{-14}$ is effectively zero, and thus
proposed band updates are effectively never accepted.  This calculation matches
with simulation tests performed in software.

%% ================================================================
\chapter{The worm algorithm}
\label{chap:worm}

This chapter, like chapter \ref{chap:worm}, presents a Markov chain for
MCMC sampling of the model of random spatial permutations (chapter
\ref{chap:model_defns}) within the framework of chapter \ref{chap:MCMC}.
The worm algorithm solves the problem of winding-number-parity conservation
within the SAR algorithm (section \ref{sec:wno_SAR}).  However, a stopping
time which is poorly bounded in the lattice size prevents the worm algorithm
from being our algorithm of choice.

%% ----------------------------------------------------------------
\section{The idea of the worm algorithm}

Worm algorithms have been used heavily in path-integral Monte Carlo (PIMC)
simulations:  see \cite{bib:bps06, bib:pst98, bib:gcl97, bib:kps, bib:nl04,
bib:pc87}.  The context is that interparticle interactions are modeled using
Brownian bridges in the Feynman-Kac approach.  A naive, pre-PIMC sampling
approach involves generating separate Brownian bridges from point A to point B.
The PIMC idea is to generate a single Brownian bridge, then modify a bit at a
time using MCMC methods.

For the random-cycle model with true Bose interactions, the Brownian bridges
implicit in the $V$ terms of equation \ref{eqn:RCM_H_gen_V} have been
integrated out in equation \ref{eqn:RCM_H_rell_V}.  In fact, this is the key
selling point of the random-cycle model in the larger context of the Bose gas
(larger meaning beyond the scope of this dissertation).  Most of the complexity
of PIMC simulations, which are an efficient approach to Brownian bridges, goes
away.  If one were to adapt a PIMC worm algorithm to the RCM, one would need to
spend significant time learning about PIMC.  Yet it is likely that most of the
complexity will also go away.  Instead, it is simpler to ask:  If one were to
have a worm algorithm for the random-cycle model, what properties would it
have?  We require the following:
\begin{itemize}
\item We have a lattice with a fixed number $N$ of points.  There is no
desire to work in the grand-canonical ensemble.
\item We want the ability to open and close permutation cycles.
(An open cycle is a ``worm''.)
\item Given that, tips of open cycles may wander around the 3-torus before
closing, permitting arbitrary winding numbers.
\end{itemize}
Thus, we want to sometimes open a cycle, then modify it with SO-like steps,
then close it again.  Following PIMC worm algorithms, all our Metropolis steps
will involve the worm.  This does touch all lattice points:  a worm is opened
at a site, then modified, then closed.  Then, a worm is opened somewhere else,
and so on.

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.6]{figures/worm1.eps}
\caption{Open cycles as permutations on $N+1$ points.
	\label{fig:worm1}}
\end{center}
\end{figure}

Question:  Can we leverage our knowledge of permutations?  To see how, consider
a closed cycle and an open cycle on $N=3$ points (see also figure
\ref{fig:worm1}):
\begin{equation*}
	\pcmatrix{
		1 & 2 & 3 \\
		2 & 3 & 1 \\
	},
	\qquad
	\pcmatrix{
		1 & 2 & 3 & {}\\
		2 & 3 & {} & 1 \\
	}
\end{equation*}
In the open cycle, $1\mapsto 2$, $2\mapsto 3$, $3 \mapsto$ nothing, and nothing
$\to 1$.  Call that nothing \emph{something} --- the \emphidx{wormhole point}.
It is an $(N+1)$st point, $w$:
$$
	\pcmatrix{
		1 & 2 & 3 & w \\
		2 & 3 & w & 1 \\
	}
$$
Henceforth, the wormhole point will be written as $w$ or $N+1$.  In diagrams,
it will be an open dot while the other $N$ points will be written with filled
dots (as in figure \ref{fig:worm1}).
Now we have permutations on $\SNp$.  Given $\pi\in\SN$, \plainidx{inject} $\pi$
into $\SNp$ via $\pi(w)=w$.

\begin{defn}
For $\pi\in\SNp$, we say $\pi$ is a \emphidx{closed permutation}
if $\pi(w)=w$.
We say $\pi$ is an \emphidx{open permutation}
if $\pi(w) \ne w$.
Likewise, a cycle of $\pi$ is said to be open or closed, respectively, if it
does or does not contain $w$.
\end{defn}

\begin{rem*}
The PIMC jargon is that closed permutations are in the \emphidx{$Z$ sector}
(for partition function), while open permutations are in the \emphidx{$G$
sector} (for Matsubara Green's function).
\end{rem*}

The goal is to invent an energy function, Gibbs distribution, and Metropolis
algorithm for these extended permutations in $\SNp$ such that the marginal
distribution on $\SNp$, conditioned on closed permutations, matches the RCM
Gibbs distribution (equation \ref{eqn:RCM_P}).  Then, random variables
will be sampled only at closed permutations.

%% ----------------------------------------------------------------
\section{Extended random-cycle model}

Recall\index{ERCM} that we inject $\pi\in\SN$ into $\SNp$ via $\pi(w)=w$.  The
$(N+1)$st point $w$ is \emphidx{non-spatial}:  it has no distance associated
with it.

\begin{defn}
The \emphidx{extended lattice} is\index{$\Lambda'$}
$$
	\Lambda' := \Lambda \cup \{w\}
$$
\end{defn}

\begin{defn}
For $\pi\in\SNp$, define
\begin{align}
\label{eqn:ERCM_H_rell_V}
	H'(\pi) &= \frac{1}{4\beta} \sum_{\stackrel{i=1}{\pi(\vecxi) \ne w}}^N
		\|\vecxi - \vecxpii\|^2 + \sum_{\ell=2}^N \alpha_\ell r_\ell(\pi)
		+ \gamma 1_{\SNpc}(\pi).
\end{align}
\end{defn}
Note that this \emphidx{extended energy} agrees with the RCM energy (equation
\ref{eqn:RCM_H_rell_V}) on closed permutations.  This is used to prove the
marginality condition below.  A small positive $\gamma$ factor is sufficient to
establish aperiodicity of the Markov chain (proposition
\ref{prop:worm_aperiodic}).  The \emphidx{extended Gibbs distribution} and
\emphidx{extended partition function} are defined in the obvious way.
\begin{defn}
Let
\begin{align}
\label{eqn:ERCM_Gibbs_P}
	P'(\pi) = \frac{e^{-H'(\pi)}}{Z'}
\end{align}
where the \plainidx{partition function} is
\begin{align}
\label{eqn:ERCM_Gibbs_Z}
	Z' &= \sum_{\pi\in\SNp} e^{-H'(\pi)}.
\end{align}
\end{defn}

%% ----------------------------------------------------------------
\section{Proof of marginality}

As long as the energy function for the ERCM and the RCM agree on closed
permutations, the desired marginality condition holds.  This means that the
interaction in section \ref{sec:prb_model} --- or any other to-be-invented
interaction model --- may use the worm algorithm as long as it agrees on closed
permutations.

\begin{prop}[Marginality condition]
\label{prop:marginality}
Let $\SN \into \SNp$ by taking $\pi(w)=w$.
Let $H,H'$ be energy functions on $\SN$ and $\SNp$, respectively,
such that for all $\pi\in\SN$,
\begin{align}
	H(\pi) = H'(\pi).
\end{align}
Let $P,P',Z,Z'$ be as above.  Then for $\pi\in\SN$,
\begin{align}
\label{eqn:marg_P}
	P'(\pi \mid \pi\in\SN) = P(\pi).
\end{align}
\begin{proof}
Let $\pi\in\SN$.  The left-hand side of equation \ref{eqn:marg_P} is, by
definition of conditional expectation,
\begin{align*}
	P'(\pi \mid \pi\in\SN)
	&= \frac{P'(\pi) \;1_{\SN}(\pi)}{P'(\SN)}.
\end{align*}
The numerator is the Gibbs probability for closed permutations, or zero for
open ones:
\begin{align*}
	P'(\pi) \;1_{\SN}(\pi)
	&= \frac{1}{Z'}e^{-H'(\pi)} \;1_{\SN}(\pi)
	= \frac{1}{Z'}e^{-H (\pi)} \;1_{\SN}(\pi)
\end{align*}
since $H$ and $H'$ agree on closed permutations.  The denominator is the total
probability of closed permutations:
\begin{align*}
	P'(\SN)
	&= \frac{1}{Z'} \sum_{\pi \in\SN} e^{-H'(\pi)}
	= \frac{1}{Z'} \sum_{\pi \in\SN} e^{-H(\pi)}.
\end{align*}
Since $\pi\in\SN$, the ratio is
\begin{align*}
	\frac{
		\frac{1}{Z'}e^{-H (\pi)} \;1_{\SN}(\pi)
	}{
		\frac{1}{Z'} \sum_{\pi \in\SN} e^{-H(\pi)}
	}
	&= \frac{
		e^{-H (\pi)} \;1_{\SN}(\pi)
	}{
		\sum_{\pi \in \SN} e^{-H(\pi)}
	} \\
	%
	&= \frac{
		e^{-H (\pi)} \;1_{\SN}(\pi)
	}{
		Z
	} \\
	%
	&= P(\pi).
\end{align*}
\end{proof}
\end{prop}

%As a sanity check, we see whether the marginal density for the $\beta=0$
%case recovers uniform permutations on $\SN$, and whether the marginal density
%for the $\beta=\infty$ case recovers the singleton density on $\SN$ supported
%at the identity permutation.
%
%[Needs to be typed up.]

%% ----------------------------------------------------------------
\section{The worm algorithm}
\label{sec:worm_algo}

Now that we have the correct Gibbs distribution for the ERCM, the next step
is to devise a Metropolis algorithm to sample from it.  Below,
we will prove correctness.

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.7]{figures/wormmoves.eps}
\caption{Metropolis moves for the worm algorithm.
	\label{fig:wormmoves}}
\end{center}
\end{figure}

The \emphidx{worm algorithm} is as follows:
\begin{itemize}

\item Start with the identity (or uniform-random closed) permutation.

\item The permutation is now closed, so $\pi(w)=w$.  Select a lattice
site $\vecx$ at uniform random.  With probability proportional to $1 \wedge
\edh$, open the permutation by swapping the arrows of $\vecx$ and $w$.
This is called an \emphidx{open} move.
(See figure \ref{fig:wormmoves}.)

%	\begin{align*}
%	\pi: &
%	\pcmatrix{
%		\ldots & \vecx      & \ldots &     w      & \ldots \\
%		       & \downarrow &        & \downarrow &        \\
%		\ldots & \pi(\vecx) & \ldots & \pi(w)=w   & \ldots \\
%	} &
%	\pi': &
%	\pcmatrix{
%		\ldots & \vecx      & \ldots & \vecy      & \ldots \\
%		       & \downarrow &        & \downarrow &        \\
%		\ldots &     w & \ldots & \pi(\vecx) & \ldots \\
%	} \\
%	\end{align*}

\item Now that the permutation is open, do a \emphidx{head swap}, \emphidx{tail
swap}, or \emphidx{close}.

\item Head swap:  Pick a lattice site $\vecx$ nearest-neighbor to the lattice
site $\piiw$.  With probability proportional to $1 \wedge \edh$, swap arrows as
in figure \ref{fig:wormmoves}.  The head swap is trivial if $\vecx = \piiw$,
which happens only if the head swap is rejected.  The head swap would be a
close if $\vecx = w$, but we choose $\vecx$ to be a lattice site.  Thus, the
permutation remains open on a head swap.

% * A trivial head swap is a stay:  if x == piiw.
% * A head swap is a close if x == w.

\item Tail swap:  Pick a lattice site $\pix$ nearest-neighbor to the lattice
site $\piw$.  With probability proportional to $1 \wedge \edh$, swap arrows as
in figure \ref{fig:wormmoves}.  The tail swap is trivial if $\pix = \piw$,
which happens only if the tail swap is rejected.  The tail swap would be a
close if $\pix = w$, but we choose $\pix$ to be a lattice site.  Thus, the
permutation remains open on a tail swap.

% * A trivial tail swap is a say:  if x == w.
% * A tail swap is a close if pi(x) == w.

\item Close: with probability proportional to $1 \wedge \edh$, swap arrows as
in figure \ref{fig:wormmoves}.  The permutation is now closed.

\item Once the permutation is closed --- after an open, some number of
head/tail swaps, and a close, or after a rejected open --- a \emphidx{worm
sweep} has been completed.  At every sweep (or after every specified number of
sweeps as noted in the autocorrelation footnote above), one may obtain a value
of random variable(s) for inclusion in computation of their sample means.

\end{itemize}

\begin{defn}
\label{defn:worm_trivial}
A head swap at $\vecx$ is \emphidx{trivial} if $\vecx = \piiw$;
a tail swap at $\vecx$ is trivial if $\pix = \piw$.
\end{defn}

%% ----------------------------------------------------------------
\section{Fibration of $\SNp$ over $\SN$}

The definitions and lemmas in this section facilitate explicit construction
of the Markov matrix, and are necessary for proving correctness of the worm
algorithm.

The key points about the structure of the fibration, formalized by the lemmas
below, are:
\begin{itemize}
\item Each open permutation is one opener move away from a base closed
permutation.  The $N$ open permutations above a base closed permutation $\pi$
are the \plainidx{fiber} over $\pi$.
\item This induces a \plainidx{disjoint partition} of the open permutations
$\SNpc$.
\item Opens and closes stay within fibers; non-trivial head swaps and tail
swaps \plainidx{cross fibers}.
\item For each open permutation, the six non-trivial head swaps and six tail
swaps result in twelve distinct permutations.
\item Head swaps and tail swaps are transitive on fibers.
\end{itemize}

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics[scale=0.9]{figures/fiber4.eps}
\caption[Fibration of $\mcS_4$ over $\mcS_3$.]
	{Fibration of $\mcS_4$ over $\mcS_3$.  Closed permutations (i.e. $\mcS_3$)
	are along the bottom row; open permutations (i.e. $\mcS_4 \setminus
	\mcS_3$) are above the bottom row.  The column, or fiber, above each closed
	permutation $\pi$ contains the open permutations obtained from $\pi$ by an
	opener move.  Arrows modified by opener moves are shown in black.
	\label{fig:fiber4}}
\end{center}
\end{figure}

We first define maps corresponding to worm Metropolis moves.

\begin{defn} The four worm Metropolis moves of figure \ref{fig:wormmoves} may
be viewed in terms of maps.  Throughout, $\vecz \in \Lambda \cup \{w\}$.

Let $O: \;\SN\times\Lambda \to \SNpc$ send $O(\pi,\vecx) = \pi'$ such that
\begin{align*}
	\pi'(\vecx) &= w, \\
	\pi'(w)     &= \pix, \\
	\pi'(\vecz) &= \pi(\vecz), \quad \vecz \ne \vecx, w.
\end{align*}
Let $C: \;\SNpc \to \SN$ send $C(\pi) = \pi'$ such that
\begin{align*}
	\pi'(\piiw) &= \piw, \\
	\pi'(w)     &= w, \\
	\pi'(\vecz) &= \pi(\vecz), \quad \vecz \ne \piiw, w.
\end{align*}
Let $S: \;\SNpc\times\Lambda \to \SNpc$ send $S(\pi,\vecx) = \pi'$ such that
\begin{align*}
	\pi'(\vecx) &= w, \\
	\pi'(\piiw) &= \pix, \\
	\pi'(\vecz) &= \pi(\vecz), \quad \vecz \ne \vecx, \piiw.
\end{align*}
Let $T: \;\SNpc\times\Lambda \to \SNpc$ send $T(\pi,\vecx) = \pi'$ such that
\begin{align*}
	\pi'(\vecx) &= \piw, \\
	\pi'(w)     &= \pix, \\
	\pi'(\vecz) &= \pi(\vecz), \quad \vecz \ne \vecx, w.
\end{align*}

\end{defn}

Throughout the proofs of the fibration-structure lemmas, we will use the
following fact.
\begin{lem}
\label{lem:pineq}
If $\vecx \ne \vecy$, then $\pix \ne \piy$ and $\piix \ne \piiy$.
\begin{proof}
%By contrapositive:  if $\pix = \piy$ then $\pii(\pi(\vecx)) = \pii(\pi(\vecy))$
%and so $\vecx = \vecy$.
If $\vecx \ne \vecy$ and $\pix = \piy$, then $\pi$ is not 1-1 which is
a contradiction since $\pi$ is a permutation.
This applies to $\pii$ as well, since $\pii$ is also a
permutation.
\end{proof}
\end{lem}

Now we may prove the fibration-structure lemmas.

\begin{lem}
Each open permutation $\pi$ is one opener move away from a base closed
permutation $\pi'$.  That is,  for all $\pi \in \SNpc$, there exists $\pi' \in
\SN$ such that $C(\pi) = \pi'$.
\begin{proof}
Let $\pi\in\SNp$.  Since $\pi$ is open, $\piw \ne w$ and $\piiw \ne w$.
Let $\veca = \piiw$ and $\vecb = \piw$.  Both are lattice points.
Applying $C$, we have $C(\pi)=\pi'$ where $\pi'(\veca)=\vecb$, $\pi'(w)=w$,
and $\pi'(\vecz)=\pi(\vecz)$ for all remaining lattice points $\vecz \ne
\veca, \vecb$.  Since $\pi'(w)=w$, $\pi'$ is closed.
\end{proof}
\end{lem}

\begin{defn}
For $\pi\in\SN$, $C^{-1}(\pi) \subset \SNpc$ is the \textbf{fiber} of open
permutations over $\pi$.
\end{defn}

\begin{lem}
\label{lem:worm_opens_and_closes}
Opens and closes stay within fibers, and each fiber has $N$ elements.
\begin{proof}
Closes stay within fibers by definition of fiber.  Next, fix $\pi\in\SN$ and
let $\vecx_1, \vecx_2 \in \Lambda$.  (These are two different ways to open the
same closed permutation.) Let
$$
	\pi'_1 = O(\pi,\vecx_1), \quad
	\pi'_2 = O(\pi,\vecx_2).
$$
Then $\pi'_1$ and $\pi'_2$ have
$$
	\vecx_1 \mapsto w \mapsto \pi(\vecx_1), \qquad
	\vecx_2 \mapsto w \mapsto \pi(\vecx_2),
$$
respectively, agreeing with $\pi$ at all other lattice points $\vecz$.
Now, $C(\pi'_1)$ and $C(\pi'_2)$ have
$$
	\vecx_1 \mapsto \pi(\vecx_1), w \mapsto w, \qquad
	\vecx_2 \mapsto \pi(\vecx_2), w \mapsto w
$$
respectively, agreeing with $\pi$ at all other lattice points $\vecz$.  But
this means $C(\pi'_1)$ agrees with $C(\pi'_2)$ agree at all points of
$\Lambda'$, so $C(\pi'_1) = C(\pi'_2)$.  Thus, $\pi'_1$ and $\pi'_2$ are in the
same fiber.

For the last claim, fix $\pi\in\SN$ and enumerate the $N$ lattice points of
$\Lambda$ as $\vecx_1, \ldots, \vecx_N$.  We claim that the $N$ permutations
$$
	\pi'_1 = O(\pi,\vecx_1), \ldots, \pi'_N = O(\pi,\vecx_N),
$$
which are all now known to be in the same fiber, are all distinct.  To see
this, fix $i \ne j$ from out of $\{1, 2, \ldots, N\}$.  Then $\pi'_i$ and
$\pi'_j$ have
$$
	\vecx_i \mapsto w \mapsto \pi(\vecx_i), \qquad
	\vecx_j \mapsto w \mapsto \pi(\vecx_j).
$$
Since $\vecx_i \ne \vecx_j$, by lemma \ref{lem:pineq} $\pi(\vecx_i) \ne
\pi(\vecx_j)$.  Since
$$
	\pi'_i(w) = \pi(\vecx_i) \ne \pi(\vecx_j) = \pi'_j(w),
$$
$\pi'_i$ and $\pi'_j$ send $w$ to different points.  Therefore, the
permutations $\pi'_i$ and $\pi'_j$ are distinct.
\end{proof}
\end{lem}

\begin{lem}
This fibration induces a disjoint partition of the open permutations $\SNpc$.
That is, for $\pi'_1, \pi'_2 \in \SN$,
\begin{align*}
	\pi'_1 \ne \pi'_2 \implies
	C^{-1}(\pi'_1) \cap C^{-1}(\pi'_2) &= \emptyset
	\qqand
	\bigcup_{\pi\in\SN} C^{-1}(\pi) = \SNpc.
\end{align*}
\begin{proof}
For the first claim, suppose the intersection is non-empty.  Let $\pi \in\SNpc$
be such that $\pi \in C^{-1}(\pi'_1)$ and $\pi \in C^{-1}(\pi'_2)$.  This means
$C(\pi) = \pi'_1$ and $C(\pi) = \pi'_2$ with $\pi'_1 \ne \pi'_2$, which is a
contradiction since the map $C$ is uniquely defined for all $\pi \in \SNpc$.

For the second claim: there are $N!$ closed permutations.  We know from the
first claim that the $N!$ fibers, one above each closed permutation, are all
disjoint.  From lemma \ref{lem:worm_opens_and_closes}, we know that each fiber
has $N$ elements.  We have accounted for all $N\cdot N! = (N+1)!-N!$ open
permutations, so we must have all of $\SNpc$.
\end{proof}
\end{lem}

%           * * * * * * * * each of height N
%           * * * * * * * *
%           * * * * * * * *
% N! fibers * * * * * * * *
% N! closed . . . . . . . .

\begin{lem}
Non-trivial head swaps and tail swaps (definition \ref{defn:worm_trivial})
cross fibers.
\begin{proof}
First consider head swaps.
Let $\pi,\pi' \in \SNpc$ differ by a non-trivial head swap, namely, there is
$\vecx \ne \piiw$ such that $\pi' = S(\pi, \vecx)$.
Then $\pi$ and $\pi'$ have
$$
	\begin{array}{lccccc}
	\pi:  & \vecx &\mapsto& \pix &\mapsto& \pi^2(\vecx), \\
	\pi': & \piiw &\mapsto& \pix &\mapsto& \pi^2(\vecx), \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	\piiw &\mapsto&    w &\mapsto& \piw, \\
	\vecx &\mapsto&    w &\mapsto& \piw,
	\end{array}
$$
respectively.  Now apply $C$ to each: $C(\pi)$ and $C(\pi')$ have
$$
	\begin{array}{lccccc}
	C(\pi):  & \vecx &\mapsto& \pix &\mapsto& \pi^2(\vecx), \\
	C(\pi'): & \piiw &\mapsto& \pix &\mapsto& \pi^2(\vecx), \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	\piiw &\mapsto&    \piw, \\
	\vecx &\mapsto&    \piw,
	\end{array}
\qquad
	\begin{array}{lccccc}
	w &\mapsto&    w, \\
	w &\mapsto&    w,
	\end{array}
$$
respectively.  Since $\vecx \ne \piiw$, $C(\pi) \ne C(\pi')$.

Next, consider tail swaps.
Let $\pi,\pi' \in \SNpc$ differ by a non-trivial tail swap, namely, there is
$\pix \ne \piw$ such that $\pi' = T(\pi, \vecx)$.
Then $\pi$ and $\pi'$ have
$$
	\begin{array}{lccccc}
	\pi:  & \piix &\mapsto& \vecx &\mapsto& \pix, \\
	\pi': & \piix &\mapsto& \vecx &\mapsto& \piw, \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	\piiw &\mapsto&    w &\mapsto& \piw, \\
	\piiw &\mapsto&    w &\mapsto& \pix, \\
	\end{array}
$$
respectively.  Now apply $C$ to each: $C(\pi)$ and $C(\pi')$ have
$$
	\begin{array}{lccccc}
	C(\pi):  & \piix &\mapsto& \vecx &\mapsto& \pix, \\
	C(\pi'): & \piiw &\mapsto& \vecx &\mapsto& \piw, \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	\piiw &\mapsto&    \piw, \\
	\piiw &\mapsto&    \pix,
	\end{array}
\qquad
	\begin{array}{lccccc}
	w &\mapsto&    w, \\
	w &\mapsto&    w,
	\end{array}
$$
respectively.  Since $\pix \ne \piw$, $C(\pi) \ne C(\pi')$.
\end{proof}
\end{lem}

\begin{lem}
\label{lem:worm_distinct_dozen}
For each open permutation, the six non-trivial head swaps and six non-trivial
tail swaps result in twelve distinct permutations.
\begin{proof}
Fix $\pi\in\SNpc$.
Let $\vecx_1, \ldots, \vecx_6$ be the six nearest-neighbor lattice sites to the
lattice site $\piiw$;
let $\vecy_1, \ldots, \vecy_6$ be the six lattice sites
such that $\pi(\vecy_1), \ldots, \pi(\vecy_6)$ are nearest-neighbor lattices
site to the
lattice site $\piw$.  (See figure \ref{fig:wormmoves}.)

First, we show that the six permutations $S(\pi,\vecx_1), \ldots,
S(\pi,\vecx_6)$ are distinct.  Let $i \ne j$ for $i,j=1, \ldots, 6$; let $\pi_i
= S(\pi,\vecx_i)$ and $\pi_j = S(\pi,\vecx_j)$.  Then $\pi$ has
$$
	\begin{array}{lccccc}
	\pi:   & \vecx_i &\mapsto& \pi(\vecx_i), \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	\vecx_j &\mapsto& \pi(\vecx_j), \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	\piiw &\mapsto& w; \\
	\end{array}
$$
$\pi_i$, and $\pi_j$ have
$$
	\begin{array}{lccccc}
	\pi_i: & \vecx_i &\mapsto& w, \\
	\pi_j: & \vecx_j &\mapsto& w, \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	\piiw &\mapsto& \pi(\vecx_i), \\
	\piiw &\mapsto& \pi(\vecx_j), \\
	\end{array}
$$
respectively.
Since $\vecx_i \ne \vecx_j$, $\pi_i \ne \pi_j$.


Second, we show that the six permutations
$T(\pi,\vecy_1), \ldots, T(\pi,\vecy_6)$ are distinct.
Let $i \ne j$ for $i,j=1, \ldots, 6$; let $\pi_i
= T(\pi,\vecy_i)$ and $\pi_j = T(\pi,\vecy_j)$.  Then $\pi$ has
$$
	\begin{array}{lccccc}
	\pi:   & \vecy_i &\mapsto& \pi(\vecy_i), \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	\vecy_j &\mapsto& \pi(\vecy_j), \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	w &\mapsto& \piw; \\
	\end{array}
$$
$\pi_i$, and $\pi_j$ have
$$
	\begin{array}{lccccc}
	\pi_i: & \vecy_i &\mapsto& \piw, \\
	\pi_j: & \vecy_j &\mapsto& \piw, \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	w &\mapsto& \pi(\vecy_i), \\
	w &\mapsto& \pi(\vecy_j), \\
	\end{array}
$$
respectively.  Since $\vecy_i \ne \vecy_j$, by lemma \ref{lem:pineq}
$\pi(\vecy_i) \ne \pi(\vecy_j)$.  Since $\pi_i,\pi_j$ send $w$ to to different
sites, $\pi_i \ne \pi_j$.

Third, we show that the head-swaps of $\pi$ are distinct from the tail-swaps
of $\pi$.  Fix $\pi\in\SNpc$ and let $i,j \in \{1, \ldots, 6\}$.  Then $\pi$
has
$$
	\begin{array}{lccccc}
	\pi:   & \vecx_i &\mapsto& \pi(\vecx_i), \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	\vecy_j &\mapsto& \pi(\vecy_j), \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	\piiw &\mapsto& w &\mapsto& \piw; \\
	\end{array}
$$
$S(\pi,\vecx_i)$ and $T(\pi,\vecy_j)$ have
$$
	\begin{array}{rccccc}
	S(\pi,\vecx_i): & \piiw &\mapsto& \pi(\vecx_i), \\
	T(\pi,\vecy_j): & \vecy_j &\mapsto& \piw, \\
	\end{array}
\qquad
	\begin{array}{rccccc}
	\vecx_i &\mapsto& w &\mapsto& \piw; \\
	\piiw   &\mapsto& w &\mapsto& \vecy_j; \\
	\end{array}
$$
respectively.  Under these two permutations, $w$ has images $\piw$ and
$\vecy_j$, respectively, and preimages $\vecx_i$ and $\piiw$.  By definition
\ref{defn:worm_trivial}, the non-trivial head swap $S(\pi,\vecx_i)$ has
$\vecx_i \ne \piiw$ and the non-trivial tail swap $T(\pi,\vecy_j)$ has $\piw
\ne \vecy_j$.  Thus, $S(\pi,\vecx_i)$ and $T(\pi,\vecy_j)$ are distinct
permutations.
\end{proof}
\end{lem}

%\begin{lem}
%Head swaps and tail swaps are transitive on fibers.
%\begin{proof}
%This needs to be proved.  This key result will establish irreducibility of the
%chain.
%
%%% This was proved different, using SO irr along with SO irr => worm irr.
%\end{proof}
%\end{lem}

%% ----------------------------------------------------------------
\section{Explicit construction of the Markov matrix}

Transition probabilities were described as being proportional to $1 \wedge
\edh$.  We put the constants of proportionality to be:
\begin{itemize}
\item $a$ for head swaps and tail swaps;
\item $b$ for closer moves;
\item $c$ for opener moves.
\end{itemize}
For SO/SAR, we chose the normalizing factor $C$ solely by the row-normalization
condition.  Here, with a more complicated algorithm, we will choose the
normalizing factors to satisfy detailed balance.

The Markov matrix is now $(N+1)! \times (N+1)!$:
\begin{itemize}

\item A closed permutation transitions only to itself, or to any of the $N$
open permutations in the fiber above it.  Thus, there are $N+1$ non-zero
entries in $\pi$'s row of $\vecM$.

\item An open permutation transitions to any of the 12 open permutations
available by head-swapping or tail-swapping, or to itself, or to the closed
permutation at the base of its fiber.  Thus, there are 14 non-zero entries in
$\pi$'s row of $\vecM$.

\end{itemize}

\begin{defn}
For open $\pi$, let
$$
	\{\vecx_1, \ldots, \vecx_6\} = \{\vecx \in \Lambda: \|\vecx,\piiw\|=1\}
$$
and
$$
	\{\vecy_1, \ldots, \vecy_6\} = \{\vecy \in \Lambda: \|\piy,\piw\|=1\}.
$$
Then define
\begin{align*}
	R_S(\pi) &= \{S(\pi,\vecx_1), \ldots, S(\pi,\vecx_6)\}, \\
	R_T(\pi) &= \{T(\pi,\vecy_1), \ldots, T(\pi,\vecy_6)\}.
\end{align*}
These are the twelve open permutations reachable from $\pi$ via head swaps and
tail swaps, respectively (lemma \ref{lem:worm_distinct_dozen}).
For closed $\pi$, define
\begin{align*}
	R_O(\pi) &= \{O(\pi,\vecx_1), \ldots, O(\pi,\vecx_N)\}.
\end{align*}
These are the $N$ open permutations reachable from $\pi$ via opener moves.
\end{defn}

With these definitions, the entries of the transition matrix are as follows.

If $\pi$ is closed:
\begin{align*}
	M(\pi,\pi') &= \begin{cases}
		\displaystyle
		c\left(1 \wedge e^{-H(\pi')+H(\pi)}\right), &
			\pi' \in R_O(\pi);
			\\
		\displaystyle
		1 - \sum_{\pi' \in R_O(\pi)}
			c\left(1 \wedge e^{-H(\pi')+H(\pi)}\right),
			& \pi'=\pi; \\
		0, & \textrm{ otherwise}.
	\end{cases}
\end{align*}

If $\pi$ is open:
\begin{displaymath}
	M(\pi,\pi') = \begin{cases}
		\displaystyle
		a\left(1 \wedge e^{-H(\pi')+H(\pi)}\right),
			& \pi' \in R_S(\pi); \\
	%
		\displaystyle
		a\left(1 \wedge e^{-H(\pi')+H(\pi)}\right),
			& \pi' \in R_T(\pi); \\
	%
		\displaystyle
		b\left(1 \wedge e^{-H(\pi')+H(\pi)}\right),
			& \pi' = C(\pi); \\
	%
		\displaystyle
		1 - t(\pi), & \pi'=\pi; \\
	%
		0, & \textrm{ otherwise}
	\end{cases}
\end{displaymath}
where
\begin{align*}
	t(\pi) =
		\sum_{\pi''\in R_S(\pi) \cup R_T(\pi)}
			a\left(1 \wedge e^{-H(\pi'')+H(\pi)}\right)
		- b\left(1 \wedge e^{-H(C(\pi))+H(\pi)}\right)
\end{align*}

Row normalization for closed $\pi$: $c(1 \wedge \edh)$ is between 0 and $c$ so
$\sum\textrm{rest}$ is between 0 and $cN$.  Take
\begin{align}
	c &= 1/N.
\end{align}
Row normalization for open $\pi$ then gives
\begin{align}
	12a + b \le 1.
\end{align}

%% ----------------------------------------------------------------
\section{Correctness}
\label{sec:worm_correctness}

As discussed in sections \ref{sec:MCMC_recipe} and \ref{sec:SO_correctness}, we
need to prove irreducibility, aperiodicity, and detailed balance for the worm
Markov chain.

\begin{prop}[Irreducibility]
\label{prop:worm_irreducible}
The worm algorithm's Markov chain is irreducible.
\begin{proof}
This follows immediately from propositions \ref{prop:SO_irreducible}
and \ref{prop:worm_irr_if_SO_irr}:  namely, the worm's chain is irreducibility
if the SO's chain is, and moreover the SO's chain is irreducible.
\end{proof}
\end{prop}

\begin{prop}
\label{prop:worm_irr_if_SO_irr}
The worm algorithm's Markov chain is irreducible if the SO algorithm's Markov
chain is irreducible.
\begin{proof}
The key point is that the composition of an open, head swap, and close are
precisely an SO swap.  Let $\vecx$ and $\vecy$ be lattice points such that
$\pix$ and $\piy$ are nearest neighbors.  Starting with $\pi$, then applying an
open at $\vecx$, a head swap at $\vecy$, and a close, we have
$$
	\begin{array}{lccccc}
	\pi:                     & \vecx &\mapsto& \pix, \\
	\pi'   = O(\pi,  \vecx): & \vecx &\mapsto& w, \\
	\pi''  = S(\pi', \vecy): & \vecx &\mapsto& \piy, \\
	\pi''' = C(\pi''):       & \vecx &\mapsto& \piy, \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	\vecy &\mapsto& \piy, \\
	\vecy &\mapsto& \piy, \\
	\vecy &\mapsto& w, \\
	\vecy &\mapsto& \pix, \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	w &\mapsto& w;    \\
	w &\mapsto& \pix; \\
	w &\mapsto& \pix; \\
	w &\mapsto& w.
	\end{array}
$$

This shows that, if the SO algorithm is irreducible on $\SN$, the worm
algorithm is irreducible on $\SN$.  But then the worm algorithm is also
irreducible on $\SNp$:  fix an initial and final permutation; close the initial
permutation, if it is open, to obtain a closed permutation; use the preceding
argument to reach the closed permutation which lies under the fiber of the
desired final open permutation; do an open move (see lemma
\ref{lem:worm_opens_and_closes}) if the final permutation is open.
\end{proof}
\end{prop}

\begin{rem}
The worm algorithm has an additional degree of freedom.  If $\vecx$ and $\vecy$
are nearest-neighbor lattice sites, then the composition of an open at $\vecx$,
a tail swap at $\vecy$, and a close results in a similar swap of the jump
targets of $\vecx$ and $\vecy$:
$$
	\begin{array}{lccccc}
	\pi:                     & \vecx &\mapsto& \pix, \\
	\pi'   = O(\pi,  \vecx): & \vecx &\mapsto& w, \\
	\pi''  = T(\pi', \vecy): & \vecx &\mapsto& w, \\
	\pi''' = C(\pi''):       & \vecx &\mapsto& \piy, \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	\vecy &\mapsto& \piy, \\
	\vecy &\mapsto& \piy, \\
	\vecy &\mapsto& \pix \\
	\vecy &\mapsto& \pix, \\
	\end{array}
\qquad
	\begin{array}{lccccc}
	w &\mapsto& w;    \\
	w &\mapsto& \pix; \\
	w &\mapsto& \piy; \\
	w &\mapsto& w.
	\end{array}
$$
\end{rem}

\begin{prop}[Aperiodicity]
The worm algorithm's Markov chain 