%% ================================================================
\section{Introduction}
\label{sec:intro}

The random-spatial-permutation model arises in the study of the Bose gas.  Its
history includes Bose-Einstein, Feynman \cite{bib:feynman}, Penrose-Onsager
\cite{bib:penons}, S\"ut\H{o} \cite{bib:suto1,bib:suto2}, and
Betz-Ueltschi-Gandolfo-Ruiz \cite{bib:gru,bib:qmath,bib:bu07,bib:bu08}.  Such
random permutations arise physically when one symmetrizes the $N$-boson
Hamiltonian with pair interactions, then applies a multi-particle Feynman-Kac
formula and a cluster expansion \cite{bib:bu07,bib:bu08}.  Specifically, given
points $\vecx_1, \ldots, \vecx_N$ in the box $[0,L]^3$ and temperature $T$,
permutations $\pi$ are given weights proportional to the Gibbs factor
$e^{-H(\pi)}$ where
\begin{align}
\label{eqn:intro_H}
	H(\pi) &= \frac{T}{4} \sum_{i=1}^N \|\vecx_i-\vecx_{\pi(i)}\|^2
\end{align}
for the non-interacting case.  (The sum is scaled by temperature rather than
reciprocal temperature.  This surprising feature is opposite that of many
models in statistical mechanics.) These energy terms involve lengths of
permutation jumps; additional interaction terms take the form
\begin{align}
\label{eqn:intro_V}
	\sum_{i<j} V(\vecx_i, \vecx_{\pi(i)}, \vecx_j, \vecx_{\pi(j)}),
\end{align}
i.e. permutation jumps from sites $\vecx_i$ and $\vecx_j$ interact pairwise.
In the above-cited papers of Betz and Ueltschi, these may be approximated and
rearranged such that one obtains interaction terms of the form
\begin{align}
\label{eqn:intro_ell}
	\sum_{\ell=1}^N \alpha_\ell r_\ell(\pi)
\end{align}
where $r_\ell(\pi)$ counts the number of $\ell$-cycles of the permutation
$\pi$, and the coefficients $\alpha_\ell$ are \emph{cycle weights}.

When the Bose gas is cooled below a critical temperature $T_c$, there is a
phase transition: a macroscopic fraction of the bosons are found in the ground
state of the external potential, and such particles are said to participate in
a Bose-Einstein condensate.  In the permutation representation, this transition
manifests itself as the onset of long permutation cycles.  Bose and Einstein
computed the critical temperature for the non-interacting Bose gas.  The
critical temperature for liquid helium, where interparticle interactions are
strong, is lower than would be expected \cite{bib:baymetal} for non-interacting
atoms of the same density.  For weakly interacting systems, however, an
emerging consensus is that interactions \emph{increase} the critical
temperature.  See \cite{bib:baymetal,bib:su09} for surveys.  Concretely, for
interactions parameterized by some $\alpha$, one defines
\begin{align*}
	\Delta T_c(\alpha) &= \frac{T_c(\alpha)-T_c(0)}{T_c(0)}.
\end{align*}
It is well accepted that
\begin{align*}
	\lim_{\alpha\to 0} \Delta T_c(\alpha) = c \rho^{1/3}\alpha
\end{align*}
where $\rho$ is the particle density, i.e. that for small $\alpha$ the shift in
critical temperature is linear in $\alpha$.  What is more contentious, as
enumerated in the surveys cited above, is the value of the constant $c$.

The interaction terms (equation \ref{eqn:intro_V}) for the permutation
representation of the Bose gas are difficult to compute.  Moreover, it is
interesting to consider the model of random spatial permutations for its own
sake.  In \cite{bib:gru}, a simulational approach is taken for points held
fixed on the cubic unit lattice in the non-interacting case.  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.  That is, interactions are
of the form of \ref{eqn:intro_ell} with $\alpha_2 = \alpha$, where $\alpha$ is
related to a hard-core scattering length, and the remaining cycle weights are
zero.  In \cite{bib:bu08}, this approach is extended to a model in which all
the $\alpha_\ell$'s may vary, but with the hypothesis that $\alpha_\ell$ goes
to zero faster than $1/\log(\ell)$.  In this paper, we take a simulational
approach to points on the lattice, with cycle weights constant in $\ell$ ---
removing the small-cycle-weight hypothesis.  The shift in critical temperature
is nonetheless found to match that predicted by Betz and Ueltschi.

An outline of the paper is as follows.  Section \ref{sec:background} provides
background necessary to understand the results of the paper:  the probability
model is defined in subsection \ref{subsec:prb_model}; qualitative and
quantitative behavior of long cycles are discussed in subsections
\ref{subsec:qual_long_cycles} and \ref{subsec:ordps}, respectively.  Known
results and conjectures are listed in subsection \ref{subsec:conj}.  In section
\ref{sec:simulational_methods}, the simulational methods are presented.  The
\emph{swap-only} and \emph{swap-and-reverse} algorithms generate simulational
data; these algorithms are proved correct in subsections
\ref{subsec:SO_algorithm} through \ref{subsec:SO_correctness}.  The
finite-size-scaling method, which reduces the raw simulational data, is
summarized in subsection \ref{subsec:fss}.  Section \ref{sec:results} presents
the data and its analysis in full detail:  estimation of critical exponents and
critical temperature in subsections \ref{subsec:rho_hat_L} through
\ref{subsec:crossing}, verification of the finite-size-scaling hypothesis in
subsection \ref{subsec:collapse}, and final results in sections
\ref{subsec:Delta_Tc_and_const} through \ref{subsec:conclusions}.

%% ================================================================
\section{The model of random spatial permutations}
\label{sec:background}

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

%% ----------------------------------------------------------------
\subsection{The probability model}
\label{subsec:prb_model}

The state space is $\Omega_{\Lambda,N} = \Lambda^N \times \SN$, where
$\Lambda=[0,L]^3$ with periodic boundary conditions; point positions are $\vecX
= (\vecx_1, \ldots, \vecx_N)$ for $\vecx_1, \ldots, \vecx_N \in \Lambda$.  The
Hamiltonian takes one of two forms.  In the first, relevant to the Bose gas, we
have
\begin{align}
\label{eqn:results_Bose_H}
	H(\vecX,\pi) &= \frac{T}{4} \sum_{i=1}^N \|\vecx_i-\vecxpii\|^2
		+ \sum_{i<j} 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 temperature scale factor $T/4$, not $\beta/4$, is surprising but correct
for the Bose-gas derivation of the Hamiltonian.) 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:results_RCM_H}
	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 interest.

Different choices of $\alpha_\ell$ result in different models: The
\emph{non-interacting model} \cite{bib:gru} has $\alpha_\ell \equiv 0$.  The
\emph{two-cycle model} \cite{bib:bu07,bib:qmath}, has $\alpha_2=\alpha$ and
other cycle weights are zero.  The \emph{general-cycle model} has no
restrictions on $\alpha_\ell$.  In \cite{bib:bu08}, the \emph{small
cycle-weight case} of the general-cycle model is considered:  the only
restriction on $\alpha_\ell$ is that $\alpha_\ell$ goes to zero in $\ell$
faster than $1/\log\ell$.  The \emph{Ewens model}, treated in this paper (see
also \cite{bib:ewens}), is another special case of the general-cycle model: it
has $\alpha_\ell \equiv \alpha$ constant in $\ell$.

\begin{figure}
\begin{center}
\psfragscanon
\includegraphics[scale=1.0]{figures/pmt1.eps}
\caption{A spatial permutation.
	\label{fig:pmt1}}
\end{center}
\end{figure}

One may hold point positions fixed, e.g. on the unit lattice; this approach has
been taken for all simulations done up to the present by the author and by
Gandolfo \cite{bib:gru}, including specifically the work described in this
paper.  One obtains a Gibbs probability distribution on $\SN$:
\begin{align}
\label{eqn:results_RCM_P}
	Y(\Lambda,\vecX) &= {\sum_{\sigma\in\SN} e^{-H(\vecX,\sigma)}}, &
	P_{\Lambda,\vecX}(\pi) &= \frac{e^{-H(\vecX,\pi)}} {Y(\Lambda,\vecX)}.
\end{align}
(Alternatively, one may integrate over all positions in $\Lambda$, with a
resulting Gibbs distribution on $\SN$.  Here, several analytical results are
available \cite{bib:bu07,bib:bu08}.)
For a random variable $S(\pi)$, we have
\begin{align}
\label{eqn:E_pi}
	\bbE_{\Lambda,\vecX}[S] = \frac{1}{Y(\Lambda,\vecX)}
		\sum_{\pi\in\SN} S(\pi) e^{-H(\vecX,\pi)},
	\qquad\textrm{i.e.}\qquad
	\bbE_\pi[S(\pi)] = \sum_{\pi\in\SN} P(\pi) S(\pi).
\end{align}

%% ----------------------------------------------------------------
\subsection{Qualitative characterization of long cycles}
\label{subsec:qual_long_cycles}

One next inquires which permutations are typical in this temperature-dependent
probability distribution on $\SN$.  In this section we develop intuition; in
the next section, 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 high $T$, medium but subcritical $T$, and low $T$.
	\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 are large whenever any jump
has non-zero length.  For large but finite $T$, there are tiny islands of
2-cycles, 3-cycles, etc.  On the other hand, 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.

For $T$ above a critical temperature $T_c$, all cycles are short: two-cycles,
three-cycles, and so on.  We find $T_c \approx 6.86$ at $\alpha=0$, and
positive $\alpha$ terms increase $T_c$.  At $T_c$, though, there is a phase
transition: for $T<T_c$ jump lengths remain short but \emph{long cycles form}.
Quantitatively, let $\ellmax$ be the length of the longest cycle in $\pi$, with
$\bbE[\ellmax]$ its mean over all permutations.  We observe that for $T>T_c$,
$\bbE[\ellmax]$ is constant as $N\to\infty$: cycles remain finite.  That is to
say, $\bbE[\ellmax]/N$ goes to zero as $N\to\infty$.  For $T<T_c$, on the other
hand, $\bbE[\ellmax]$ scales with $N$, i.e.  $\bbE[\ellmax]/N$ approaches a
temperature-dependent constant as $N\to\infty$: there are arbitrarily long
cycles, or infinite cycles, in the infinite-volume limit.  See figure
\ref{fig:dot_plots} for depictions of typical permutations at high $T$,
subcritical $T$, and low $T$; 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$.
Note in particular that higher alpha shifts the order-parameter curve to the
right, with resulting upward shift in critical temperature $T_c$.

\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$ for finite systems, with
	$\alpha=0.000,0.001$.  Interactions increase the critical temperature.
	\label{fig:E_ellmax_N}}
\end{center}
\end{figure}

Feynman's claim for the Bose gas is that Bose-Einstein condensation occurs if
and only if there are infinite cycles in the infinite-volume limit.  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 the venerable question: how does the critical
temperature of Bose-Einstein condensation depend on inter-particle interaction
strength?

Obtaining a full answer to this notoriously difficult question is a long-term
project.  As an intermediate step, we here consider the Ewens cycle-weight
Hamiltonian with point positions on the unit lattice.  Through careful use of
MCMC algorithms, statistical analysis, and finite-size scaling, we are able to
quantify the dependence of critical temperature on interaction strength.

%% ----------------------------------------------------------------
\subsection{Quantitative characterization of long cycles}
\label{subsec:ordps}

Various order parameters may be defined; all of them may be used to locate the
critical temperature $T_c(\alpha)$.
%
The fraction $\bbE[\ellmax]/N$ discussed above will, for brevity, be hereafter
referred to as $f_M$.
%
The \emph{fraction of sites in long cycles}, $f_I$, is described in detail in
\cite{bib:gru}.
%
The \emph{correlation length} $\xi(T)$ is defined to be the spatial length of
the cycle containing a given point $\vecx$: for $T<T_c$, it blows up as $L$
increases.  Namely, we define
\begin{align*}
	s_\vecx(\pi) &= \|\pi(\vecx)-\vecx\| &
	&\tand&
	s(\pi) &= \frac{1}{N} \sum_{i=1}^N s_{\vecx_i}(\pi).
\end{align*}
The expectation over all $\pi$ of $s_\vecx$ is the same as $s$, of course;
in a Monte Carlo simulation, however, the latter yields a larger sample
size and thus a smaller error bar.  We use $\xi(T) = \bbE[s]$.
%
\emph{Winding numbers} count the number of $x,y,z$ wraps around the 3-torus
($\Lambda$ with periodic boundary conditions).  Specifically,
the winding number of a permutation $\pi$ is the triple
\begin{align}
\label{eqn:def_wno}
	\vecW &= (W_x, W_y, W_z) = \sum_{i=1}^N \tilde{d}(\vecx_{\pi(i)} - \vecx_i).
\end{align}
We also write
\begin{align*}
	\vecW^2 &= \vecW \cdot \vecW = W_x^2 + W_y^2 + W_z^2.
\end{align*}
Equation \ref{eqn:def_wno} 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 distances
$\tilde{d}$ in equation \ref{eqn:def_wno} are to be interpreted in terms of
lifts to a separate Euclidean chart for each site $\vecx_i$.  The \emph{scaled
winding number} \cite{bib:pc87} is
$$
	f_S=\frac{\pig{\vecW^2} L^2}{3\beta N}.
$$
%
Lastly, the order parameter $f_W$ is the fraction of sites which participate in
winding cycles.

The order parameters $f_I(T)$, $f_S(T)$, and $f_W(T)$ show behavior similar to
$f_M := \bbE[\ellmax]/N$ (figure \ref{fig:E_ellmax_N}):  asymptotically as
$N\to\infty$, they are zero to the left of $T_c$ and non-zero to the right.
For finite $N$, the curves remain analytic:  finite-size effects persist.  The
inverse correlation length $1/\xi(T)$, on the other hand, is zero to the left
of $T_c$ and non-zero to the right.

Our goal is to quantify the dependence of $T_c$ on $\alpha$, where
\begin{align}
\label{eqn:Delta_Tc}
	\Delta T_c(\alpha) = \frac{T_c(\alpha)-T_c(0)}{T_c(0)}.
\end{align}
Known results and conjectures are formulated quantitatively in terms of
$
	\lim_{\alpha\to 0} \Delta T_c(\alpha).
$

%% ----------------------------------------------------------------
\subsection{Known results and conjectures}
\label{subsec:conj}

Known results for point locations averaged over the continuum 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}.  (This taps into a long and controversial history in the
physics literature: see \cite{bib:baymetal} or \cite{bib:su09} for surveys.)
The critical $(\rho,T,\alpha)$ manifold relates $\rho_c$ to $T_c$.
Specifically,
\begin{align}
\label{eqn:rho_c}
	\rho_c(\alpha)
	&\approx \sum_{\ell \ge 1} e^{-\alpha_\ell}
		\int_{\R^3} e^{-\ell \, 4\pi^2 \beta \|\veck\|^2}\,d\veck %\\
	= \frac{1}{(4\pi\beta)^{3/2}}
		\sum_{\ell \ge 1} e^{-\alpha_\ell} \ell^{-3/2} \\
	%
	%r_2: \rho_c(\alpha)
	%	&= \rho_c(0) - \frac{(1-e^{-\alpha})}{(8\pi\beta)^{3/2}} \\
	%
	\Delta T_c(\alpha) &\approx c \rho^{1/3} \alpha, \quad
	\textrm{for}\; \alpha \approx 0.
\end{align}
Using this formula for constant cycle weights $\alpha_\ell \equiv \alpha$ and
for lattice density $\rho=1$, we have
\begin{align}
\label{eqn:theo_Tc}
	\rho_c \approx \frac{\zeta(3/2) e^{-\alpha}}{(4\pi\beta)^{3/2}},
	\qquad
	T_c \approx \frac{4 \pi e^{2\alpha/3}}{\zeta(3/2)^{2/3}}
		\approx 6.626 e^{2\alpha/3}, \qquad
	c \approx 0.667.
\end{align}
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.

For $\alpha_\ell \equiv 0$ (the non-interacting model), $\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 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$ is that
same number for the non-interacting case $\alpha_\ell \equiv 0$.  The intuition
is that long cycles are ``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.

We suspect that the fine details of point positions are unimportant for the
shift in critical temperature.  Thus, $\Delta T_c(\alpha)$ on the lattice
should be similar to that on the continuum, if small cycle weights are used.
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.

%% ================================================================
\section{Simulational methods}
\label{sec:simulational_methods}

We run Markov chain Monte Carlo experiments for various values of $L$, $T$, and
interaction strength $\alpha$.  For each parameter combination, we generate $M$
typical permutations $\pi_1, \ldots, \pi_M$ from the stationary distribution,
using MCMC algorithms described below, and we compute random variables $X_i =
X(\pi_i)$.  (The values of $M$ used are $10^5$ away from $T_c$, and $10^6$ near
$T_c$ where sample variance is higher.) We find the sample mean and estimate
the variance of the sample mean.  The correlation of the $X_i$'s complicates
the latter.  Finite-size scaling compensates for finite-size effects:
mathematically, we are interested in estimating infinite-volume quantities
based on finite-volume numerical experiments.

%% ----------------------------------------------------------------
\subsection{The swap-only algorithm}
\label{subsec:SO_algorithm}

Recall from equation \ref{eqn:E_pi} that the expectation of a random variable
$S$ (such as $\xi$, $f_M$, $f_W$, $f_I$, $f_S$) is
$$
	\bbE_\pi[S(\pi)] = \sum_{\pi\in\SN} P(\pi) S(\pi).
$$
The number of permutations, $N!$, grows intractably in $N$.  As is typical in
Markov chain Monte Carlo methods \cite{bib:berg,bib:landaubinder}, one contents
oneself with a smaller number of samples: the expectation is instead estimated
by summing over some number $M$ ($10^5$ or $10^6$) of typical permutations.

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

The \emph{swap-only algorithm} for sampling from the Gibbs distribution
(equation \ref{eqn:results_RCM_P}) is as follows:
\begin{itemize}
\item Start with the identity or uniform-random permutation.
\item Sweep through sites $\vecx$ of the lattice in either lexical or
uniform-random order.
\item For each site $\vecx$, do a \emph{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:gkswap}.)
	\item With probability proportional to $\min\{1, \edh\}$
	where $\Delta H = H(\pi') - H(\pi)$, accept the change.
	(If the change is rejected, $\pi'=\pi$.)
	\end{itemize}

\item After each sweep, obtain a value of each random variable for inclusion in
computation of its sample mean.

\end{itemize}

One starts accumulating data only after a suitable number of
\emph{thermalization sweeps}.  The idea is that the initial, identity
permutation is not typical, nor are the first few afterward.  The integrated
autocorrelation time \cite{bib:berg} of system energy $H$ gives an idea of how
many Metropolis sweeps should be discarded before the permutations become
typical.  Also, one may examine $H$ to ensure that it has reached its long-term
average value.  This is explained in detail in \cite{bib:kerl}.  We next prove
correctness of this algorithm.

%% ----------------------------------------------------------------
\subsection{Explicit construction of the Markov matrix}
\label{subsec:SAR_Markov_matrix}

For section \ref{subsec:SO_correctness} we will need an explicit construction of
the Markov matrix corresponding to the swap-only algorithm as described in
section \ref{subsec:SO_algorithm}.  The Markov perspective on the algorithm is
that the distribution $P^{(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*}
	P^{(k+1)}(\pi') &= \sum_{\pi \in \SN} P^{(k)}(\pi) M(\pi, \pi')
\end{align*}
or, in matrix/vector notation,
\begin{align*}
	\vecP^{(k+1)} &= \vecP^{(k)} \vecM.
\end{align*}
In this section we precisely describe the matrix $\vecM$; in section
\ref{subsec:SO_correctness} we show that $\vecP^{(k)}$ approaches the
Gibbs distribution (equation \ref{eqn:results_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}
\label{defn:SO_trivial}
A swap on sites $\vecx$ and $\vecy$ is \emph{trivial} if $\vecx = \vecy$.
\end{defn}

\begin{defn}
Lattice sites $\vecx,\vecy$ are \emph{nearest-neighbor} if
$\|\vecx-\vecy\| = 1$.
\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 the normalization factor $C$ is, we observe the following:
\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}

We require rows sums to be 1 for stochastic matrices so we 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 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}

%% ----------------------------------------------------------------
\subsection{Correctness of the swap-only algorithm}
\label{subsec:SO_correctness}

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

We note the following terminology: \emph{detailed balance} is the same as
\emph{reversibility}\index{reversible}.  Also, an irreducible, aperiodic chain
on a finite state space is called \emph{ergodic}.  Also note from Markov-chain
theory 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{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 numbers, and the empirical fact that the
swap-only algorithm changes them only rarely.  The chain is irreducible but
various non-zero transition probabilities can still be very small.
\end{rem*}

\begin{defn}
\label{defn:period}
A permutation $\pi$ has \plainidx{period} $p$ if any return to $\pi$ must occur
at multiples of $p$ steps.  Precisely, let $\Pi_n$ be the random variable which
is the $n$th value of the Markov chain.  Then
\begin{align*}
	p = \gcd\{n : P(\Pi_n = \pi \mid \Pi_0 = \pi) > 0\}.
\end{align*}
Furthermore, the chain is said to be aperiodic if all states are aperiodic.
\end{defn}

\begin{prop}[Aperiodicity]
\label{prop:SO_aperiodic}
The swap-only algorithm's Markov chain is aperiodic.
\begin{proof}
From Markov-chain theory, we know that all states in a recurrence class have
the same 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 is in fact
the case:  The identity has zero energy (equations \ref{eqn:results_RCM_H}).
Any state proposed as a transition from the identity is necessarily a
two-cycle, for which the distance terms in equation \ref{eqn:results_RCM_H} are
positive.  The cycle-weight terms or interaction terms are non-negative.  Thus
$\Delta H > 0$, and the 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:results_RCM_P}) and the 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 swap-only moves.

\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:gkcycles}}
\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 \bfidx{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{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:gkcycles}).
\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
\qand
	\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
\qand
	\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 &= \ellypi = \ellxpi + \ellypi.
\end{align*}
\end{rem}

%% ----------------------------------------------------------------
\subsection{Winding cycles and the swap-and-reverse algorithm}
\label{subsec:wno_SAR}

The propositions of section \ref{subsec:SO_correctness} showed that the
swap-only algorithm is correct --- 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 that
the swap-only algorithm almost always generates permutations with zero winding
number.

This problem, and a partial solution, is explained intuitively by figure
\ref{fig:gkwinding} and rigorously in \cite{bib:kerl}.  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 (with full details in
\cite{bib:kerl}), winding numbers of even parity can be generated.

\begin{figure}
\begin{center}
\psfragscanon
\includegraphics[scale=0.9]{figures/gkwinding.eps}
\caption{Conservation of winding number in the swap-only algorithm.
	\label{fig:gkwinding}}
\end{center}
\end{figure}

Our current best algorithm (swap-and-reverse) has two types of sweeps:  (1) For
each lattice site, do a Metropolis step as above.  (2) For each cycle in the
permutation, reverse the direction of the cycle with probability $1/2$.  This
permits winding numbers of even parity in each of the three axes.

We have experimented with various methods to obtain winding numbers of all
parities.  The creation or destruction of a winding cycle is a non-local
update; one is reminded of the Swendsen-Wang algorithm for the Ising model.
However, 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.

We have also created a \emph{worm algorithm}, inspired by approaches to this
same winding-number problem in path-integral Monte Carlo methods
\cite{bib:bps06,bib:pst98}.  That is, a permutation loop is selected at random
and then cut open at a randomly selected point.  The resulting worm is allowed
to move around $\Lambda$ via Metropolis moves; eventually, it closes again.
This worm algorithm has an elegant theory and correctness proof
\cite{bib:kerl}; yet, it has an unacceptably long stopping time for loop
closure, and none of our attempts to remedy the stopping-time problem have
satisfied detailed balance.

% At present, we content ourselves with the swap-and-reverse algorithm.  We
% next make an entropy argument suggesting that sampling permutations with even
% winding number does not harm our final result.
%
% \begin{prop}
% [statement goes here.]
% \begin{proof}
% [proof goes here.  Type up the handwritten notes.]
%
% I do not remember exactly what we discussed. But it seems to me that the
% critical temperature is not affected by restricting to even winding numbers.
%
% The log of the partition function is roughly the same, up to an additive
% constant. Dividing by the volume, we get the same in the infinite volume
% limit. This means that the thermodynamic potentials are exactly the same. One
% expects that the critical temperature corresponds to non-analyticity of the
% thermodynamic potential, so it does not change.  Nor does the difference
% between critical temperatures, obviously.
%
% Let $W = \#\Omega$.  $S = \log (\Omega/8) = \log\Omega - \log 8$.
% Also $N \sim \log\Omega$.
%
% $W \sim \exp(a N + b N^{2/3} + c N^{1/3})$.  ``bulk'', vol, surface:
% finite-size effects.
%
% $S = \log W \sim a N + c N^{2/3}$.  Shift by $\log 8$ is less than
% finite-size effects.
%
% \end{proof}
% \end{prop}
%
% We also note that 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.

At present, we content ourselves with the swap-and-reverse algorithm.  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; results
obtained in section \ref{sec:results} using each of the five order parameters
are compatible.

%% ----------------------------------------------------------------
\subsection{Finite-size scaling}
\label{subsec:fss}

Finite-size scaling takes the form of a hypothesis, or rather a set of
hypotheses, which is tested against the data.  See also \cite{bib:cggp} for a
nice survey.

We have an infinite-volume random variable $S(T)$, e.g. any of the order
parameters defined in section \ref{subsec:ordps}.  The finite-volume quantity
is $S_L(T)$.  Define $t=(T-T_c)/T_c$.  Examine, say, $0.99 < t < 1.01$.
The first hypothesis is that the correlation length $\xi(T)$ follows a power
law
$$
	\xi(T) \sim |t|^{-\nu}, \quad T \to T_c
$$
For the infinite-volume quantity, we expect a power-law behavior
$$
	S(T) \sim t^{\rho}, (-t)^{\rho}, \qor |t|^\rho.
$$
(The domain of validity is $t<0$ or $t>0$ depending on whether the order
parameter $S$ is left-sided or right-sided, respectively.) One moreover
hypothesizes that for $T$ near $T_c$, $S_L(T)$ and $S(T)$ are related by a
\emph{universal function} $Q_S$ which depends on $T$ only through the ratio
$L/\xi$:
\begin{align}
\label{eqn:fss}
	S_L(T) &= L^{-\rho/\nu} Q_S(L^{1/\nu} t)
	\sim L^{-\rho/\nu} Q_S((L/\xi)^{1/\nu}).
\end{align}

%% ================================================================
\section{Results}
\label{sec:results}

Here we complete the steps sketched in section \ref{subsec:fss}.  The flow of
data and uncertainties are as follows:

\begin{itemize}

\item Markov chain Monte Carlo simulations, with error bars determined using
the method of integrated autocorrelation time \cite{bib:berg}, yield $S_L(T,
\alpha)$ data points.  There are five order parameters $S$, six values of $L$
(30, 40, 50, 60, 70, 80), nine values of $\alpha$, and a few dozen
values of $T$ for each $\alpha$.

\item For each $S$, $L$, and $\alpha$, we use $S_L(T,\alpha)$ values for all
available values of $T$ and $\alpha$ to estimate $\rhohat_S(L)$.  (Critical
exponents are assumed to be independent of $\alpha$ for small $\alpha$, or with
weak enough dependence on $\alpha$ that that dependence is lost in the noise.)
Error bars may be propagated from the MCMC simulations, or computed from
regression uncertainties.

\item Extrapolating $\rhohat_S(L)$ in $L \to \infty$ results in the five
estimated critical exponents $\rhohat_S$.  Uncertainties are computed from the
regression analysis.

% Uncertainties may be propagated from the previous step using the method of
% synthetic-data \cite{bib:nr}.

\item Once the critical exponents are estimated, we obtain
$\hat{T}_{c,S}(\alpha)$ for each of the five order parameters $S$ and for each
$\alpha$.  Uncertainties are computed by visual inspection of the crossing
plots discussed in section \ref{subsec:crossing}.

\item Once the critical exponents and $T_c$ are known, one should be able to
obtain plots of the universal function $Q_S$ which is, up to sampling
variability, independent of $L$, $T$, and $\alpha$.  This verifies the
finite-size-scaling hypothesis.

\item The shift in reduced critical temperature is as in equation
\ref{eqn:Delta_Tc}.  Error bars are computed from regression uncertainties.

% Uncertainties are computed by the method of synthetic data \cite{bib:nr}.

\end{itemize}

%% ----------------------------------------------------------------
\subsection{Determination of $L$-dependent critical exponents}
\label{subsec:rho_hat_L}

For each of order parameter $S$, interaction parameter $\alpha$, and box length
$L$, we examine all $S(L,T,\alpha)$ data for which $S > \veps$, with $\veps$
taken from the plots to ensure that we examine the portions of the curves
corresponding to non-zero order parameter in the infinite limit (see figure
\ref{fig:f_M_and_recip_xi}).  For $1/\xi$, this means $T>T_c$; for the other
four order parameters, this means $T<T_c$.  From plots such as those in figure
\ref{fig:f_M_and_recip_xi}, we choose $\veps$ to be
$0.1$  for $1/\xi$,
$0.01$ for $f_M$,
$0.01$ for $f_I$,
$0.05$ for $f_S$, and
$0.01$ for $f_W$.
For each $S$, $\alpha$, and $L$, we then apply linear regression to
$S(L,T)^{1/\rho_S}$ for varying $\rho_S$.  We find $\rhohat_S(L)$ which
optimizes the correlation coefficient \cite{bib:young} of the linear
regression.  Results are shown in figure \ref{fig:smiley_and_rhoestalpha}.
Given $\rhohat_S(L)$ along with its corresponding linear-regression parameters
$m$ and $b$, we may plot a power-law fit to the simulational data.  One such
comparison plot is shown in figure \ref{fig:ordp_fit_comparison}.

% Results are shown in table \ref{table:rhohat_S_L} and figure
% \ref{fig:smiley_and_rhoestalpha}.

\begin{figure}
\begin{center}
\psfragscanon
\includegraphics[scale=0.50]{../tests/fss3a/resultsfigs/graph_rell_alpha0_0.000_E_ellmax_N_small.eps}
\includegraphics[scale=0.50]{../tests/fss3a/resultsfigs/graph_rell_alpha0_0.001_E_ellmax_N_small.eps}
%
\includegraphics[scale=0.50]{../tests/fss3a/resultsfigs/grapheb_rell_alpha0_0.000_recip_xi_small.eps}
\includegraphics[scale=0.50]{../tests/fss3a/resultsfigs/grapheb_rell_alpha0_0.001_recip_xi_small.eps}
\caption{Order parameters $f_M$ and $1/\xi$ for $L=40,60,80$ and $\alpha=0.000$
and $0.001$.  The remaining order parameters $f_S$, $f_W$, and $f_I$ behave
similarly to $f_M$ but with not all with the same critical exponents.
	\label{fig:f_M_and_recip_xi}}
\end{center}
\end{figure}

% \begin{figure}
% \begin{center}
% \psfragscanon
% \includegraphics[scale=0.50]{../tests/fss3a/data/smiley_fS_L_80_alpha_0.000.eps}
% \caption[Determination of critical exponent $\hat{\rho}_S(L,\alpha)$ for $f_S$.]
% 	{Determination of critical exponent $\hat{\rho}_S(L,\alpha)$ for order
% 	parameter $f_S$, as the value which minimizes linear-regression error for
% 	$S_L(T,\alpha)^{1/\rho}$.  Visually, one sees $\hat{\rho}_S(L=80,
% 	\alpha=0.0) \approx 0.59$.
% 	\label{fig:smiley}}
% \end{center}
% \end{figure}
%
% \begin{table}
% \input{../tests/fss3a/rhoestalpha.tex}
% \caption{Estimated critical exponents for varying $L$.
% 	\label{table:rhohat_S_L}}
% \end{table}
% 
% \begin{figure}
% \begin{center}
% \psfragscanon
% %\includegraphics[scale=0.7]{figures/rhoestalpha.eps}
% \includegraphics[scale=0.5]{../tests/fss3a/rhoestalpha.eps}
% \caption{Estimated critical exponents for $L=30,40,50,60,70,80$.
% 	\label{fig:rhoestalpha}}
% \end{center}
% \end{figure}

\begin{figure}
\begin{center}
\psfragscanon
\includegraphics[scale=0.50]{../tests/fss3a/data/smiley_fS_L_80_alpha_0.000.eps}
\includegraphics[scale=0.50]{../tests/fss3a/rhoestalpha.eps}
\caption[Determination of critical exponent.]
	{On the left:  determination of critical exponent
	$\hat{\rho}_S(L,\alpha)$ for order parameter $f_S$, as the value which
	minimizes linear-regression error for $S_L(T,\alpha)^{1/\rho}$.  Visually,
	one sees $\hat{\rho}_S(L=80, \alpha=0.0) \approx 0.59$.
	%
	On the right:  estimated critical exponents for $L=30,40,50,60,70,80$.
	%
	\label{fig:smiley_and_rhoestalpha}}
\end{center}
\end{figure}

\begin{table}
\input{../tests/fss3a/fMfIvalues.tex}
\caption{$f_M/f_I$ as a function of $\alpha$.  An upward trend is
	visible, but it is not pronounced.
	\label{table:fMfI_alpha}}
\end{table}

\begin{figure}
\begin{center}
\psfragscanon
\includegraphics[scale=0.5]{../tests/fss3a/resultsfigs/fS_alpha_0.000_comp.eps}
\caption[Power-law fit vs. simulational order-parameter data.]
	{Power-law fit vs. simulational data for order parameter $f_S$,
	$\alpha=0.000$
	\label{fig:ordp_fit_comparison}}
\end{center}
\end{figure}

%% ----------------------------------------------------------------
\subsection{Extrapolation of critical exponents for the infinite-volume limit}
\label{subsec:rhohat_extrap}

Next, for each $S$, given estimates $\rhohat_S(L)$ for increasing values of
$L$, we plot $\rhohat_S(L)$ versus $1/L$.  The vertical intercept of this plot
estimates the infinite-volume exponent $\rhohat_S(\alpha)$.  (See figure
\ref{fig:smiley_and_rhoestalpha}.) Results are shown in table \ref{table:rhohat_S}.

\begin{table}
\input{../tests/fss3a/rhoest.tex}
\caption{Extrapolated estimates of the infinite-volume critical exponents,
	found from the vertical intercept of figure \ref{fig:smiley_and_rhoestalpha}.
	\label{table:rhohat_S}}
\end{table}

%% ----------------------------------------------------------------
\subsection{Determination of critical temperature}
\label{subsec:crossing}

% $fss3a/crossing.py

Given the above estimators of the critical exponents, the \emphidx{crossing
method} \cite{bib:cggp} estimates $T_c(\alpha)$.  Namely, we plot
$L^{\rhohat/\nuhat} S_L(T)$ as a function of $T$.  At $T=T_c$ we have $t=0$ and
$L^{\rho/\nu} S_L(T)  = Q_S(0)$, regardless of $L$ (equation \ref{eqn:fss}).
Thus, these curves will cross (approximately, due to sampling variability) at
$T=T_c$.  If they do not, the finite-size-scaling hypothesis is not verified.
(Note in particular that for order parameter $1/\xi$ whose critical exponent is
$\nu$, we apply the crossing method to $L S_L(T)$ as a function of $T$:  thus,
the $T_c(\alpha)$ estimate using $1/\xi$ is independent of $\nuhat$.) See for
example figure \ref{fig:fS_crossing}.
Results are in figure \ref{fig:Tc_alpha}.

%Results are in table \ref{table:Tc_alpha}.

\begin{figure}
\begin{center}
\psfragscanon
\includegraphics[scale=0.53]{../tests/fss3a/data/crossing_raw_fI_alpha_0.000.eps}
\includegraphics[scale=0.53]{../tests/fss3a/data/crossing_raw_zoom_fI_alpha_0.000.eps}

\includegraphics[scale=0.53]{../tests/fss3a/data/crossing_fit_fI_alpha_0.000.eps}
\includegraphics[scale=0.53]{../tests/fss3a/data/crossing_fit_zoom_fI_alpha_0.000.eps}
\caption
	[Estimation of $T_c(\alpha)$ for order parameter $f_I$.]
	{The crossing method to estimate $T_c(\alpha)$ for order parameter $f_I$,
	with $\rhohat$ and $\nuhat$ as above: $T_c(\alpha)$ corresponds to the
	horizontal coordinate of the intersection point of the plots.
	The upper-right-hand plot is a close-up of the upper-left-hand plot.
	The lower pair have raw MCMC data replaced with power-law-fit data.
	\label{fig:fS_crossing}}
\end{center}
\end{figure}

% \begin{table}
% \input{../tests/fss3a/Tcest.tex}
% \caption{Estimated $T_c(\alpha)$'s.
% 	\label{table:Tc_alpha}}
% \end{table}

%% ----------------------------------------------------------------
\subsection{Verification of finite-size-scaling hypothesis}
\label{subsec:collapse}

Now that we have estimated $\rho_S$, $\nu$, and $T_c(\alpha)$ for each of the
five order parameters $S$, we  may plot $L^{\rho_S/\nu} S_L(T,\alpha)$ as a
function of $L^{1/\nu}t$.  This is a plot of the scaling function $Q_S$.  If
the hypothesis is correct, the curves for all $L$ should coincide, or collapse,
to within sampling error --- which they do (e.g. figure \ref{fig:fS_collapse}).

\begin{figure}
\begin{center}
\psfragscanon
\includegraphics[scale=0.6]{../tests/fss3a/data/collapse_raw_fS_alpha_0.000.eps}
\caption
	{Collapse plot for order parameter $f_S$.
	\label{fig:fS_collapse}}
\end{center}
\end{figure}

%% ----------------------------------------------------------------
\subsection{Determination of the shift in critical temperature}
\label{subsec:Delta_Tc_and_const}

% @@@ update manually

As discussed in section \ref{subsec:conj}, we are seeking a linear relationship
between $\Delta T_c(\alpha)$ and $\alpha$, with constant $c$.  This can be
visualized in figure \ref{fig:Delta_Tc_alpha}, which is obtained from the
$T_{c,S}(\alpha)$ data of figure \ref{fig:Tc_alpha} using equation
\ref{eqn:Delta_Tc}.  Performing a linear regression with error bars
\cite{bib:young} on all the $(\alpha,\Delta T_c(\alpha))$ data points, we find
\begin{align*}
	c &= 0.665 \pm 0.067.
\end{align*}
This is the slope of the linear fit; the vertical intercept is $0.00026 \pm
0.00012$ which is consistent with equation \ref{eqn:Delta_Tc}.  Within
experimental uncertainty, this result, for points on the lattice with Ewens
cycle-weights, matches the $c$ value of equation \ref{eqn:theo_Tc} for point
positions varying on the continuum with small-cycle-weight interactions.

\begin{figure}
\begin{center}
\psfragscanon
\includegraphics[scale=0.6]{../tests/fss3a/Tcest.eps}
\caption
	{Critical temperature as function of $\alpha$.
	\label{fig:Tc_alpha}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\psfragscanon
\includegraphics[scale=0.6]{../tests/fss3a/Delta_Tcest.eps}
\caption
	[Shift in critical temperature as function of $\alpha$.]
	{Shift in critical temperature as function of $\alpha$.
	Recall from equation \ref{eqn:Delta_Tc} that
	$\Delta T_c(\alpha) = \frac{T_c(\alpha)-T_c(0)}{T_c(0)}$.
	\label{fig:Delta_Tc_alpha}}
\end{center}
\end{figure}

%% ----------------------------------------------------------------
\subsection{Constancy of the macroscopic-cycle quotient}
\label{subsec:mcq_const}

As discussed in section \ref{subsec:conj}, we hypothesize that the
macroscopic-cycle quotient $f_M / f_I$ in the infinite-volume limit is
dependent on $\alpha$ but is constant in $T$ where it is defined, i.e.  for
$T<T_c$ since $f_I=0$ for $T>T_c$.  This may be visualized by comparing figures
such as \ref{fig:E_ellmax_N}:  one sees that $f_M$ and $f_I$ appear to have the
same critical exponent.  Alternatively, one may plot the ratio $f_M/f_I$
(figure \ref{fig:fMfI_curves}).  In the infinite-volume limit, $f_I$ is zero
for $T>T_c$ and so we are interested only in the values of the quotient for
$T<T_c$.  In that region, the quotient does indeed appear to be constant in
$T$.

We test this constancy hypothesis as follows.  The respective critical
exponents are $\rho_M$ and $\rho_I$.  The estimators are $\rhohat_M$ and
$\rhohat_I$, computed by averaging over several different values of $L$ and
$\alpha$ as described in section \ref{subsec:rhohat_extrap}.  Treating these
estimators as normally distributed (as justified by the raw data), we obtain
the standard deviations of the $\rhohat_{M,I}(L,\alpha)$ samples, along with the
standard deviations of the means $\rhohat_{M,I}$:
\input{../tests/fss3a/fMfIseparate.tex}
The difference $\rhohat_M - \rhohat_I$ is also normally distributed about the
true mean $\rho_M - \rho_I$, but $\rhohat_M$ and $\rhohat_I$ are not
independent since they are sample means of random variables computed from the
same Markov chain Monte Carlo sequence of permutations.  Thus
\begin{align*}
	\Var(\rhohat_M - \rhohat_I) &= \Var(\rhohat_M) + \Var(\rhohat_I)
		- 2 \Cov(\rhohat_M, \rhohat_I).
\end{align*}
Computing the sample covariance of the $\rhohat_M(L,\alpha)$ and
$\rhohat_I(L,\alpha)$ data series, we obtain the covariance and resulting
standard error $s_d$ of the difference
\input{../tests/fss3a/fMfIcov.tex}
Normalizing, we find
\input{../tests/fss3a/fMfInormalize.tex}
We hypothesize $\rho_M -  \rho_I = 0$; the estimated value $\rhohat_M -
\rhohat_I$ lies comfortably within a standard deviation of this.  We note,
moreover, that the value of $f_M/f_I$, while constant in $T$, trends upward
with $\alpha$ (see table \ref{table:fMfI_alpha} and figure
\ref{fig:fMfI_alpha}).  This merits further investigation.

\begin{figure}
\begin{center}
\psfragscanon
%\includegraphics[scale=0.50]{../tests/fss3a/data/graph_rell_alpha0_0.000_fMfI.eps}
%\includegraphics[scale=0.50]{../tests/fss3a/data/graph_rell_alpha0_0.002_fMfI.eps}
%\includegraphics[scale=0.50]{../tests/fss3a/data/graph_rell_alpha0_0.004_fMfI.eps}
\includegraphics[scale=0.55]{../tests/fss3a/resultsfigs/graph_rell_alpha0_0.000_fMfI.eps}
\includegraphics[scale=0.55]{../tests/fss3a/resultsfigs/graph_rell_alpha0_0.002_fMfI.eps}
%\includegraphics[scale=0.50]{../tests/fss3a/resultsfigs/graph_rell_alpha0_0.004_fMfI.eps}
\caption
	{Macroscopic-cycle quotient $f_M/f_I$ for $\alpha=0.000, 0.002$.
	\label{fig:fMfI_curves}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\psfragscanon
\includegraphics[scale=0.50]{../tests/fss3a/fMfIvalues.eps}
\caption{$f_M/f_I$ as a function of $\alpha$.
	\label{fig:fMfI_alpha}}
\end{center}
\end{figure}

%% ----------------------------------------------------------------
\subsection{Conclusions}
\label{subsec:conclusions}

% @@@ update manually

(1) For annealed point positions, equation \ref{eqn:theo_Tc} gives $T_c(0)
\approx 6.625$.  Our result $T_c(0) = 6.863 \pm 0.003$ unambiguously shows that
the lattice structure modifies the critical temperature, even in the
non-interacting ($\alpha=0$) case.

% @@@ update manually

(2) As detailed in section \ref{subsec:Delta_Tc_and_const}, we find that the
reduced shift in critical temperature as a function of interaction parameter
$\alpha$ is
\begin{align*}
	\Delta T_c(\alpha) &\approx \frac{T_c(\alpha) - T_c(0)}{T_c(0)}
		= c \alpha
\end{align*}
with
\begin{align*}
	c &= 0.665 \pm 0.067.
\end{align*}
This is compatible (section \ref{subsec:conj}) with the related result of
\cite{bib:bu08}.  Even though the lattice structure changes the critical
temperature (conclusion 1), the \emph{shift} in critical temperature is
unaffected.

% @@@ update manually

(3) As described in section \ref{subsec:conj}, Shepp and Lloyd
\cite{bib:shepplloyd} find that $\bbE[\ellmax]/N \approx 0.6243$ for
uniform-random (non-spatial) permutations.  For spatial permutations, we define
a macroscopic-cycle quotient $\bbE[\ellmax]/N f_I$ which is the ratio of mean
maximum cycle length as a fraction of the number of sites in long cycles.  Our
result (table \ref{table:fMfI_alpha}) is compatible with that of Shepp and
Lloyd for the non-interacting case, with an increase which appears to be linear
as a function of interaction parameter $\alpha$.

%% ================================================================
\section{Future work}
\label{sec:future_work}

\textbf{Macroscopic-cycle quotient:} Now that the $\alpha$-dependence of the
macroscopic-cycle quotient's constant upon $\alpha$ has been found empirically,
one would next like to explain that dependence analytically.

\textbf{Non-asymptotic algorithm correctness:} The detailed-balance correctness
proof of the swap-and-reverse algorithm shows that there is a non-zero
transition probability between all pairs of permutations.  However, those
non-zero transition probablities can be quite small.  As the number of
Metropolis steps goes to infinity, asymptotically all permutations can be
reached; more interesting is the question of which permutations are actually
reachable in reasonable simulation time.  One answers this question empirically
simply by running simulations; perhaps this suffices.  For the system discussed
in this paper, as well as for other systems studied using MCMC methods, it
would be useful to have a non-asymptotic correctness theory.

\textbf{Winding numbers of all parities:} Ideally, one would have an algorithm
to permit odd winding numbers, as discussed in section \ref{subsec:wno_SAR}.

\textbf{Bose-gas Hamiltonian:} Sampling from the true Bose-gas distribution
using the random-cycle model requires three changes.  First, one needs to
conduct simulations using the Bose-gas interaction (equation
\ref{eqn:results_Bose_H}) rather than the cycle-weight interaction (equation
\ref{eqn:results_RCM_H}).  The interaction term $V$ is a CPU-intensive
Brownian-bridge computation \cite{bib:bu07}; unpublished work of Ueltschi and
Betz shows that it may be approximated in the weak-interaction case by a
simpler Riemann integral.  Precomputed tables and interpolation may make use of
this integral feasible.

Second, point positions must be allowed to vary on the continuum.  This entails
a second type of Metropolis step, in addition to that shown in section
\ref{subsec:SO_algorithm}.  Namely, one picks a point and moves it to a new
position nearby, using the detailed-balance condition to choose the acceptance
probability \cite{bib:kerl}.

Third, software efficiency requires a hierarchical partitioning of $\Lambda$.
The Metropolis step of section \ref{subsec:SO_algorithm} relies on picking
$\pi(\vecy)$ near to $\pi(\vecx)$.  For points on the lattice, this is easy:
each site has six nearest neighbors.  For freely placed points, though, one
must remember which sites are close to which.  The most naive implementation
involves computing the distances between all $N(N-1)/2$ pairs of points; the
$O(N^2)$ computation time is overwhelming.  Instead, the lattice may be
partitioned into subcubes.  Distances need to be computed only between each
given point $\vecx$ and those in $\vecx$'s subcube and the 26 nearest-neighbor
subcubes.

The second and third points simply require a software effort.  Implementing
them will be worthwhile only if the interaction terms can be simplified to the
point that they are computationally feasible.  This is a mathematical effort.

%% ================================================================
\section{Acknowledgements}
\label{sec:acknowledgements}

The author's doctoral dissertation work was done with co-advisors Daniel
Ueltschi and Tom Kennedy.  Daniel Gandolfo contributed several helpful
discussions.  The author was supported for three semesters by the National
Science Foundation, through NSF grant DMS-0601075 as well as the University of
Arizona Department of Mathematics VIGRE grant.
