%% ----------------------------------------------------------------
\section{Acknowledgements}

We implement Monte Carlo techniques for the following:

\begin{itemize}

\item \texttt{arXiv:cond-mat/0703315} (Gandolfo, Ruiz, Ueltschi): describes
the non-interacting model.  Referred to herein as the \emph{GRU paper}.
(Monte Carlo results were obtained by Gandolfo and Ruiz; we have reproduced
their results.)

\item \texttt{arXiv:0711.1188} (Betz, Ueltschi):  Describes the interacting
model in detail.  The U07 paper (next) summarizes much of the content of this
longer paper.

\item \texttt{arXiv:0712.2443v3} (Ueltschi):  Describes the interacting model.
Referred to herein as the \emph{U07 paper}.

\end{itemize}

This work (spring 2008) is supported by NSF grant DMS-0601075.

Many thanks to the following people for multiple insights:
Daniel Ueltschi,
Tom Kennedy,
Janek Wehr,
Ben Dyhr.

%% ================================================================
\newpage
\section{Review of spatial random permutations}

%% ----------------------------------------------------------------
\subsection{Hamiltonians for bosons and permutations}

As described in \cite{bib:bu} and \cite{bib:qmath},
the Hamiltonian for $N$ interacting bosons in a domain $\Lambda$ is
$$
	\bsH = -\sum_{i=1}^N \Delta_i
	+ \sum_{i<j} U(x_i-x_j) \qquad \text{in } L^2_{\rm sym}(\Lambda^N).
$$
Then, with inverse temperature $\beta$,
\begin{align*}
	 \Tr \e^{-\beta \bsH} =
	\sum_\pi \frac1{N!}
		\int \dd x_1 \dots \dd x_N &
		\int \dd W^{2\beta}_{x_1 x_{\pi(1)}}(w_1) 
		\dots
		\dd W^{2\beta}_{x_N x_{\pi(N)}}(w_N) \\
	&\exp\Bigl\{
		-\tfrac12
		\sum_{i<j}
		\int_0^{2\beta} U
		\bigl( w_i(s) - w_j(s) \bigr) \dd s
	\Bigr\}
\end{align*}
where $w_i(s)$ is a Brownian bridge running from $x_i$ to $x_{\pi(i)}$ in time
$2\beta$.  Write this as
$$
	\Tr \e^{-\beta \bsH}
	= \frac1{N!}
	\int_{\Lambda^N} \dd\bsx
	\sum_\pi \e^{-H(\bsx,\pi)}
$$
where
$$
\e^{-H(\bsx,\pi)} = \Bigl[
		\prod_{i=1}^N \dd W_{x_i x_{\pi(i)}}^{2\beta}(\omega_i)
	\Bigr]
	\exp\Bigl\{
		-\tfrac12
		\sum_{i<j}
		\int_0^{2\beta}
		U \bigl( w_i(s) - w_j(s) \bigr) \dd s
	\Bigr\}.
$$
After cluster expansion (a highly non-trivial step, as yet lacking rigrous
justification), one obtains
$$
H(\bsx,\pi) = \frac1{4\beta}
	\sum_{i=1}^N |x_i-x_{\pi(i)}|^2
	+ \sum_{i<j} V(x_i, x_{\pi(i)}, x_j, x_{\pi(j)}) + \text{higher orders}.
$$
Note that the bosonic Hamiltonian $\bsH$ has been converted to a Hamiltonian
$H$ on permutations.

The interaction between jumps $x \mapsto y$ and $x' \mapsto y'$ is
$$
	V(x,y,x',y') = \int \bigl[
		1 - \e^{-\frac14 \int_0^{4\beta} U(\omega(s)) \dd s}
	\bigr]
	\dd\widehat W_{x-x',y-y'}^{4\beta}(\omega).
$$

If $U$ is a hard-core potential with radius $a$ (i.e. $U(r)=\infty$ for $r < a$
and $U(r)=0$ for $r \ge a$), then $V(\cdot)$ is the probability that a Brownian
bridge from $x-x'$ to $y-y'$ hits the ball of radius $a$ centered at the
origin.

Is there a simple expression involving special functions?  Apparently not.

%% ----------------------------------------------------------------
\subsection{Models}

We simulate three models for spatial random permutations.  The first two have
been completely coded; the third is in progress.

\begin{itemize}

\item The \emph{non-interacting model} (\cite{bib:gru}):
	$$
		H(\bsx,\pi) = \frac1{4\beta}
			\sum_{i=1}^N |x_i-x_{\pi(i)}|^2.
	$$

\item The \emph{$r_2$ interacting model} (\cite{bib:qmath}):
	$$
		H(\bsx,\pi) = \frac1{4\beta}
			\sum_{i=1}^N |x_i-x_{\pi(i)}|^2
			+ \alpha r_2(\pi).
	$$

\item The \emph{interacting model} (\cite{bib:qmath}):
	$$
		H(\bsx,\pi) = \frac1{4\beta}
			\sum_{i=1}^N |x_i-x_{\pi(i)}|^2
			+ \sum_{i<j} V(x_i, x_{\pi(i)}, x_j, x_{\pi(j)}).
	$$

\end{itemize}

%% ----------------------------------------------------------------
\subsection{Conceptualization}

\begin{itemize}

\item The distance-dependent term $\e^{-\frac{1}{4\beta} \sum_x
\|x-\pi(x)\|^2}$ makes a permutation $\pi$ with a long jump (i.e. $\pi(x)$ far
from $x$) less probable.

\item The $e^{-\alpha r_2(\pi)}$ term discourages permutations with 2-cycles.

\item The interacting term discourages permutations with $x_i$ close to $x_j$
and $\pi(x_i)$ close to $\pi(x_j)$, regardless of jump lengths
$\|x_i-\pi(x_i)\|$ or $\|x_j-\pi(x_j)\|$.  The permutation is favored even less
if the two black arrows cross (i.e.  larger $\theta$ as discussed below).

\end{itemize}

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

%% ----------------------------------------------------------------
\subsection{Context}

The critical temperature $T_c$ for Bose-Einstein condensation is a (mostly
unknown) function of scattering length $a$.  Even the sign of the slope of
$T_c(a)$ near zero is contested.  It is believed that
$$
	\frac{T_c(a) - T_c(0)}{T_c(0)} = c \rho^{1/d} a + o(\rho^{1/3}a).
$$
Currently, it is thought that $c \approx 1.3$.  The Monte Carlo simulations
described here will permit tighter estimation of $c$.

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

%% ----------------------------------------------------------------
\subsection{Physics literature}

\begin{itemize}

\item 1964:
\emph{Huang}:
$\frac{\Delta T}{T_{\rm c}} \sim (a \rho^{1/3})^{3/2}$, increases

\item 1971:
\emph{Fetter \& Walecka}:
$\frac{\Delta T}{T_{\rm c}}$ decreases

\item 1982:
\emph{Toyoda}:
$\frac{\Delta T}{T_{\rm c}}$ decreases

\item 1992:
\emph{Stoof}:
$$
\frac{\Delta T}{T_{\rm c}} = c \, a \rho^{1/3} + o(a \rho^{1/3}), \quad c>0
$$

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

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

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

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

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

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

\item 2001:
\emph{Arnold, Moore}:
$c=1.32$

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

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

\end{itemize}

%% ----------------------------------------------------------------
\subsection{Critical temperature}

We define $\phi(\alpha)$ to be the probability that the origin is in an
infinite cycle.  (Here, $\alpha=1/4\beta$; this figure is from \cite{bib:gru}.)
At the critical temperature $\alpha$, $\phi(\alpha)$ goes to zero.  Monte Carlo
simulations undertaken in this project will discover how this graph changes in
the presence of interactions.

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

%% ================================================================
\newpage
\section{The computational project}

%% ----------------------------------------------------------------
\subsection{Density of sites in infinite cycles}

Given a random variable $\theta(\pi)$, compute its expected value.  The random
variable of interest for this project is the density of sites in cycles of
specified length:
$$
	\rho_{mn}(\pi) =
	\frac1V \#\bigl\{ i = 1,\dots,N : m \leq \ell_i(\pi) \leq n \bigr\}
$$
% xxx figure here with $\ell_x(\pi)$'s.

The usual prescription in probability is
$$
	E[\rho_{mn}]
	= \sum_{\pi\in \caS_N} \rho_{mn}(\pi) P(\pi)
	= \sum_{\pi\in \caS_N} \rho_{mn}(\pi) \frac{e^{-H(\bsx,\pi)}}{Y}.
$$

The computational burden splits into three main components:
\begin{itemize}
\item[(1)] Finding $H$, especially its $V$ term.  (For Metropolis,
	$\Delta H$ including $\Delta V$.)
\item[(2)] Sampling (via Metropolis) from a non-uniform probability
	distribution on $N!$ permutations for $N$ as big as $50^3$.
	%It turns out that for Metropolis, we don't need to compute $Y$:
	%the theory involves ratios of probabilities and the $Y$'s cancel.
\item[(3)] Visualizing the results.
\end{itemize}


%% ----------------------------------------------------------------
\subsection{Tools}

\begin{itemize}
\item Linux environment, although in principle everything should be portable
	to other operating systems.
\item Optimizing compiler: \texttt{gcc -O3}.
\item Build tool: \texttt{make} and automatic makefile generation.
\item Performance analyzer: \texttt{gprof}.  This shows where a program
	is spending most of its time.
\item Error detector: \texttt{valgrind}.  Finds many (but not all!)
	common errors, e.g. \texttt{malloc} without \texttt{free}.
\item Code navigation:  \texttt{ctags}.  Allows a smart editor
	(\texttt{vim}, \texttt{emacs}) to jump directly to a subroutine body.
\item Graphing utility:  \texttt{xgr}.  Nice plots in these slides
	were made in Matlab; quick-and-dirty plots without axis labels were done
	using \texttt{xgr}.
\end{itemize}

Sample \texttt{gprof} output:
\begin{verbatim}
%     cumul.  self
time  seconds seconds calls    name
64.74 1.67    1.67     3729508 pmt_send_x_to_y_n2_delta
10.08 1.93    0.26    13047273 get_distance_squared
 7.37 2.12    0.19     3617278 get_Delta_H
 5.81 2.27    0.15     3617278 x_to_uniform_y
 4.26 2.38    0.11     3617278 metro_step
 3.49 2.47    0.09     3617278 pmt_send_x_to_y
 3.10 2.55    0.08       10545 pmt_get_cycle_counts
 0.78 2.57    0.02       10000 vector_accumulate
 0.39 2.58    0.01       10546 metro_sweep
 0.00 2.58    0.00       10545 get_rho_L_pi
  ...
\end{verbatim}

%% ================================================================
\newpage
\section{Visualization}

There are three main plots:

\begin{itemize}
\item[(1)] Plots of the system energy $H$.
\item[(2)] Dot plots of the cycles.
\item[(3)] $E[\rho_{0,k}]$ as a function of $k$ from $0$ to $N$.
\end{itemize}

%% ----------------------------------------------------------------
\subsection{$H$ plots}

Here is a plot of system energy $H$ for $L=10$, $d=3$, no interactions, and
$\beta=1$:

\begin{center}
\includegraphics[scale=0.75]{figures/H_plot_beta_1.eps}
\end{center}

This plot is typical for various parameter values; only one such plot is shown
here.

\begin{itemize}
\item The horizontal axis counts Metropolis sweeps.
\item The system was found to be thermalized (as described below) after 559
	steps; $\rho$ values were accumulated over 10,000 sweeps.
\item The system energy $H$ is shown in blue.
\item In red is $H$ smoothed out over a sliding window of 100 sweeps.
\item In green is the same smoothed system energy, multiplied by 0 before
	thermalization and 1 after.  Thus, the plot ``goes green'' when
	thermalization has occurred.
\end{itemize}

%% ----------------------------------------------------------------
\subsection{Dot plots}

A \emph{dot plot} of the points $\{x_1, \ldots, x_N\}$ and a permutation $\pi$
has a dot for each point $x$, along with a line from $x$ to $\pi(x)$ for each
point $x$.

% xxx note these are qualitative only -- not quantitative.  For that, we
% use the rho plots.

Key points:

\begin{itemize}

\item For infinite $\beta$, the permutation weight
$e^{-\frac{1}{4\beta}\sum_{x} \|x-\pi(x)\|^2}$ becomes uniform: individual
permutation jumps can be arbitrarily long.

\item For $\beta=0$, only the identity permutation is possible.

\item For moderate $\beta$, long jumps are discouraged.  Nonetheless, a long
cycle can occur when short jumps chain together.

\end{itemize}

%% ----------------------------------------------------------------
\subsection{Plots for the non-interacting model}

Here is $L=10$, $d=3$, point positions uniformly distributed on the cube
of width 10 but not metropolized, no interactions, varying $\beta$:

% xxx annotated plot:  for identity and uniform; arrow for decreasing $T$
% xxx show T_c and \phi(\beta).
% xxx replace GRU alpha with beta throughout.

\begin{center}
\includegraphics[scale=0.8]{figures/dot_plot_beta_0.001.eps}
\includegraphics[scale=0.8]{figures/dot_plot_beta_0.01.eps}

\includegraphics[scale=0.8]{figures/dot_plot_beta_0.05.eps}
\includegraphics[scale=0.8]{figures/dot_plot_beta_0.5.eps}
\includegraphics[scale=0.8]{figures/dot_plot_beta_10.eps}
\end{center}

The $E[\rho]$ plots are much as in the GRU paper.

\begin{itemize}

\item The horizontal axis is $k/N$ for $k$ from 0 to $N$.

\item In blue on the vertical axis is $\rho_{0,k}$ for the permutation
realized on the last Metropolis sweep.

\item In green on the vertical axis is $\rho_{k,k}$ for the permutation
realized on the last Metropolis sweep.

\item In red on the vertical axis is $E[\rho_{0,k}]$ over 10,000 Metropolis
sweeps.

\item In yellow on the vertical axis is $E[\rho_{k,k}]$.

\end{itemize}

Here are $E[\rho]$ plots for the same parameter values as the dot plots:

\begin{center}
\includegraphics[scale=0.8]{figures/rho_plot_beta_0.001.eps}
\includegraphics[scale=0.8]{figures/rho_plot_beta_0.01.eps}

\includegraphics[scale=0.8]{figures/rho_plot_beta_0.05.eps}
\includegraphics[scale=0.8]{figures/rho_plot_beta_0.5.eps}
\includegraphics[scale=0.8]{figures/rho_plot_beta_10.eps}
\end{center}

%% ----------------------------------------------------------------
\subsection{Plots for the $r_2$ model}

Here we fix $\beta=0.5$ and vary $\alpha$.  Note that $\alpha=0$ recovers the
non-interacting case.  The dot plots are indistinguishable from the
non-interacting case.  The $E[\rho_{0,k}]$ plots are similar, so they are
superimposed.  Blue is $\alpha=0$, red is $\alpha=5$, and green is $\alpha=20$.

\begin{center}
\includegraphics[scale=0.8]{figures/rho_plot_beta_0.5_alpha_0_5_20.eps}
\end{center}

%% ----------------------------------------------------------------
\subsection{Plots for the interacting model}

This is recent work --- more are to be obtained.

Here is $\beta=0.15626$ (just below non-interacting critical temperature), with
$a=0.0$ and $a=0.1$:

\begin{center}
\includegraphics[scale=0.8]{figures/rho_plot_gru_alpha_1.6_a_0.0.eps}
\includegraphics[scale=0.8]{figures/rho_plot_gru_alpha_1.6_a_0.1.eps}
\end{center}

The value $\phi(\beta)$ is the probability that the origin is in an
``infinite'' cycle.  It may be read off the $E[\rho]$ plots as the distance
from the upper left corner of the $\rho$ plot to the first leftward lean of the
red curve.  Critical $\beta_c$ has $\phi(\beta)=0$.

\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline $\beta$ & $\phi_0(\beta)$ & $\phi_{\alpha=4}(\beta)$ &
	$\phi_{a=0.1}(\beta)$ \\
\hline
\hline 0.227273 & 0.5203 & 0.5824 & 0.8081 \\
\hline 0.208333 & 0.4373 & 0.5114 & 0.8057 \\
\hline 0.192308 & 0.3703 & 0.4440 & 0.7835 \\
\hline 0.178571 & 0.2625 & 0.3097 & 0.7868 \\
\hline 0.166667 & 0.1517 & 0.2148 & 0.7769 \\
\hline 0.161290 & 0.1133 & 0.1637 & 0.7663 \\
\hline 0.156250 & 0.0824 & 0.1220 & 0.7693 \\
\hline 0.147059 & 0.0311 & 0.0351 & 0.7645 \\
\hline
\end{tabular}
\end{center}

Conclusion:  interactions lower critical $\beta$.  More simulations are needed.

%% ================================================================
\newpage
\section{Computation of $V$}

%% ----------------------------------------------------------------
\subsection{Brownian bridges}

Write $\xhat$ and $\yhat$ for $x-x'$ and $y-y'$ respectively.  Simply generate
$N_b$ Brownian bridges from $\xhat$ to $\yhat$, with $N_p$ mesh points per
bridge, and see what fraction of them intersects the ball of radius $a$
centered at the origin.

\begin{itemize}
\item Start with a unit-uniform pseudorandom number generator (RNG).
\item Use a Box-Muller transform (cf. \emph{Numerical Recipes}) to
	get standard-normal deviates.
\item Brownian motion for $t$ from 0 to 1 in steps of $\Delta t$:
	$B_0=0$ and $B_{t+1} = B_t + \Delta B$ where $\Delta B$ is normal
	with mean zero and variance $\Delta t$.
\item Brownian bridge from $\xhat=0$ to $\yhat=0$ for $t$ from 0 to 1:
	$R_t = B_t - t B_1$.
\item Brownian bridge from $\xhat$ to $\yhat$ for $t$ from 0 to $T$:
	$
		\sqrt{T} R_t + \xhat + \frac{t}{T}(\yhat-\xhat).
	$
\end{itemize}

The plot on the left shows, for $d=1$, $N_b=20$ bridges with $N_p=1000$ points
per bridge, bridged from $x=-1$ to $y=2$ ($d=1$) in time $T=1$, with $R_t$
plotted against $t$.

The plot on the right shows, for $d=3$, the trajectory of a single bridge from
$x=(-1,0,0)$ to $y=(2,0,0)$ in time $T=1$, with the first two components of
$R_t$ plotted.

\begin{center}
\psfragscanon
\includegraphics[scale=0.8]{figures/bridges_vs_t.eps}
\includegraphics[scale=0.8]{figures/bridge.eps}
\end{center}

%% ----------------------------------------------------------------
\subsection{Software testing}

Incrementally test the subroutines for Brownian motion,
zero-to-zero bridges, and general bridges.

\begin{itemize}
\item Generate $N_b$ bridges of $N_p$ points each, for $d=1,2,3$.
\item Select time slices $s$ and $t$.
\item Compute sample means, sample variances, and sample covariances for
	those time slices and compare against theoretically expected results.
\end{itemize}

Expected results for Brownian motion:
$$
	E[B_t] = 0,
	\quad
	\Var[B_t] = t,
	\quad
	\Cov[B_s, B_t] = s \wedge t.
$$

Expected results for zero-to-zero Brownian bridge ($s<t$ to simplify notation):
$$
	E[R_t] = 0,
	\quad
	\Var[R_t] = t(1-t),
	\quad
	\Cov[R_s, R_t] = s (1-t).
$$

Expected results for $\xhat$-to-$\yhat$ Brownian bridge ($s<t$):
$$
	E[R_t] = \xhat + \frac{t}{T}(\yhat-\xhat),
	\quad
	\Var[R_t] = \frac{t(T-t)}{T},
	\quad
	\Cov[R_t] = \frac{s(T-t)}{T}.
$$

Example with $T=$, $\xhat=(2,0,0)$, $\yhat=(-2,0,0)$, $N_b=1000$, $N_p=1000$,
$s=0.004$ (i.e. index 1 of 1000), $t=2.0$ (i.e. index 500 of 1000):

\begin{verbatim}

Actual     E[b(s)]        =  1.9976  -0.0024 -0.0014
Expected   E[b(s)]        =  1.9960  -0.0000  0.0000
Difference E[b(s)]        =  0.0016  -0.0024 -0.0014

Actual     E[b(t)]        =  0.0086   0.0542 -0.0208
Expected   E[b(t)]        =  0.0000  -0.0000  0.0000
Difference E[b(t)]        =  0.0086   0.0542 -0.0208

\end{verbatim}

\begin{verbatim}
...

Actual     Var[b(s)]      =  0.0037   0.0040  0.0040
Expected   Var[b(s)]      =  0.0040   0.0040  0.0040
Difference Var[b(s)]      = -0.0003   0.0000  0.0000

Actual     Var[b(t)]      =  1.0582  0.9712  0.9175
Expected   Var[b(t)]      =  1.0000  1.0000  1.0000
Difference Var[b(t)]      =  0.0582 -0.0288 -0.0825

Actual     Cov[b(s),b(t)] =  0.0021  0.0020  0.0027
Expected   Cov[b(s),b(t)] =  0.0020  0.0020  0.0020
Difference Cov[b(s),b(t)] =  0.0001  0.0000  0.0007

\end{verbatim}

%% ----------------------------------------------------------------
\subsection{Bridge results}

Experimental results are discouraging.  Performance requirements are too
stiff for generation of Brownian bridges during Metropolis steps.
To help this, one can (1) compute a database of zero-to-zero $N_b$ Brownian
bridges of $N_p$ points each, and re-use this database for different
$\xhat,\yhat$.  (2) Tabulate $V$ off-line and interpolate at runtime.

\begin{itemize}

\item Dependence on $N_b$: Increasing $N_b$ decreases sampling variability of
$V$.

\item Dependence on $N_p$: For small $N_p$, increasing $N_b$ only decreases
sampling variablity, but non-zero bias remains (vs. the integral and exact
expressions, shown next).  For the test case $r_1=1$, $r_2=1$, $\theta=\pi$,
one needs $N_p$ on the order of 500,000 before $V$ begins to stabilize.

\end{itemize}

Interpretation:  Note that $\Delta t = T/N_p$.  Standard deviation of bridge
steps is on the order of $\sqrt{T/N_p}$.  For smaller $N_p$, bridges are too
``hoppy'' and miss the $a$-ball at the origin.


%% ----------------------------------------------------------------
\subsection{Integral expression for $V$}

Ueltschi and Betz have recently found an approximation which is valid to low
order in $a$:
\begin{align*}
	V_2(\xhat,\yhat) &= \frac{a}{\sqrt{8\pi\beta}}
		e^{\frac{+\|\xhat-\yhat\|^2}{8\beta}}
		\int_0^1 \frac{1}{[s(1-s)]^{3/2}}
			e^{-\frac{\|\xhat\|^2}{8\beta s}}
			e^{-\frac{\|\yhat\|^2}{8\beta (1-s)}}
			\dd s.
\end{align*}
where, for notational convenience, we write
$$
	\xhat = x-x',
	\quad
	\yhat = y-y',
	\quad
	V_2(\xhat,\yhat) = V(x,y,x',y').
$$
If $\|\xhat\|=\|\yhat\|$ then we have the exact expression
\begin{align*}
	V_2(\xhat,\yhat) &= \frac{2a}{\|\xhat\|}
		e^{\frac{+\|\xhat-\yhat\|^2}{8\beta}}
		e^{\frac{-\|\xhat\|^2}{2\beta}}.
\end{align*}
This can be written in terms of the five real variables
$r_1=\|x\|$,
$r_2=\|y\|$,
$\theta = \cos^{-1}(\langle x, y \rangle/\|x\|\,\|y\|)$,
$\beta$, and
$a$.

%% ----------------------------------------------------------------
\subsection{Argument reduction}

The potential $V$ depends on $d^4$ real variables:  we have $V(x,y,x',y')$
where $x,x',y,y' \in \R^d$.  Since only $x-x'$ and $y-y'$ appear in the
formula, we can reduce to $d^2$ real variables: we have $V_2(\xhat,\yhat)$ as
above.

Using rotation and translation invariance of $V$, we can write down $V$ in
terms of $r_1=\|\xhat\|$, $r_2=\|\yhat\|$, and angle $\theta$.  Using the Law
of Cosines, we have
$$
	\|\xhat-\yhat\|^2 = r_1^2 + r_2^2 - 2 r_1 r_2 \cos(\theta).
$$
In particular, when $r_1=r_2=r$, we have
$$
	\|\xhat-\yhat\|^2
	= 2 r^2 (1 - \cos(\theta))
	= 4 r^2 \sin^2(\theta/2).
$$

Now $V$ depends only on three real variables.  The integral expression is
\begin{align*}
	V_2(r_1, r_2, \theta) = \frac{a}{\sqrt{8\pi\beta}}
		e^{\frac{r_1^2+r_2^2-2r_1r_2\cos\theta}{8\beta}}
		\int_0^1 \frac{1}{[s(1-s)]^{3/2}}
			e^{-\frac{r_1^2}{8\beta s}}
			e^{-\frac{r_2^2}{8\beta (1-s)}}
			\dd s
\end{align*}
and the exact expression, for $r_1=r_2=r$, is
\begin{align*}
	V_2(r,r,\theta) = \frac{2a}{r}
	e^{\frac{-r^2(1-\sin^2(\theta/2))}{2\beta}}.
\end{align*}

%% ----------------------------------------------------------------
\subsection{Visualization}

Here is a surface plot of $V(r,r,\theta)$ for $r$ from 1 to 4, $\theta$ from 0
to $\pi$, $\beta=1$, and $a=0.1$.  Note that probability of intersecting the
$a$-sphere decays as $r$ increases, and grows as $\theta$ runs from
$0^\circ$ to $180^\circ$, as expected.
\begin{center}
\includegraphics[scale=0.8]{figures/vrrt.eps}
\end{center}

%% ================================================================
\newpage
\section{Metropolis-Hastings}
\label{sec:rcm_metro}

%% ----------------------------------------------------------------
\subsection{Overview}
The \emph{Metropolis-Hastings algorithm} is a special case of Monte Carlo
Markov chain (MCMC).  It is best introduced by example: consider the 1D
$N$-point Ising model.

\begin{itemize}

\item One has a \emph{system} with multiple possible \emph{configurations}.
In the Ising model, 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}

\item A state is described by
$
	\vecomega = (\omega_1, \ldots, \omega_n).
$
The configuration space $\Omega$ has $2^N$ possible configurations.

\item The system is endowed with an \emph{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.

\item One picks an \emph{initial configuration}.  Typically, there are three
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$.

\item There is a \emph{temperature}-related parameter $\beta$.

\item One selects a \emph{site} $i$ and decides whether to flip $\omega_i$ to
$-\omega_i$.

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

\item This decision is made using the \emph{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 \emph{Metropolis step}\index{step}.

\item Looping through all $n$ sites from $i=1$ to $i=n$, performing a
Metropolis step at each site $i$, is called a \emph{Metropolis
sweep}\index{sweep}.

\item If one realizes a random variable $\theta(\pi)$ at each of $M$ sweeps,
averaging $\theta$ over the $M$ sweeps, one obtains an approximation
$\overline\theta$ for the expectation $E[\theta]$.

\item 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 \emph{thermalization
phase}; the $M$ sweeps are called the \emph{accumulation phase}.

There is no general method to determine whether the system has thermalized
(Kennedy); the underlying concern is the convergence of the Metropolis
probability distribution to the stable distribution of its implicit Markov
chain.

\end{itemize}

%% ----------------------------------------------------------------
\subsection{Metropolis for the random-cycle model}

The random-cycle model is metropolized in a manner analogous to the 1D Ising
model:

\begin{itemize}

\item The \emph{state space} is $\caS_N$, the permutations on $N$ elements.
It has size $N!$.

\item The \emph{energy function} is $H$ as described above.

\item The \emph{initial configuration} is found in one of two ways: (1) Start
with an identity permutation.  (2) Start with a permutation $\pi$ selected from
a uniform probability distribution on $\caS_N$.

\item There is a \emph{temperature}-related parameter $\beta$.

\item The analogue of conditionally flipping one of $N$ Ising spins is the
following.  One selects two of the $N$ points $x$ and $y$ and decides whether
to send the old permutation $\pi$ to new permutation $\pi'$ via
$$
	\pi:
	\left(
	\begin{matrix}
		x & \iy & u & \cdots\\
		\downarrow & \downarrow & \downarrow\\
		\ix & y & \iu & \cdots
	\end{matrix}
	\right)
	\qquad
	\pi':
	\left(
	\begin{matrix}
		x & \iy & u & \cdots\\
		\downarrow & \downarrow & \downarrow\\
		y & \ix & \iu & \cdots
	\end{matrix}
	\right)
$$
where $\pi'(u) = \pi(u)$ for all $u \ne x, \iy$.

\item The energy change $\Delta H$ may be expensively computed by finding
$H(\pi')$ and $H(\pi)$ and subtracting the two, but again, one may do some
algebra to take advantage of the fact that most of the terms are identical.
This computation is shown below.

\item A single \emph{Metropolis step} selects a site $x$ uniformly from the
lattice.  The site $y$ is selected uniformly from a \emph{Metropolis window}:
consider only sites $y$ within the radius $r$ such that
$e^{-\frac{1}{4\beta}r^2} \approx e^{-10}$, i.e.  $r \approx 6 \sqrt{\beta}$.
(Selecting $y$ uniformly from the entire lattice gives too many rejected
Metropolis steps, harming performance.)

\item Looping through all $N$ sites $x$, performing a Metropolis step at each
site $i$, is a \emph{Metropolis sweep}\index{sweep}.

\item The principal random variable of interest is $\rho_{mn}$ as described
above.

\item Thermalization is detected as follows:  smooth out $H$ over a sliding
window of 100 Metropolis sweeps.  Consider the system thermalized when this
smoothed $H$ has reached 30 turning points.  Heuristically, this is overly
conservative (which is fine).  It is better to run too many thermalization
sweeps than too few.

\end{itemize}

%% ----------------------------------------------------------------
\subsection{$\Delta H$ for the non-interacting model}

Recall that
\begin{align*}
	H(\pi) &= \frac{1}{4\beta} \sum_{u \in \Lambda} \|u-\pi(u)\|^2.
\end{align*}

For an implementation of a Metropolis algorithm, one wishes to compute the
change in potential energy when one sends $\pi$ to $\pi'$ via
$$
	\pi:
	\left(
	\begin{matrix}
		x & \iy & u & \cdots\\
		\downarrow & \downarrow & \downarrow\\
		\ix & y & \iu & \cdots
	\end{matrix}
	\right)
	\qquad
	\pi':
	\left(
	\begin{matrix}
		x & \iy & u & \cdots\\
		\downarrow & \downarrow & \downarrow\\
		y & \ix & \iu & \cdots
	\end{matrix}
	\right)
$$
where $\pi'(u) = \pi(u)$ for all $u \ne x, \iy$.

We have
\begin{align*}
	H_{\pi'} - H_\pi
	&=	\frac{1}{4\beta}
	\sum_{u \in \Lambda} \bigg(
		\|u-\pi'(u)\|^2 -
		\|u-\pi (u)\|^2
	\bigg).
\end{align*}
Since $\pi'(u) = \pi(u)$ for all $u \ne x, \iy$, we may rewrite this in terms
of $\pi$ only as
$$
	\Delta H =
	\frac{
	\|x   - y \|^2 - \|x    - \ix\|^2
	+\|\iy - \ix\|^2  - \|\iy - y  \|^2}
	{4\beta}.
$$


%% ----------------------------------------------------------------
\subsection{$\Delta H$ for the $r_2$ model}
The distance-dependent term is the same as in the non-interacting case.
Additionally, we need to compute $\Delta r_2 = r_2(\pi') - r_2(\pi)$.  There
are three cases:

\begin{itemize}
\item[(1)] $x \ne y$ and $x,y$ are in different cycles;
\item[(2)] $x \ne y$ and $x,y$ are in the same cycle;
\item[(3)] $x = y$.
\end{itemize}

In each case, it is easy to track the change $\Delta r_2$, without having to
compute $r_2(\pi)$ and $r_2(\pi')$.


%
%%%% ----------------------------------------------------------------
%%\frametitle{\sectitle: $\Delta H$ for the $r_2$ model}
%
%%  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%% Case (1): x != y and x,y are in different cycles.
%% The distinct cycles are merged.
%% * lx  = length of (old) cycle containing x.
%% * ly  = length of (old) cycle containing y.
%% * lxy = length of (new) cycle containing x and y.
%%   We have lxy = lx + ly.
%% * Updates:
%%   if (lx  == 2) n2--
%%   if (ly  == 2) n2--
%%   if (lxy == 2) n2++
%%
%
%%   Before:
%%          o-<--o                     o-<--o
%%         /      \                   /      \.
%%        /        \                 /        \.
%%       .          x               y          .
%%       .          |               ^          .
%%       .          |               |          .
%%       .          |               |          .
%%       .          v               |          .
%%       .        pi(x)          piinv(y)      .
%%        \        /                 \        /
%%         \      /                   \      /
%%          o-->-o                     o-->-o
%%
%%   After:
%%          o-<--o                     o-<--o
%%         /      \                   /      \.
%%        /        \                 /        \.
%%       .          x ------------> y          .
%%       .                                     .
%%       .                                     .
%%       .                                     .
%%       .                                     .
%%       .        pi(x) <------- piinv(y)      .
%%        \        /                 \        /
%%         \      /                   \      /
%%          o-->-o                     o-->-o
%
%%  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%% Case (2): x != y and x,y are in the same cycle.
%% The cycle is split.
%% * lxy  = length of (old) cycle containing x and y.
%% * lx   = length of (new) cycle containing y.
%% * lpix = length of (new) cycle containing pi(x).
%% * Updates:
%%   if (lxy  == 2) n2--
%%   if (lx   == 2) n2++
%%   if (lpix == 2) n2++
%%
%%   Before:
%%                     o-< ... -<o
%%                    /           \.
%%                   /             \.
%%                  x               y
%%                  |               ^
%%                  |               |
%%                  |               |
%%                  v               |
%%                pi(x)          piinv(y)
%%                   \             /
%%                    \           /
%%                     o-> ... ->o
%%
%%
%%   After:
%%                     o-< ... -<o
%%                    /           \.
%%                   /             \.
%%                  x ------------> y
%%
%%
%%
%%                pi(x) <------- piinv(y)
%%                   \             /
%%                    \           /
%%                     o-> ... ->o
%%
%%
%%  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%% Case (3): x == y.
%% The singleton x is excised from the cycle.
%% * lxy = length of (old) cycle containing x.
%% * ln  = length of (new) cycle excluding x.
%%   We have ln = lxy - 1.
%% * Updates:
%%   if (lxy == 2) n2--
%%   if (ln  == 2) n2++
%%
%%   Before:
%%                     o--<--o
%%                    /        \.
%%                   /          \.
%%               piinv(y)        .
%%                  |            .
%%                  |            .
%%                  v            .
%%                 x=y           .
%%                  |            .
%%                  |            .
%%                  |            .
%%                  v            .
%%                pi(x)          .
%%                   \          /
%%                    \        /
%%                     o-->--o
%%
%%   After:
%%                     o--<--o
%%                    /        \.
%%                   /          \.
%%               piinv(y)        .
%%                  |      __    .
%%                  |     /  \   .
%%                  |    |    |  .
%%                  |   x=y   ^  .
%%                  |    |    |  .
%%                  |     \__/   .
%%                  |   1-cycle  .
%%                  v            .
%%                pi(x)          .
%%                   \          /
%%                    \        /
%%                     o-->--o
%

%% ----------------------------------------------------------------
\subsection{$\Delta H$ for the interacting model}

The distance-dependent term is again the same as in the non-interacting case.
Additionally, we need to compute $\Delta V$.  Recall that
$$
	V(\pi) = \sum_{u<v} V(u, \pi(u), v, \pi(v))
$$
where $u$, $v$ are lattice sites.  When we send $\pi$ to $\pi'$ as above
($x \mapsto y$ replaces $x \mapsto \ix$ etc.), we have
\begin{align*}
	V_{\pi'} - V_\pi
	&=	\sum_{u<v} V(u, \pi'(u), v, \pi'(v))
	-	\sum_{u<v} V(u, \pi(u), v, \pi(v)).
\end{align*}

Since $\pi'(u) = \pi(u)$ for all $u \ne x, \iy$, and since
the interaction potential satisfies the symmetry condition
$$
	V(x,y,x',y') = V(x',y',x,y),
$$
we have simply
\begin{align*}
	\Delta V
&= \sum_{{u\ne x,} \atop{u\ne \iy}} \bigg(
		V(u,   \iu, x,   y)
	-	V(u,   \iu, x,   \ix) \bigg) \\
&+ \sum_{{u\ne x,} \atop{u\ne \iy}} \bigg(
		V(u,   \iu, \iy, \ix)
	-	V(u,   \iu, \iy, y) \bigg) \\
%
	&+	V(x,   y,    \iy,  \ix)
	-	V(x,  \ix, \iy,  y).
\end{align*}


%%% ================================================================
\section{Conclusions and further directions}

%% ----------------------------------------------------------------
\subsection{Conclusions}

\begin{itemize}

\item The $r_2$ model is easy to simulate.  The $r_2$ term raises the critical
temperature.  One can quantify this dependence and verify it against the result
of Betz and Ueltschi.

\item Preliminary results show that in the full-interaction model, the
critical temperature is also raised.  Software optimization is currently in
progress, so that more simulations may be done in a timely manner.  Then,
$T_c(a)$ may be plotted with confidence.

\end{itemize}

%% ----------------------------------------------------------------
\subsection{Further directions}

\begin{itemize}

\item The cluster expansion is non-rigorous and needs further justification,
in particular for non-lattice point distributions where inter-particle spacing
can be small.

\item Examine random variables other than $\rho_{mn}$.

\item Use non-Gaussian weights for $d=2$.

\item Place the points not on a cubic lattice but distributed according to a
point process; metropolize point positions as well as permutations.  The
correct point process for Bose-Einstein condensation is not known; it is known
\emph{not} to be Poisson.

\item We can greatly increase system size by using parallelization: on a
multiprocessor system, partition $\Lambda$ into subcubes.  When $x,y$ are in
the same subcube, computation is local; when $x$ is in one subcube and $y$ is
in a neighbor, use message-passing.

\item See what people come up with as $T_c(a)$ becomes better known \ldots .
Stay tuned for this as well!

\end{itemize}
