%% ================================================================
\section{Introduction}

My comprehensive examination exposits the paper \cite{bib:qmath}.  The problem
at hand is to determine the effects of interparticle interactions on the
critical temperature of Bose-Einstein condensation; a so-called random-cycle
model is employed to this end.

A large fraction of this paper is occupied by the appendices.  This reflects
the fact that much of my work on this project has been in filling the gap
between my coursework and the research papers being presented.  The experienced
reader may wish to skip the appendices entirely.

The plan of this paper is as follows:

\begin{itemize}

\item A non-technical discussion of the historical context of the project; a
mention of alternative approaches.

\item The random-cycle approach to the BEC problem requires a model of spatial
permutations.  This model is developed mathematically --- it is of
probabilistic interest on its own --- without reference to the Bose gas.  (Main
results about the latter are summarized in the appendices.)

\item The model of spatial permutations is connected to the Bose gas:

	\begin{itemize}

	\item One begins with a Hamiltonian for particles with two-body
	interactions.

	\item Using a multi-body Feynman-Kac approach involving permutation
	symmetry of bosonic wave functions, one obtains a Hamiltonian in which
	permutation jumps rather than particles interact.

	\item A cluster expansion, to first order in the scattering length of the
	particles, yields a Hamiltonian with only jump-pair interactions.  At this
	point, one may employ the random-cycle model of spatial permutations.

	\end{itemize}

\item A simplified random-cycle model, called the two-cycle-interaction model,
is described.  In particular, this simplified model is amenable to Monte Carlo
simulation and will play a key role in my dissertation work.  Brownian bridges
are the theoretical workhorse of this paper, but are computationally expensive
to simulate.  The two-cycle-interaction model bypasses Brownian bridges
entirely.

\end{itemize}

%% ================================================================
\newpage
\section{Historical context}
\label{sec:history}

%% ----------------------------------------------------------------
\subsection{Theory}

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

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

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

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

%% ----------------------------------------------------------------
\subsection{Experiments}

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

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

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

\begin{itemize}

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

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

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

\end{itemize}

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

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

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

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

Ueltschi makes the following summary of the theoretical work on this question.
(See also \cite{bib:arnoldmoore,bib:kps,bib:nholandau} for a review on the
widely varying experimental results on $T_c^{(a)}$; see \cite{bib:baymetal} for
a thorough listing of the progress up to 2001.)

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

The work of Ueltschi et al.  \cite{bib:gru,bib:bu,bib:ducycles,bib:qmath}
extends the permutation point of view originated by Feynman, Penrose, and
Onsager, drawing on the work of S\"ut\H{o}, Buffet, and Pul\`e
\cite{bib:feynman,bib:penons,bib:suto1,bib:suto2,bib:buffetpule}.  The main
goal of the project is to quantify $\Delta T/T_c$ for non-ideal Bose gases in
the small-scattering-length regime.  As is often the case in statistical
mechanics, the study of this interacting system necessitates the use of
computational methods.  The papers \cite{bib:bu,bib:qmath}, which are exposited
in my comprehensive examination, develop the theoretical basis for the model
which my disseration work will explore computationally and statistically.

% The main goal of this project is to determine
% the constant of linearity with a high degree of certainty.

% Order parameter:  Penrose-Onsager $\sigma_{\vecx\vecy}$ and Feynman's
% $\rho_\infty$.  xref to DU paper about the connection between the two order
% parameters.

%%% ----------------------------------------------------------------
%\subsection{Alternative approaches}
%
%The following works are cited by Ueltschi as alternative computational
%approaches to the interacting-Bose-gas problem:
%
%\begin{itemize}
%
%\item $O(2) \phi^4$ field theory: Arnold and Moore 2001 \cite{bib:arnoldmoore}.
%
%\item Worm algorithm: Kashurnikov et al. \cite{bib:kps}.
%
%\item Path integral\index{path integral Monte Carlo} Monte Carlo methods:
%Ceperly, Nho, and Landau \cite{bib:ceperley,bib:nholandau}.
%
%\end{itemize}
%
%I will have more to say about these papers in my dissertation; for the present,
%I can only draw the reader's attention to them.

%% ================================================================
\newpage
\section{Models of spatial permutations}
\label{sec:rcm}

Here we define and describe two configuration models of \plainidx{spatial
permutations} from a mathematical point of view.  These models will be related
to the physics of the Bose gas in section \ref{sec:connect}.  The exposition
here follows section 2 of \cite{bib:qmath}, with additional insights drawn from
section 2 of \cite{bib:bu}.

%% ----------------------------------------------------------------
\subsection{Definitions}
\label{sec:model_defns}

\textbf{State space:}
Let $\Lambda \subset \R^d$ be a cube of width $L$ and volume $V=L^d$.  Let $N
\in \Z^+$ and $\vecX = (\vecx_1, \ldots, \vecx_N)$ for $\vecx_1, \ldots,
\vecx_N \in \Lambda$.  The \plainidx{state space} of the model of spatial
permutations is
$$
	\Omega_{\Lambda,N} = \Lambda^N \times \SN
$$
where $\SN$ is the group of \plainidx{permutations} of $N$ points.  (See figure
\ref{fig:piham}.)

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

\textbf{Hamiltonian:}
The \plainidx{probability measure} on this state space will be constructed
through a Hamiltonian.  The background probability measure is discrete
(uniform) in $\pi$ and continuous (Lebesgue) in $\vecX$.  The Hamiltonian is as
follows:
\begin{align}
\label{eqn:rcm_H}
	\HP(\vecX,\pi) &= \sum_{i=1}^N \frac{1}{4\beta}\|\vecx_i-\vecxpii\|^2
		+ \sum_{1\le i<j \le N}
			V(\vecx_i, \vecxpii, \vecx_j, \vecxpij).
\end{align}
(The $1/4\beta$ factor\index{$\beta$} is the reciprocal of what is more
commonly encountered in statistical mechanics.  For now, think of this as an
\plainidx{ansatz}.  In section \ref{sec:U0}, this choice of scale factor will
be found to be appropriate.)

The operator $H$ is unbounded\index{unbounded operator}, but it is
symmetric\index{symmetric operator} so we consider its \plainidx{self-adjoint
extension}.  We take its domain to be $f$ in $C^2(\Lambda^N)$ w/
\plainidx{Dirichlet boundary conditions}.

There are two contributions to the \plainidx{energy} of a configuration
$(\vecX,\pi)$:

\begin{itemize}

\item The first\footnote{The papers \cite{bib:bu} and \cite{bib:qmath}
generalize from $\|\vecx_i-\vecxpii\|$ to $\xi(\vecx_i,\vecxpii)$
where $\xi$ is a spherically symmetric non-negative-valued function on $\R^d$
having integrable $e^{-\xi}$.  This generalization is not of interest in the
current paper, nor will it be of interest in my dissertation to follow.}
contribution to the energy is the sum of squares of permutation jump lengths.
(One jump is marked with a heavy arrow in figure \ref{fig:piham}.)  This makes
permutations with long jumps disfavored; permutations with many short jumps
will be less strongly disfavored.

\item The second contribution to the energy is the double sum over interactions
between \plainidx{permutation jumps}.  (One such \plainidx{jump-pair
interaction} is marked with wavy grey lines in figure \ref{fig:piham}.)

\end{itemize}

\textbf{Jump-interaction potentials:}
The specific form of the jump-pair interaction is left unspecified at present,
but we will require the following properties: $V$ is\index{translation
invariant} translation-invariant i.e.
$$
	V(\vecx,\vecy,\vecx',\vecy')
	= V(\vecx+\veca,\vecy+\veca,\vecx'+\veca,\vecy'+\veca)
$$
for all $\veca \in \Lambda$, and
$$
	V(\vecx,\vecy,\vecx',\vecy') = V(\vecx',\vecy',\vecx,\vecy)
$$
for all $\vecx,\vecy \in \Lambda$.

We will consider three different \plainidx{jump-interaction potentials}:
\begin{itemize}
\item One may certainly take the jump interaction $V$ to be zero.
\item In section \ref{sec:tji}, the jump interaction $V$ will involve the
probability of the intersection of \plainidx{Feynman paths}.
\item In section \ref{sec:n2}, we will treat all jump pairs as non-interacting
unless they participate in a \plainidx{two-cycle}.
\end{itemize}

\textbf{Partition functions:}
We consider two \index{partition function}partition functions, for a fixed
point configuration $\vecX$ and for an average over point configurations,
respectively:
$$
	Y(\Lambda,\vecX) = \sum_{\sigma\in\SN} e^{-\HP(\vecX,\sigma)}
	\qand
	Z(\Lambda,N) = \frac{1}{N!} \int_{\Lambda_N} Y(\Lambda,\vecX)\,d\vecX.
$$

\textbf{Probability measures:}
Then we have two corresponding discrete probability measures\index{probability
measure} on the finite set $\SN$, for a fixed point configuration $\vecX$ and
for an average over point configurations, respectively:
$$
	P_{\Lambda,\vecX}(\pi)
	= \frac{e^{-\HP(\vecX,\pi)}}{\sum_{\sigma\in\SN} e^{-\HP(\vecX,\sigma)}}
	\qand
	P_{\Lambda,N}(\pi)
	= \frac
		{\frac{1}{N!} \int_{\Lambda^N}d\vecX e^{-\HP(\vecX,\pi)}}
		{\frac{1}{N!} \int_{\Lambda^N}d\vecX \sum_{\sigma\in\SN}
			e^{-\HP(\vecX,\sigma)}}.
$$
These may also be written in terms of the partition functions as
$$
	P_{\Lambda,\vecX}(\pi)
	= \frac{e^{-\HP(\vecX,\pi)}}{Y(\Lambda,\vecX)}
	\qand
	P_{\Lambda,N}(\pi)
	= \frac{\int_{\Lambda^N}d\vecX e^{-\HP(\vecX,\pi)}}{Z(\Lambda,N)N!}.
$$

Note that for the non-interacting $V=0$ case, we have the following heuristic:
\begin{itemize}
\item As $\beta \to 0$, the probability measure becomes supported only on the
\plainidx{identity permutation}.
\item As $\beta \to \infty$, the probability measure approaches the
\plainidx{uniform distribution} on $\SN$.
\end{itemize}

\textbf{Random variables:} Recall from section \ref{sec:history} that Feynman
claimed, and S\"ut\H{o} proved for the ideal gas, that Bose-Einstein
condensation occurs if and only if there  are infinite cycles in an
\plainidx{infinite-volume extension} of the above probability model.
Furthermore, note that cycles depend on the permutation only and not on the
geometric placement of $\vecx_1, \ldots, \vecx_N$.  Thus we will be interested
in random variables $\theta(\pi)$ rather than $\theta(\vecX,\pi)$.  In this
section we confine ourselves to mathematical statements only; the connection
with the physics will be made in section \ref{sec:connect}.  Yet, we have
motivated the choice of the following random variable.

The particular \plainidx{random variable} of interest, $\varrho_{m,n}(\pi)$, is
as follows.  First define
$$
	\ell_i(\pi)
$$
to be the length of the permutation cycle containing the point $\vecx_i$.
For example, for the point configuration $\vecX$ and the permutation $\pi$
in figure \ref{fig:piham}, we have
$$
	\ell_1(\pi) =
	\ell_2(\pi) =
	\ell_3(\pi) =
	\ell_4(\pi) = 4,
	\quad
	\ell_5(\pi) =
	\ell_6(\pi) =
	\ell_7(\pi) = 3,
	\qand
	\ell_8(\pi) = 1.
$$
Also let
$$
	\rho = \frac{N}{V},
$$
i.e. $\rho$ is the \plainidx{particle density}.  For $1 \le m \le n \le N$,
define
$$
	\varrho_{mn}(\pi)
	= \frac{1}{V} \;
	\# \left\{ i=1, \ldots, N: m \le \ell_i(\pi) \le n \right\}
$$
This random variable, taking values between 0 and $\rho$, is the \bfidx{density
of sites} in cycles of specified length.

One may also consider the related random variable
$$
	f_{m,n}
	= \frac{1}{N} \;
	\# \left\{ i=1, \ldots, N: m \le \ell_i(\pi) \le n \right\}
$$
which is $\varrho_{m,n}/\rho$.  This runs from 0 to 1 and is the
\bfidx{fraction of sites} in cycles of specified length.  For figure
\ref{fig:piham}, we have $f_{2,3}(\pi) = 3/8$.

\textbf{Expectations:}
For\index{expectation} a random variable $\theta(\pi)$, we have
$$
	\bbE_{\Lambda,\vecX}(\theta) = \frac{1}{Y,\vecX(\Lambda)}
		\sum_{\pi\in\SN} \theta(\pi) e^{-\HP(\vecX,\pi)}
	\qand
	\bbE_{\Lambda,N}(\theta) = \frac{1}{Z(\Lambda,N) N!}
	\int_{\Lambda^N} d\vecX \sum_{\pi\in\SN} \theta(\pi) e^{-\HP(\vecX,\pi)}.
$$

\textbf{Models:}
The probability measure $P_{\Lambda,\vecX}(\pi)$ and the expectation
$\bbE_{\Lambda,\vecX}(\rho_{m,n})$ is treated in the paper \cite{bib:gru} where
$\vecX$ is a\index{unit lattice} \plainidx{cubic unit lattice};
$P_{\Lambda,N}(\pi)$ and $\bbE_{\Lambda,N}(\rho_{m,n})$ are treated in
\cite{bib:bu} and \cite{bib:qmath}.  The former will be referred to as the
\bfidx{lattice-configuration model}; the latter will be referred to as the
\bfidx{point-process-configuration model}.

%% ----------------------------------------------------------------
\subsection{Existence of infinite cycles}
\label{sec:inf_cyc_existence}

Here\index{infinite cycles} we describe theorem 1 of \cite{bib:qmath}, which is
proved in section 1 of \cite{bib:bu}.  This applies to the $U\equiv 0$ case.

%Here we describe theorem 1 of \cite{bib:qmath}, deferring [XXX?] a sketch
%of the proof to section [XXX].

\textbf{Thermodynamic limit:}  We\index{thermodynamic limit} inquire about the
fraction of sites participating in short and long cycles (as quantified below)
in the \plainidx{infinite-volume limit}.  Namely, we let $V,N \to \infty$ with
fixed ratio $\rho = N/V$, and we ask about the\index{density} cycle-length
distribution as a function of $\rho$.

One does not need to construct an infinite-volume model, although this is done
in section 3, for pure interest:  We will examine limits of expectations of
random variables, where the limit is taken as the number of points $N$ of the
model goes to infinity.  The limits then are in $\R$.

\textbf{Critical density:}
We define the \bfidx{critical density}, $\rho_c{(0)}$, by the following formula.
Here, this is simply an \plainidx{ansatz}.  (See also section
\ref{sec:N2_rho_c} where the critical density is discussed from a thermodynamic
point of view.)
\begin{equation}
\label{eqn:ideal_rho_c}
	\rho_c{(0)} = \int_{\R^d} \frac{d\veck}{e^{4\beta\pi^2\|\veck\|^2}-1}.
\end{equation}

\begin{rem}
One may compute
\begin{align}
\label{eqn:ideal_rho_c_zeta}
	\rho_c{(0)} &= \int_{\R^d} \frac{d\veck}{e^{4\beta\pi^2\|\veck\|^2}-1}
	= \frac{\zeta(3/2)}{(4\pi\beta)^{3/2}}
\end{align}
where $\zeta(z)$ is the \plainidx{Riemann zeta function}.
\end{rem}

% XXX Note that $\exp(-4\beta\pi^2|\veck|^2)$ is the \plainidx{Fourier
% transform} of $\exp(\beta|\vecx|^2/4)$.

% XXX wik bec article:
% N = V int d3k/(2pi)^3 p(k)/(1-p(k))
%   = V int d3k/(2pi)^3 1 / (exp(k^2/2mT) - 1)
% w/ p(k) = exp(-k^2/2mT).
% "maximum number of excited particles p/(1-p)".

\textbf{Statement of the theorem:}
For any $0<a<b<1$ (nominally, $a$ is just above 0 and $b$ is just below 1)
and any $s \ge 0$,
\begin{equation}
\label{eqn:thm_1_rho}
\begin{aligned}
%
	\lim_{V \to \infty} \bbE_{\Lambda,N}(\varrho_{1,N^a}) &=
		\begin{cases}
		\rho,   & \rho \le \rho_c{(0)}\\
		\rho_c{(0)}, & \rho_c{(0)} \le \rho
		\end{cases} \\
%
	\lim_{V \to \infty} \bbE_{\Lambda,N}(\varrho_{N^a,N^b}) &=
		0 \\
%
	\lim_{V \to \infty} \bbE_{\Lambda,N}(\varrho_{N^b,sN}) &=
		\begin{cases}
		0,           & \rho \le \rho_c{(0)} \\
		\rho-\rho_c{(0)}, & \rho_c{(0)} \le \rho \le s+\rho_c{(0)} \\
		s,           & s+\rho_c{(0)} \le \rho.
		\end{cases}
\end{aligned}
\end{equation}
(See figure \ref{fig:limErho}.)

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics{figures/limErho.eps}
\caption{Distribution of cycle lengths as function of density in the
	infinite-volume limit.
	\label{fig:limErho}}
\end{center}
\end{figure}

This statement, in terms of the density of sites in cycles of specified
lengths, may be more intuitively rephrased in terms of the fraction of sites in
cycles of specified length:
\begin{align*}
%
	\lim_{V \to \infty} \bbE_{\Lambda,N}(f_{1,N^a}) &=
		\begin{cases}
		1,   & \rho \le \rho_c{(0)}\\
		\rho_c{(0)}/\rho, & \rho_c{(0)} \le \rho
		\end{cases} \\
%
	\lim_{V \to \infty} \bbE_{\Lambda,N}(f_{N^a,N^b}) &=
		0 \\
%
	\lim_{V \to \infty} \bbE_{\Lambda,N}(f_{N^b,sN}) &=
		\begin{cases}
		0,             & \rho \le \rho_c{(0)} \\
		1-\rho_c{(0)}/\rho, & \rho_c{(0)} \le \rho \le s+\rho_c{(0)} \\
		s/\rho,        & s+\rho_c{(0)} \le \rho.
		\end{cases}
\end{align*}

The theorem has the following interpretation:
\begin{itemize}
\item At density below $\rho_c{(0)}$, all sites are in \plainidx{short cycles}
and no sites are in \plainidx{long cycles}.
\item As density increases past $\rho_c{(0)}$ and $\rho_c{(0)}+s$,
fewer sites are in short cycles and more sites are in long cycles.
In particular, a strictly positive fraction of sites are in long cycles.
\item Asymptotically, all sites are in long cycles.
\end{itemize}

\emph{Late note:}  The recent paper \cite{bib:bu08} produces an expression for
$\rho_c$, as well as an analogue of the previous theorem, for the weakly
interacting ($U \not\equiv 0$) case.

% XXX get the godd*mned units correct!!  units of s?  rescale factor?!?

%%% ================================================================
%\newpage
%\section{Proof of existence of infinite cycles}
%
%XXX \cite{bib:bu}
%
%Description of the steps for proving the existence of infinite cycles. It
%involves certain lemmas that could be stated without proof.
%
%%% ----------------------------------------------------------------
%\subsection{Step \ldots}
%
%%% ----------------------------------------------------------------
%\subsection{Step \ldots}
%
%%% ----------------------------------------------------------------
%\subsection{Step \ldots}
%
%%% ----------------------------------------------------------------
%\subsection{Step \ldots}


%% ----------------------------------------------------------------
\newpage
\section{Bosonic Feynman-Kac formulas}
\label{sec:fkboson}

In this section, the heart of this paper, we recast a familiar
\plainidx{Hamiltonian} involving \textbf{pair-interacting particles} as a new
Hamiltonian involving \textbf{pair-interacting permutations}.  We use
the\index{partition function} \plainidx{canonical partition function} as the
vehicle for this transformation:
$$
	\textrm{Hamiltonian for particles}
	\longto
	\textrm{partition function}
	\longto
	\textrm{Hamiltonian for permuations}.
$$
A Feynman-Kac formula for $N$ interacting \plainidx{bosons} allows the
partition function to be transformed in the middle step.

%% ----------------------------------------------------------------
\subsection{Outline}
\label{sec:fkboson_outline}

As in section \ref{sec:model_defns}, we write
$
	\vecX = (\vecx_1, \ldots, \vecx_N)
$
for $\vecx_1, \ldots, \vecx_N$ in a $d$-dimensional \plainidx{cube} $\Lambda$
of width $L$.  The Hamiltonian for $N$ pair-interacting particles is
\begin{align}
\label{eqn:particle_Hamiltonian}
	H(\vecX) &=
	-\sum_{i=1}^N \nabla_i^2 + \sum_{1\le i,j \le N} U(\vecx_i - \vecx_j).
\end{align}
The \plainidx{$U$} considered in this paper is either identically zero (for the
non-interacting case), or a \bfidx{hard-core potential} with radius $a$, i.e.
$U(\vecx_i-\vecx_j)$ is infinite for $|\vecx_i-\vecx_j| \le a$ and $0$ for
$|\vecr| > a$.  (This is an approximation to the true pair potential between
\plainidx{helium atoms}.  See figure \ref{fig:ceperley_potential}
\cite{bib:ceperley}.)  The hard-core radius $a$ is also known as the
\bfidx{scattering length}.

\begin{figure}[!htb]
\begin{center}
%\psfragscanon
\includegraphics[scale=0.7]{figures/ceperley_potential.eps}
\caption{Pair potential between helium atoms (Ceperley, 1995).
	\label{fig:ceperley_potential}}
\end{center}
\end{figure}

As discussed in appendix \ref{sec:canon_ens_rho}, the partition function for
$N$ distinguishable particles\footnote{For a particle Hamiltonian, the $\beta$
factor is in the expected place.  This is in contrast to the permutation
expression in sections \ref{sec:model_defns}, where the $\beta$ factor is,
surprisingly, \plainidx{reciprocated}.  We will see in section \ref{sec:U0} why
the reciprocated $\beta$ is correct for the permutation Hamiltonian.} is
$\Tr(e^{-\beta H})$.  (Section \ref{sec:canon_ens_rho} explains why $\TrebH$ is
of importance.)  Symmetrizing\index{symmetrizing} the partitition function,
since our particles are bosons, the \plainidx{trace} is
$$
	\Trsym(e^{-\beta H})
	= \Tr_{L^2} \left( P_+ e^{-\beta H} \right)
	= \Tr_{L^2} \left( e^{-\beta H} P_+ \right)
$$
where
$$
	P_+\, f(\vecx_1, \ldots, \vecx_N)
	:= \frac{1}{N!} \sum_{\pi \in \Sn} M_\pi f(\vecx_1, \ldots, \vecx_N)
$$
and
$$
	M_\pi(f\vecx_1, \ldots, \vecx_N)
	:= f(\vecx_{\pi(1)}, \ldots, \vecx_{\pi(N)}).
$$
That is,
$$
	\Trsym(e^{-\beta H})
	= \frac{1}{N!}\;
	\sum_{\pi\in\SN} \Tr_{L^2} \left( e^{-\beta H} M_\pi \right).
$$
(The $e^{-\beta H}$ is bounded\index{bounded operator} and
compact\index{compact operator}, but this fact is not needed.)

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics{figures/wavy5.eps}
% [Short caption for list of figures]
% {Long caption for the body of the document}
\caption[Feynman-Kac representation of a gas of 5 bosons.]{
	Feynman-Kac representation of a gas of 5 bosons. The horizontal plane
	represents the $d$ spatial dimensions, and the vertical axis is the
	imaginary time dimension. The picture shows five particles and two cycles,
	of respective length 4 and 1.
	\label{fig:wavy5}}
\end{center}
\end{figure}


The following steps remain to develop a bosonic \plainidx{Feynman-Kac formula}.
The first three steps closely parallel the steps used to construct the familiar
single-particle Feynman-Kac formula in appendix \ref{sec:fk1}.

\begin{itemize}

\item
Section \ref{sec:fkboson_exp}:
Interpret $e^{-\beta H} M_\pi$ as an expectation over Brownian motions, as in
proposition \ref{prop:ebHG} for the single-particle case.

\item
Section \ref{sec:fkboson_ker}:
Write $e^{-\beta H} M_\pi$ as an integral operator, and find the kernel.

\item
Section \ref{sec:fkboson_tr}:
Compute $\Tr(e^{-\beta H} M_\pi)$ in terms of Brownian bridges.

\item
Section \ref{sec:fkboson_sum}:
Sum over $\pi \in \SN$ to obtain $Z=\Trsym(e^{-\beta H})$.  Importantly, we
will express $Z$ as sum over permutations $\pi$ of $e^{-\HP(\vecX,\pi)}$,
where this new $\HP$ will be viewed a Hamiltonian for a single permuation
$\pi$.  At this point, the permutation Hamiltonian is found inside
$e^{-\HP(\vecX,\pi)}$; we lack an expression for its logarithm.

\item
Sections \ref{sec:U0} and \ref{sec:Unon0} decouple
the non-interacting from the interacting terms in the permutation Hamiltonian,
so that we may write $e^{-\HPzero(\vecX,\pi)-\HPone(\vecX,\pi)}$.

\item
The bosonic Feynman-Kac formula now contains terms for 2-jump interactions,
3-jump interactions, and so on.  In section \ref{sec:cluexp}, we discuss the
cluster expansion which allows us to drop all but 2-jump interactions.  The
cluster expansion furthermore allows us to take the logarithm of
$e^{-\HP(\vecX,\pi)}$, with an explicit expression for $\HP(\vecX,\pi)$.
We recognize the random-cycle model from section \ref{sec:rcm}, with an
explicit 2-jump interaction $V$.

\end{itemize}

%% ----------------------------------------------------------------
\subsection{$e^{-\beta H} M_\pi$ as expectation}
\label{sec:fkboson_exp}

Here\index{expectation} we parallel the development for the single-particle
case in section \ref{sec:fk1_exp}.

\begin{prop}
\label{prop:diffusion_FK_bosonic}
With $H$ as in equation \ref{eqn:particle_Hamiltonian}, we have
\begin{align*}
	e^{-\beta H} f(\vecx_{\pi(1)}, \ldots, \vecx_{\pi(N)})
	&= \bbE_{0}^{\vecx_{\pi(1)}, \ldots, \vecx_{\pi(N)}} \left[
		\exp\left\{
			-\frac{1}{2}
			\sum_{i<j}
				\int_0^{2\beta} U\left(\vecwi_s - \vecwj_s\right)\,ds
		\right\}
		f \left( \vecw_{2\beta}^{(1)}, \ldots, \vecw_{2\beta}^{(N)} \right)
	\right].
\end{align*}
\begin{proof}
Using the \plainidx{Trotter product formula} (equation \ref{eqn:trotter_AB}),
namely
$$
	e^{\beta(A+B)}
	= \lim_{n\to\infty} \left(e^{\beta A/n} e^{\beta B/n}\right)^n
$$
with
$
	A = \sum_{i=1}^N \nabla_i^2
$
and
$
	B = -\sum_{i<j} U(\vecx_i-\vecx_j),
$
we have
\begin{align*}
	&e^{-\beta H} f(\vecx_{\pi(1)},\ldots,\vecx_{\pi(N)}) \\
	&= \lim_{n\to\infty} \left(
		e^{\frac{\beta}{n} \sum_i \nabla_i^2}
		e^{-\frac{\beta}{n} \sum_{i<j} U(\vecx_i-\vecx_j)}
		\right)^n
		f(\vecx_{\pi(1)},\ldots,\vecx_{\pi(N)}) \\
	&= \lim_{n\to\infty}
		e^{\frac{\beta}{n} \sum_i \nabla_i^2}
		e^{-\frac{\beta}{n} \sum_{i<j} U(\vecx_i-\vecx_j)}
	\left(
		e^{\frac{\beta}{n} \sum_i \nabla_i^2}
		e^{-\frac{\beta}{n} \sum_{i<j} U(\vecx_i-\vecx_j)}
		\right)^{n-1}
		f(\vecx_{\pi(1)},\ldots,\vecx_{\pi(N)}).
\end{align*}
Recall that $e^{-\frac{\beta}{n} \sum_{i<j} U(\vecx_i-\vecx_j)}$ is simply a
scalar.  Using the result of section \ref{sec:ebH0} to write
$e^{\frac{\beta}{n} \sum_i \nabla_i^2}$ as an integral operator (since the sum
of $N$ Laplacians, each on $d$ dimensions, is a single $(Nd)$-dimensional
\plainidx{Laplacian}), and writing
$$
	\vecZ^{(k)} = (\vecz_1^{(k)}, \ldots, \vecz_N^{(k)}),
$$
we have
\begin{align*}
	e^{-\beta H} f(\vecx_{\pi(1)},\ldots,\vecx_{\pi(N)}) 
	&= \lim_{n\to\infty}
		\int_{\R^{Nd}}
		g_{2\beta/n}\left( \vecX - \vecZ^{(1)} \right)
		e^{-\frac{\beta}{n} \sum_{i<j} U(\vecz^{(1)}_i-\vecz^{(1)}_j)} \\
	&\qquad \left(
		e^{\frac{\beta}{n} \sum_i \nabla_i^2}
		e^{-\frac{\beta}{n} \sum_{i<j} U(\vecx_i-\vecx_j)}
		\right)^{n-1}
		f(\vecx_{\pi(1)},\ldots,\vecx_{\pi(N)})
		\;d\vecZ^{(1)}.
\end{align*}
Repeating $n-1$ more times yields
\begin{align*}
	& \lim_{n\to\infty}
		\intRdntimes{Nd}{n}
		g_{2\beta/n}\left( \vecX- \vecZ^{(1)} \right)
		\cdots
		g_{2\beta/n}\left( \vecZ^{(n-1)} - \vecZ^{(n)} \right) \\
	& \qquad \left(\prod_{k=1}^n
		e^{-\frac{\beta}{n} \sum_{i<j} U(\vecz^{(k)}_i-\vecz^{(k)}_j)}
	\right)
	f(\vecx_{\pi(1)},\ldots,\vecx_{\pi(N)})
		\;d\vecZ^{(1)} \cdots d\vecZ^{(n)}.
\end{align*}
We have an integrand in the form of remark \ref{rem:ExFomega}, with
$\beta_k = 2k\beta/n$, so we can write
\begin{align*}
%
	& \lim_{n\to\infty}
		\intRdntimes{Nd}{n}
		g_{2\beta/n}\left( \vecX - \vecZ^{(1)} \right)
		\cdots
		g_{2\beta/n}\left( \vecZ^{(n-1)} - \vecZ^{(n)} \right) \\
	& \qquad \exp\left\{
		\frac{2\beta}{n} \left(-\frac{1}{2}\right)
		\sum_{i<j} \sum_{k=1}^n U\left(\vecz^{(k)}_i-\vecz^{(k)}_j\right)
		\right\}
%
	f(\vecx_{\pi(1)},\ldots,\vecx_{\pi(N)})
		\;d\vecZ^{(1)} \cdots d\vecZ^{(n)} \\
%
	&= \bbE_0^{\vecx_{\pi(1)},\ldots,\vecx_{\pi(N)}} \left[
		\exp\left\{
			-\frac{1}{2}
			\sum_{i<j} \int_0^{2\beta}
				U\left(\vecwi_s - \vecwj_s\right) \,ds
		\right\}
	f\left(\vecw_{2\beta}^{(1)},\ldots,\vecw_{2\beta}^{(N)}\right)
	\right].
\end{align*}
\end{proof}
\end{prop}

%% ----------------------------------------------------------------
\subsection{$e^{-\beta H} M_\pi$ as an integral operator}
\label{sec:fkboson_ker}

\begin{prop}
\label{prop:ebHG_bosonic}
If\index{integral operator}
\begin{equation*}
	H = -\sum_i \nabla_i^2 + \sum_{i<j} U(\vecx_i-\vecx_j),
\end{equation*}
then
\begin{equation}
\label{eqn:temp1_bosonic}
\begin{aligned}
	e^{-\beta H} f(\vecx_{\pi(1)}, \ldots, \vecx_{\pi(N)})
	= \int G_{2\beta,U}&(\vecx_{\pi(1)}, \ldots, \vecx_{\pi(N)},
		\vecy_1, \ldots, \vecy_N) \\
		&f(\vecy_1, \ldots, \vecy_N) \; d\vecy_1 \cdots d\vecy_N
\end{aligned}
\end{equation}
where
\begin{equation}
\label{eqn:temp2_bosonic}
\begin{aligned}
	&G_{2\beta,U}(\vecx_{\pi(1)}, \ldots, \vecx_{\pi(N)},
		\vecy_1, \ldots, \vecy_N) = \\
	& \bbE_0^{\vecx_{\pi(1)}, \ldots, \vecx_{\pi(N)}} \left[
		\exp\left\{
			-\frac{1}{2}
			\sum_{i<j}
				\int_0^{2\beta} U\left(\vecwi_s - \vecwj_s\right)\,ds
		\right\}
		\prod_{i=1}^N
		\delta\left(\vecwi_\beta-\vecy^{(i)}\right)
	\right].
\end{aligned}
\end{equation}
\begin{proof}
Insert equation \ref{eqn:temp2_bosonic} into the right-hand side of
\ref{eqn:temp1_bosonic}, interchange expectation and integral, and integrate
out the \plainidx{delta function} as in proposition \ref{prop:ebHG}.
\end{proof}
\end{prop}

%% ----------------------------------------------------------------
\subsection{$\Tr(e^{-\beta H} M_\pi)$ using Brownian bridges}
\label{sec:fkboson_tr}

\begin{prop}
\label{prop:TrebH_bridge_bosonic}
The \plainidx{trace} may be computed using Brownian bridges\index{Brownian
bridge} as follows:
\begin{align*}
	\Tr(e^{-\beta H} M_\pi)
	&=
	\int \,d\vecX
	\int
		\left(
			\prod_{k=1}^N
			d\vecW^{\vecx_k,\vecxpik}_{0,2\beta}\left(\vecwk\right)
		\right)
		\left[
			\exp\left\{
				-\frac{1}{2}
				\sum_{i<j} \int_0^{2\beta}
					U\left(\vecwi_s - \vecwj_s\right)\,ds
			\right\}
		\right].
\end{align*}
\begin{proof}
Using proposition \ref{prop:operator_trace}, we have
$$
	\Tr(e^{-\beta H} M_\pi)
	=
	\int G_{2\beta,U}(
		\vecx_{\pi(1)},\ldots,\vecx_{\pi(N)},
		\vecx_1,\ldots, \vecx_N)
	\,d\vecX.
$$
Equation \ref{eqn:temp2_bosonic} of proposition \ref{prop:ebHG_bosonic} gives
us an expression for $G$.  Then
\begin{align*}
	\Tr(e^{-\beta H} M_\pi)
	=
	\int
	& \bbE_0^{\vecx_{\pi(1)}, \ldots, \vecx_{\pi(N)}} \left[
		\exp\left\{
			-\frac{1}{2}
			\sum_{i<j}
				\int_0^{2\beta} U\left(\vecwi_s - \vecwj_s\right)\,ds
		\right\}
		\prod_{i=1}^N
		\delta\left(\vecwi_{2\beta}-\vecx^{(i)}\right)
	\right]
	\,d\vecX.
\end{align*}
Using proposition \ref{prop:Ex_Exy}, we may convert this \plainidx{expectation}
over \plainidx{Brownian motion} into an expectation over Brownian bridges to
obtain
\begin{align*}
	\Tr&(e^{-\beta H} M_\pi) = \\
	&\prod_{i=1}^{N} g_{2\beta}\left(\vecx_i-\vecxpii\right)
	\int
		\bbE_{0,2\beta}^{\vecx_1,\vecx_{\pi(1)}; \ldots; \vecx_N,\vecx_{\pi(N)}}
		\left[
			\exp\left\{
				-\frac{1}{2}
				\sum_{i<j} \int_0^{2\beta}
					U\left(\vecwi_s - \vecwj_s\right)\,ds
			\right\}
		\right]
	\,d\vecX.
\end{align*}
Using definition \ref{defn:E_exp_U_bb}, we may write this using
$d\vecW$ notation:
\begin{align*}
	\Tr(e^{-\beta H} M_\pi) &=
	\int \,d\vecX
	\int
		\left[
			\prod_{k=1}^N
			d\vecW^{\vecx_k,\vecxpik}_{0,2\beta}\left(\vecwk\right)
		\right]
		\exp\left\{
			-\frac{1}{2}
			\sum_{i<j} \int_0^{2\beta}
				U\left(\vecwi_s - \vecwj_s\right)\,ds
		\right\}.
\end{align*}
\end{proof}
\end{prop}

%% ----------------------------------------------------------------
\subsection{Sum over $\pi \in \SN$}
\label{sec:fkboson_sum}

Applying\index{permutations} proposition \ref{prop:TrebH_bridge_bosonic}, we
continue the plan laid out in section \ref{sec:fkboson}, namely:
\begin{equation*}
\begin{aligned}
	\Trsym(e^{-\beta H})
	&= \frac{1}{N!}\;
		\sum_{\pi\in\SN} \Tr_{L^2} \left( e^{-\beta H} M_\pi \right) \\
	&= \frac{1}{N!}\;
		\int \,d\vecX
		\sum_{\pi\in\SN}
			\left[
				\prod_{k=1}^N
				\int
				d\vecW^{\vecx_k,\vecxpik}_{0,2\beta}
					\left(\vecwk\right)
			\right]
			\exp\left\{
				-\frac{1}{2}
				\sum_{i<j} \int_0^{2\beta}
					U\left(\vecwi_s - \vecwj_s\right)\,ds
			\right\}.
\end{aligned}
\end{equation*}
Notationally, we may split this up as
\begin{equation}
\label{eqn:Tr_L2_ebH}
\begin{aligned}
	\Trsym(e^{-\beta H})
	&= \frac{1}{N!}\;
		\int \,d\vecX
		\sum_{\pi\in\SN}
		e^{-\HP(\vecX, \pi)} \\
	e^{-\HP(\vecX, \pi)}
	&=
		\left[
			\prod_{k=1}^N
			\int
			d\vecW^{\vecx_k,\vecxpik}_{0,2\beta}
				\left(\vecwk\right)
		\right]
		\exp\left\{
			-\frac{1}{2}
			\sum_{i<j} \int_0^{2\beta}
				U\left(\vecwi_s - \vecwj_s\right)\,ds
		\right\}.
\end{aligned}
\end{equation}
Here is the pivotal point of this paper --- the original partition function
appears as a sum over $\pi$ of an $\vecX$-averaged quantity.  That quantity is
non-negative so we may write it as the exponential of something which we call
$\HP$.  The sum over permutations of $e^{-\HP}$ is precisely what we would want
for a \plainidx{partition function} involving\index{energy} energies, not of
\emph{particles}, but of individual \emph{permutations}.  The remaining steps
of this section involve separating out the \plainidx{non-interacting terms}
from the \plainidx{interacting terms} in $e^{-\HP}$, and finding an expression
for the logarithm of $e^{-\HP}$.

%% ----------------------------------------------------------------
\subsection{Extraction of the non-interacting terms}
\label{sec:U0}

In\index{non-interacting terms} this section we obtain
$e^{-\HPzero(\vecX,\pi)}$ as well as its \plainidx{logarithm}.

Let $U \equiv 0$ in equation \ref{eqn:Tr_L2_ebH}.  Then we have
\begin{align*}
	e^{-\HP(\vecX, \pi)}
	&= \left[
		\prod_{k=1}^N
		\int
		d\vecW^{\vecx_k,\vecxpik}_{0,2\beta}
			\left(\vecwk\right)(1)
	\right].
\end{align*}
Recall from proposition \ref{prop:dW_normalization} and definition
\ref{defn:g_t} that
\begin{equation}
\label{eqn:g2bH0}
	\int d\vecW^{\vecx_k,\vecxpik}_{0,2\beta}\left(\vecwk\right)(1)
	= g_{2\beta}(\vecx_k-\vecxpik)
	= \frac{1}{(4\pi\beta)^{d/2}}
		\exp\left\{ -\frac{1}{4\beta} \|\vecx_k-\vecxpik\|^2 \right\}.
\end{equation}
Then
\begin{align*}
	e^{-\HP(\vecX, \pi)}
	&= \left[
		\frac{1}{(4\pi \beta)^{dN/2}}
		\prod_{k=1}^N
		\exp\left\{ -\frac{1}{4\beta} \|\vecx_k-\vecxpik\|^2 \right\}
	\right] \\
	&= \left[
		\frac{1}{(4\pi \beta)^{dN/2}}
		\exp\left\{ -\frac{1}{4\beta}
			\sum_{k=1}^N
			\|\vecx_k-\vecxpik\|^2 \right\}
	\right].
\end{align*}
We write
\begin{align}
\label{eqn:H0pre}
	e^{-\HP(\vecX, \pi)}
	&= \frac{1}{(4\pi \beta)^{dN/2}} e^{-\HPzero(\vecX,\pi)}
\end{align}
where
\begin{align}
\label{eqn:H0}
	\HPzero(\vecX,\pi) &= \frac{1}{4\beta} \sum_{k=1}^N
		\|\vecx_k-\vecxpik\|^2.
\end{align}
(We ignore the \plainidx{prefactor} in equation \ref{eqn:H0pre} since it
cancels out in the computation of expectations of random variables.)

\begin{rem}
Here is the key point where we discover that the $\beta$\index{$\beta$} in the
permutation Hamiltonian is indeed \plainidx{reciprocated} --- in contrast to
our experience with particle Hamiltonians.
\end{rem}

%% ----------------------------------------------------------------
\subsection{Extraction of the interacting terms}
\label{sec:Unon0}

Here\index{interacting terms} we use the result of the previous section to
decompose the exponentiated permutation \plainidx{energy} as
$$
	e^{-\HP(\vecX,\pi)} = e^{-\HPzero(\vecX,\pi)} e^{-\HPone(\vecX,\pi)}.
$$
Unlike section \ref{sec:U0} where we found the \plainidx{logarithm} of
$e^{-\HPzero}$, though, at this point we have only $e^{-\HPone}$.
Finding the logarithm of the latter requires the \plainidx{cluster expansion}
in section \ref{sec:cluexp}.

Equation \ref{eqn:Tr_L2_ebH} is
\begin{align*}
	e^{-\HP(\vecX, \pi)} &=
		\left[
			\prod_{k=1}^N
			\int
			d\vecW^{\vecx_k,\vecxpik}_{0,2\beta}\left(\vecwk\right)
		\right]
		\exp\left\{
			-\frac{1}{2}
			\sum_{i<j} \int_0^{2\beta}
				U\left(\vecwi_s - \vecwj_s\right)\,ds
		\right\}.
\end{align*}
Recall from definition \ref{defn:dWhat} that
\begin{equation*}
	\,d\vecW_{0,2\beta}^{\vecx_i,\vecxpii}\left(\vecwk\right)
	=
	g_{2\beta}(\vecx_i-\vecxpii)
	\,d\vecWhat_{0,2\beta}^{\vecx_i,\vecxpii}\left(\vecwk\right).
\end{equation*}
Combining the two, we have
\begin{align*}
	e^{-\HP(\vecX, \pi)}
%
	&=
		\left[
			\prod_{k=1}^N
			\int
			g_{2\beta}(\vecx_i-\vecxpii)
			\,d\vecWhat_{0,2\beta}^{\vecx_i,\vecxpii}
				\left(\vecwk\right)
		\right]
		\exp\left\{
			-\frac{1}{2}
			\sum_{i<j} \int_0^{2\beta}
				U\left(\vecwi_s - \vecwj_s\right)\,ds
		\right\} \\
%
	&=
		\left[
			\prod_{k=1}^N
			g_{2\beta}(\vecx_i-\vecxpii)
		\right]
		\left[
			\prod_{k=1}^N
			\int \,d\vecWhat_{0,2\beta}^{\vecx_i,\vecxpii}
				\left(\vecwk\right)
		\right]
		\exp\left\{
			-\frac{1}{2}
			\sum_{i<j} \int_0^{2\beta}
				U\left(\vecwi_s - \vecwj_s\right)\,ds
		\right\}.
\end{align*}
Since
$$
	g_{2\beta}(\vecx_k-\vecxpik)
	= \frac{1}{(4\pi\beta)^{d/2}}
		\exp\left\{ -\frac{1}{4\beta} \|\vecx_k-\vecxpik\|^2 \right\}
$$
(equation \ref{eqn:g2bH0}), we obtain
\begin{align*}
	e^{-\HP(\vecX, \pi)}
	&= \frac{1}{(4\pi\beta)^{dN/2}}
		e^{-\HPzero(\vecX,\pi)}
		\left[
			\prod_{k=1}^N
			\int \,d\vecWhat_{0,2\beta}^{\vecx_i,\vecxpii}
				\left(\vecwk\right)
		\right]
		\exp\left\{
			-\frac{1}{2}
			\sum_{i<j} \int_0^{2\beta}
				U\left(\vecwi_s - \vecwj_s\right)\,ds
		\right\}.
\end{align*}
Thus we define
\begin{align}
\label{eqn:expH1}
	e^{-\HPone(\vecX, \pi)} &=
		\left[
			\prod_{k=1}^N
			\int
			d\vecWhat^{\vecx_k,\vecxpik}_{0,2\beta}(\vecwk)
		\right]
		\exp\left\{
			-\frac{1}{2}
			\sum_{i<j} \int_0^{2\beta}
				U\left(\vecwi_s - \vecwj_s\right)\,ds
		\right\}.
\end{align}
This appears very similar to equation \ref{eqn:Tr_L2_ebH}; indeed, the only
notational difference on the right-hand side is the replacement of $d\vecW$ by
$d\vecWhat$.

%% ----------------------------------------------------------------
\subsection{Expansion of multi-jump interactions}

In preparation for the cluster expansion in section \ref{sec:cluexp}, we
reorganize $e^{-\HPone}$.  We develop a heuristic for \plainidx{multi-jump
collision probabilities}, then in section \ref{sec:inclexcl} apply the
\plainidx{inclusion/exclusion principle} to identify $k$-jump interactions for
increasing values of $k$.

These  manipulations may be done formally, with little intuition.
To develop a useful intuition to guide understanding of the manipulations, we
first interpret $e^{-\HPone}$ in terms of \plainidx{collision probabilities}.

\textbf{Result:} The result of this section is that $e^{-\HPone}$ (equation
\ref{eqn:expH1}) may be interpreted as the probability that all $N(N-1)/2$ jump
pairs avoid one another.

The justification for this result is as follows.

\begin{itemize}
\item Recall that $U(\vecr)=\infty$ for $\vecr \le a$; $U(\vecr)=0$ for
$\vecr > a$.

\item If the $i$th and $j$th Brownian-bridge paths ever come within radius $a$
of one another at any \plainidx{Feynman time} between $0$ and $2\beta$, for a
particular realizations of $\vecwi$ and $\vecwj$, then
$$
	\int_0^{2\beta} U\left( \vecwi_s - \vecwj_s \right) ds
$$
is infinite.  If they do not come within radius $a$ of one another, then
that integral is 0.

\item Thus,
$$
	\exp\left\{
		-\frac{1}{2}
		\int_0^{2\beta} U\left( \vecwi_s - \vecwj_s \right) ds
	\right\}
$$
is 0 if the $i$th and $j$th bridges collide, and 1 if they avoid one another.

\item Recall from definition \ref{defn:E_exp_U_bb} that
\begin{align*}
	\int \exp\left\{ \int_0^{2\beta} f(\vecw_s)\,ds \right\}
		\,d\vecW_{0,2\beta}^{\vecx,\vecy}(\vecw)
	&:=
	\bbE_{0,2\beta}^{\vecx,\vecy}
		\left[\exp\left\{\int_0^{2\beta} f(\vecw_s)\,ds\right\} \right].
\end{align*}
The difference of the Brownian bridges $\vecwi$ and $\vecwj$ is twice another
Brownian bridge (proposition \ref{prop:two_bridges}) so definition
\ref{defn:E_exp_U_bb} applies to $\vecwi-\vecwj$.

\item Thus, the expectation
$$
	\left[
	\int d\vecWhat^{\vecx_i,\vecxpii}_{0,2\beta} \left(\vecwi\right)
	\int d\vecWhat^{\vecx_j,\vecxpij}_{0,2\beta} \left(\vecwj\right)
	\right]
	\exp\left\{
		-\frac{1}{2}
		\int_0^{2\beta}
			U\left(\vecwi_s - \vecwj_s\right)\,ds
	\right\}
$$
is the \bfidx{avoidance probability} for the bridges $\vecwi$ and
$\vecwj$.

\item There are $N$ permutation jumps and $N(N-1)/2$ distinct jump pairs, so
\begin{align*}
	e^{-\HPone(\vecX, \pi)} &=
		\left[
			\prod_{k=1}^N
			\int
			d\vecWhat^{\vecx_k,\vecxpik}_{0,2\beta}
				\left(\vecwk\right)
		\right]
		\prod_{i<j}
		\exp\left\{
			-\frac{1}{2}
			\int_0^{2\beta} U\left(\vecwi_s - \vecwj_s\right)\,ds
		\right\}
\end{align*}
(equation \ref{eqn:expH1}) is the probability that all $N(N-1)/2$ jump pairs
avoid one another.

\end{itemize}

% XXX need a figure here.

%% ----------------------------------------------------------------
\subsection{Inclusion/exclusion and collision probabilities}
\label{sec:inclexcl}

The goal of this section is to reinterpret the \plainidx{avoidance probability}
of the previous section as an \plainidx{alternating sum} of \plainidx{collision
probabilities}.

% XXX summary here.

\begin{defn}
Define
$$
	\Upsilon\left(\vecwi -\vecwj \right)
	=
	1 -
	\exp\left\{
		-\frac{1}{2}
		\int_0^{2\beta}
			U\left(\vecwi_s - \vecwj_s\right)\,ds
	\right\}.
$$
\end{defn}

We will frequently abbreviate this\index{$\Upsilon$} as $\Upsilon_{ij}$.  In
light of the discussion in the previous section, the expected value of this is
the \bfidx{collision probability} for $\vecwi$ and $\vecwj$.

Then equation \ref{eqn:expH1} becomes
\begin{align}
\label{eqn:Upsilon_pre_FOIL}
	e^{-\HPone(\vecX, \pi)} &=
		\left[
			\prod_{k=1}^N
			\int
			d\vecWhat^{\vecx_k,\vecxpik}_{0,2\beta}
				\left(\vecwk\right)
		\right]
		\prod_{i<j} \left( 1 - \Upsilon\left(\vecwi-\vecwj\right) \right).
\end{align}

Using \plainidx{de Morgan's law}, we have
$$
	P(\textrm{all avoid}) = 1 - P(\textrm{any collide}).
$$
Let $C_\ell$ be the event that the $\ell$th pair $(\vecwi,\vecwj)$
collides.  There are $N(N-1)/2$ such events.
The latter term is a union of events, so we may use the
\plainidx{inclusion/exclusion principle} to write
\begin{align*}
	1 - P(\cup C_\ell)
	&= 1 - \sum_i P(C_i)
	+ \sum_{ij} P(C_i \cap C_j)
	- \sum_{ijk} P(C_i \cap C_j \cap C_k)
	+ \ldots \\
	&+ (-1)^{m} \sum_{m-\textrm{tuples of jump pairs}}
		P(m \textrm{ bridges collide})
	+ \ldots \\
	&+ (-1)^{N(N-1)/2} P(\textrm{all bridges collide}).
\end{align*}
Collecting summation symbols, i.e. ranking the terms of $e^{-\HPone}$ by
increasing $m$ for probabilities of $m$-tuples of jump pairs colliding,
we have
\begin{align}
\label{eqn:Upsilon_post_FOIL}
	e^{-\HPone(\vecX,\pi)}
	&=
	\left[
		\prod_{k=1}^N
		\int d\vecWhat^{\vecx_k,\vecxpik}_{0,2\beta} \left(\vecwk\right)
	\right]
	\sum_{m=0}^{N(N-1)/2}
	(-1)^m
	\sum_{(i_1,j_1),\ldots,(i_m,j_m)}
	\prod_{\ell=1}^{m} \Upsilon\left(\vecwil-\vecwjl\right)
\end{align}
The first sum is over sizes of subsets of the $N(N-1)/2$ jump pairs;
the second sum is over all possible ways of selecting $m$ pairs.

For example, with $N=3$, there are $N(N-1)/2=3$ jump pairs:  $1,2$;
$1,3$; $2,3$.  Then
\begin{align*}
	e^{-\HPone(\vecX,\pi)}
	&=
	1  \\
	&-P(1,2 \textrm{ collide})
	-P(1,3 \textrm{ collide})
	-P(2,3 \textrm{ collide}) \\
	&+P(1,2 \textrm{ and } 1,3 \textrm{ collide})
	+P(1,2 \textrm{ and } 2,3 \textrm{ collide})
	+P(1,3 \textrm{ and } 2,3 \textrm{ collide}) \\
	&-P(1,2 \textrm{ and } 1,3 \textrm{ and } 2,3 \textrm{ collide}).
\end{align*}
Here there is one $m$ value per line; the different $P$'s on that line
are indexed by $\ell$.  In terms of $\Upsilon$'s, we have
\begin{align*}
	e^{-\HPone(\vecX,\pi)}
	&=
	\left[
	\int d\vecWhat^{\vecx_1,\vecx_{\pi(1)}}_{0,2\beta} \left(\vecw^{(1)}\right)
	\int d\vecWhat^{\vecx_2,\vecx_{\pi(2)}}_{0,2\beta} \left(\vecw^{(2)}\right)
	\int d\vecWhat^{\vecx_3,\vecx_{\pi(3)}}_{0,2\beta} \left(\vecw^{(3)}\right)
	\right] \\
%
	&
	\left(
	\underbrace{1}_{m=0}
	- \underbrace{(\Upsilon_{12} + \Upsilon_{13} + \Upsilon_{23})}_{m=1}
	+ \underbrace{(\Upsilon_{12} \Upsilon_{13}
		+ \Upsilon_{12} \Upsilon_{23}
		+ \Upsilon_{13} \Upsilon_{23})}_{m=2}
	- \underbrace{\Upsilon_{12} \Upsilon_{13} \Upsilon_{23}}_{m=3}
	\right).
\end{align*}
This collection of $m$-size subsets of the set of jump pairs is the difference
between equations \ref{eqn:Upsilon_pre_FOIL} and \ref{eqn:Upsilon_post_FOIL}.

This expression has been obtained using the inclusion/exclusion principle.
Notice, however, that one obtains the same expression by distributing the
product
$
	(1-\Upsilon_{12}) (1-\Upsilon_{13}) (1-\Upsilon_{23})
$
(for $N=3$) or, in general,
$
	\prod_{i<j} (1 - \Upsilon_{ij}).
$

%% ----------------------------------------------------------------
\subsection{Heuristic for cluster expansion}
\label{sec:precluexp}

Moving the integrals through the sums in equation \ref{eqn:Upsilon_post_FOIL},
we obtain
\begin{align}
\label{eqn:precluexp_almost}
	e^{-\HPone(\vecX,\pi)}
	&=
	\sum_{m=0}^{N(N-1)/2}
	(-1)^m
	\sum_{(i_1,j_1),\ldots,(i_m,j_m)}
	\left[
		\prod_{k=1}^N
		\int d\vecWhat^{\vecx_k,\vecxpik}_{0,2\beta} \left(\vecwk\right)
	\right]
	\prod_{\ell=1}^{m}
	\Upsilon\left(\vecwil-\vecwjl\right).
\end{align}
The heuristic (formalized by the \plainidx{cluster expansion} as described in
section \ref{sec:cluexp}) is that one may form an approximation by also moving
the integrals through the leftmost product.

(For example, let $N=3$.  For shorthand, let
\begin{align*}
	\int_k =
	d\vecWhat^{\vecx_k,\vecx_{\pi(k)}}_{0,2\beta} \left(\vecw^{(k)}\right).
\end{align*}
Then
\begin{align*}
	e^{-\HPone(\vecX,\pi)}
%
	&=
	\left[
	\int_1
	\int_2
	\int_3
	\right]
	\left(
		1
		- \Upsilon_{12} - \Upsilon_{13} - \Upsilon_{23}
		+ \Upsilon_{12} \Upsilon_{13}
		+ \Upsilon_{12} \Upsilon_{23}
		+ \Upsilon_{13} \Upsilon_{23}
		- \Upsilon_{12} \Upsilon_{13} \Upsilon_{23}
	\right) \\
%
	&=
		1
		- \int_1 \int_2 \Upsilon_{12}
		- \int_1 \int_3 \Upsilon_{13}
		- \int_2 \int_3 \Upsilon_{23} \\
	&	+ \int_1 \int_2 \int_3
		\Upsilon_{12} \Upsilon_{13}
		+ \int_1 \int_2 \int_3
		\Upsilon_{12} \Upsilon_{23}
		+ \int_1 \int_2 \int_3
		\Upsilon_{13} \Upsilon_{23}
		- \int_1 \int_2 \int_3
		\Upsilon_{12} \Upsilon_{13} \Upsilon_{23}.
\end{align*}
The integrals move trivially through the product for $m=0,1$;
for $m\ge 2$, the approximations are of the form
\begin{align}
\label{eqn:precluexp_approx}
	 \left[\int_1 \int_2 \int_3 \Upsilon_{12} \Upsilon_{13}\right]
	 &\approx
	 \left[\int_1 \int_2 \Upsilon_{12}\right]
	 \left[\int_1 \int_3 \Upsilon_{13}\right]
\end{align}
and so on.  Dropping the $m\ge 2$ terms depends on the smallness of this
modification.  This step is formalized by the cluster expansion in section
\ref{sec:cluexp}.  Intuitively, recall that $\Upsilon_{12}$ is the probability
that bridges $1$ and $2$ collide, while $\Upsilon_{12} \Upsilon_{13}$ is the
probability of a collision between bridges $1$ and $2$ \emph{and} bridges $1$
and $3$.  The approximation condition in equation \ref{eqn:precluexp_approx} is
that, for weak interactions, these events are weakly correlated and hence
nearly independent.
This is always true when $ij$'s do not overlap.  E.g. for $N=4$,
$
	\left[ \int_1\int_2\int_3\int_4 \Upsilon_{12}\Upsilon_{34}\right]
	=
	\left[ \int_1\int_2 \Upsilon_{12}\right]
	\left[ \int_3\int_4 \Upsilon_{34}\right].
$
The approximation applies to overlapping $ij$'s.)

Once the integrals have been moved through the leftmost product of equation
\ref{eqn:precluexp_almost}, the expectation of $\Upsilon(\vecw^{(i_\ell)} -
\vecw^{(j_\ell)})$ depends only on the Brownian bridges for $\vecw^{(i_\ell)}$
and $\vecw^{(j_\ell)}$.

We define
\begin{align}
\label{eqn:precluexp_V}
	V_{ij} &=
	V(\vecx_i,\vecxpii,\vecx_j,\vecxpij) =
	\left[
	\int d\vecWhat^{\vecx_i,\vecxpii}_{0,2\beta} \left(\vecwi\right)
	\int d\vecWhat^{\vecx_j,\vecxpij}_{0,2\beta} \left(\vecwj\right)
	\right]
	\Upsilon\left(\vecwi-\vecwj\right).
\end{align}
Since we assume \plainidx{small interactions} $V_{ij}$,
\begin{equation}
\label{eqn:precluexp_exp_sum}
\begin{aligned}
	e^{-\HPone(\vecX,\pi)}
	&\approx
	\prod_{i<j} \left(1 - V_{ij}\right) \\
%
	&\approx
	\prod_{i<j} \left(
	1 - V_{ij} + \frac{V_{ij}^2}{2} - \frac{V_{ij}^3}{6} + \ldots\right)
%
	= \prod_{i<j} e^{-V_{ij}}
%
	= \exp\left\{-\sum_{i<j} V_{ij}\right\}.
\end{aligned}
\end{equation}
Here, the first $\approx$ comes from moving the stochastic integrals through
the product as discussed at the start of this section, and the second $\approx$
arises from the insertion of negligible \plainidx{higher-order terms} in the
$V_{ij}$'s.

For example, with $N=3$, we have
\begin{align*}
	(1 - V_{12})
	(1 - V_{13})
	(1 - V_{23})
	&=
	1
	-(V_{12} +V_{13} +V_{23})
	+ (V_{12} V_{13} +V_{12} V_{23} +V_{13} V_{23})
	-V_{12} V_{13} V_{23}
\end{align*}
which we compare to
\begin{align*}
&
\left(1 -V_{12} +\frac{V_{12}^2}{2} - \ldots \right)
\left(1 -V_{13} +\frac{V_{13}^2}{2} - \ldots \right)
\left(1 -V_{23} +\frac{V_{23}^2}{2} - \ldots \right) \\
%
&=1
%
-\left(V_{12} +V_{13} +V_{23} \right)
%
+\left(V_{12} V_{13} + V_{12} V_{23} + V_{13} V_{23} \right)
%
-V_{12} V_{13} V_{23} \\
%
&+\left( \frac{V_{12}^2} {2}
+\frac{V_{13}^2} {2}
+\frac{V_{23}^2} {2} \right)
%
-\left(
 \frac{V_{12}^2 V_{13}} {2}
+\frac{V_{12} V_{13}^2} {2}
+\frac{V_{12}^2 V_{23}} {2}
+\frac{V_{12} V_{23}^2} {2}
+\frac{V_{13}^2 V_{23}} {2}
+\frac{V_{13} V_{23}^2} {2}
\right) \\
%
&+\left(
 \frac{V_{12}^2  V_{13}^2} {4}
+\frac{V_{12}^2  V_{23}^2} {4}
+\frac{V_{13}^2  V_{23}^2} {4}
\right) \\
%
&+\left(
 \frac{V_{12}^2 V_{13}    V_{23}}   {2}
+\frac{V_{12}   V_{13}^2  V_{23}}   {2}
+\frac{V_{12}   V_{13}    V_{23}^2} {2}
\right)
%
-\left(
 \frac{V_{12}   V_{13}^2  V_{23}^2}  {4}
+\frac{V_{12}^2 V_{13}^2  V_{23}}    {4}
+\frac{V_{12}^2 V_{13}    V_{23}^2 } {4}
\right)
%
+\frac{V_{12}^2  V_{13}^2  V_{23}^2} {8}
+ \ldots.
\end{align*}

We have now achieved one of the goals outlined in section
\ref{sec:fkboson_outline}, namely, to find the \plainidx{logarithm} of
$e^{\HPone}$.  Equation \ref{eqn:precluexp_exp_sum} becomes
\begin{align*}
	\HPone(\vecX,\pi)
	&= \sum_{1 \le i < j \le N} V(\vecx_i,\vecxpii,\vecx_j,\vecxpij)
\end{align*}
where the \plainidx{jump-pair interaction} is as in equation
\ref{eqn:precluexp_V}.

%% ----------------------------------------------------------------
\subsection{Cluster expansion and jump-pair interactions}
\label{sec:cluexp}

This section formalizes the heuristic in section \ref{sec:precluexp}.  This
topic is addressed in section 4 of \cite{bib:qmath}.  Given time constraints
leading up to my comprehensive examination, I will first finish other sections
of this paper, then complete this section if time permits.

%% ----------------------------------------------------------------
\subsection{Simplified jump-pair interactions}
\label{sec:2w_to_1w}

When one computes the jump-pair interaction, it is possible to replace the
double Brownian bridge by a single \plainidx{Brownian bridge}.  (See also
figure \ref{fig:2w_to_1w}.)

\begin{figure}[!htb]
\begin{center}
\psfragscanon
\includegraphics{figures/2w_to_1w.eps}
\caption{Replacement of double Brownian bridge by single Brownian bridge
	for jump-pair interactions.
	\label{fig:2w_to_1w}}
\end{center}
\end{figure}

\begin{prop}
\label{prop:2w_to_1w}
The jump-pair interaction $V_{ij} = V(\vecx_i,\vecxpii,\vecx_j, \vecxpij)$
(equation \ref{eqn:precluexp_V}) satisfies
\begin{align*}
	&
	\int d\vecWhat^{\vecx_i,\vecxpii}_{0,2\beta}(\vecwi)
	\int d\vecWhat^{\vecx_j,\vecxpij}_{0,2\beta}(\vecwj)
	\left[
	1 - 
		\exp\left\{
			-\frac{1}{2}
			\int_0^{2\beta}
				U\left(\vecwi_s - \vecwj_s\right)\,ds
		\right\}
	\right] \\
%
	=& \;
	\int d\vecWhat^{\vecx_i-\vecx_j,\vecxpii,\vecxpij}_{0,4\beta} (\vecwij)
	\left[
	1 - 
		\exp\left\{
			-\frac{1}{4}
			\int_0^{4\beta}
				U\left(\vecw^{(ij)}_s\right)\,ds
		\right\}
	\right].
\end{align*}
\begin{proof}
From the definition of $\int d\vecW$ and $\int d\vecWhat$ (definitions
\ref{defn:E_exp_U_bb} and \ref{defn:dWhat}, respectively), we have
\begin{align}
\label{eqn:2w_to_1w_normalization}
	\int d\vecWhat^{\vecx_i,\vecxpii}_{0,2\beta}(\vecwi)
	\int d\vecWhat^{\vecx_j,\vecxpij}_{0,2\beta}(\vecwj)
	&=
	\frac{
		\int d\vecW^{\vecx_i,\vecxpii}_{0,2\beta}(\vecwi)
	}{g_{2\beta}(\vecx_i-\vecxpii)}
	\frac{
		\int d\vecW^{\vecx_j,\vecxpij}_{0,2\beta}(\vecwj)
	}{g_{2\beta}(\vecx_j-\vecxpij)}
\end{align}
and
\begin{equation}
\label{eqn:2w_to_1w_dWdW}
\begin{aligned}
	&
	\int d\vecW^{\vecx_i,\vecxpii}_{0,2\beta}(\vecwi)
	\int d\vecW^{\vecx_j,\vecxpij}_{0,2\beta}(\vecwj)
	\left[
	1 - 
		\exp\left\{
			-\frac{1}{2}
			\int_0^{2\beta}
				U\left(\vecwi_s - \vecwj_s\right)\,ds
		\right\}
	\right] \\
	&=
	\lim_{n\to\infty}
		\intRdntimes{d}{n-1}
		\intRdntimes{d}{n-1}
	\; d\vecz^{(i)}_1 \cdots d\vecz^{(i)}_{n-1}
	\; d\vecz^{(j)}_1 \cdots d\vecz^{(j)}_{n-1} \\
	&
		g_{2\beta/n}\left(\vecx_i-\vecz^{(i)}_1\right)
		g_{2\beta/n}\left(\vecz^{(i)}_1-\vecz^{(i)}_2\right)
		\cdots
		g_{2\beta/n}\left(\vecz^{(i)}_{n-2}-\vecz^{(i)}_{n-1}\right)
		g_{2\beta/n}\left(\vecz^{(i)}_{n-1}-\vecxpii\right) \\
	&
		g_{2\beta/n}\left(\vecx_j-\vecz^{(j)}_1\right)
		g_{2\beta/n}\left(\vecz^{(j)}_1-\vecz^{(i)}_2\right)
		\cdots
		g_{2\beta/n}\left(\vecz^{(j)}_{n-2}-\vecz^{(j)}_{n-1}\right)
		g_{2\beta/n}\left(\vecz^{(j)}_{n-1}-\vecxpij\right)
		\\
	&
		\left[
		1 -
		\exp\left\{
			\frac{2\beta}{n}
			\left(-\frac{1}{2}\right)
			\left(
				U(\vecx_i-\vecx_j) + \sum_{k=1}^{n-1}
				U\left(\vecz^{(i)}_k - \vecz^{(j)}_k\right)
			\right)
		\right\}
		\right].
\end{aligned}
\end{equation}
For brevity, let
\begin{align*}
	\vecz^{(i)}_0 &= \vecx_i, &\vecz^{(i)}_n &= \vecxpii \\
	\vecz^{(j)}_0 &= \vecx_j, &\vecz^{(j)}_n &= \vecxpij.
\end{align*}
Then the right-hand side of equation \ref{eqn:2w_to_1w_dWdW} becomes
\begin{equation*}
\begin{aligned}
	&
	\lim_{n\to\infty}
		\intRdntimes{d}{n-1}
		\intRdntimes{d}{n-1}
	\; d\vecz^{(i)}_1 \cdots d\vecz^{(i)}_{n-1}
	\; d\vecz^{(j)}_1 \cdots d\vecz^{(j)}_{n-1} \\
	&
		\prod_{k=0}^{n-1}
		g_{2\beta/n}\left(\vecz^{(i)}_{k}-\vecz^{(i)}_{k+1}\right)
		\prod_{k=0}^{n-1}
		g_{2\beta/n}\left(\vecz^{(j)}_{k}-\vecz^{(j)}_{k+1}\right)
	\left[
		1 -
		\exp\left\{
			\frac{2\beta}{n}
			\left(-\frac{1}{2}\right)
			\sum_{k=0}^{n-1} U\left(\vecz^{(i)}_k - \vecz^{(j)}_k\right)
		\right\}
	\right].
\end{aligned}
\end{equation*}

We now make a change of variable.  For $k=0, \ldots, n$, let
$$
	\vecy^{(ij)}_k = \vecz^{(i)}_k - \vecz^{(j)}_k.
$$
The linear map sending
$\left(\vecz^{(i)}_k, \vecz^{(j)}_k\right)$ to
$\left(\vecz^{(i)}_k, \vecy^{(ij)}_k\right)$
is
$$
	\prmatrix{\vecz^{(i)}_k \\ \vecy^{(ij)}_k}
	=
	\prmatrix{1 & 0 \\ 1 & -1}
	\prmatrix{\vecz^{(i)}_k \\ \vecz^{(j)}_k}
$$
which has Jacobian determinant of absolute value 1, so we may replace
$$
	d\vecz^{(i)}_1 \cdots d\vecz^{(i)}_{n-1}
	d\vecz^{(j)}_1 \cdots d\vecz^{(j)}_{n-1}
$$
with
$$
	d\vecz^{(i)}_1 \cdots d\vecz^{(i)}_{n-1}
	d\vecy^{(ij)}_1 \cdots d\vecy^{(ij)}_{n-1}.
$$
Lemma \ref{lem:g2tg2t_g4t_gt_ab} tells us
\begin{align*}
	&
	g_{2\beta/n}\left(\vecz^{(i)}_{k} - \vecz^{(i)}_{k+1}\right)
	g_{2\beta/n}\left(\vecz^{(j)}_{k} - \vecz^{(j)}_{k+1}\right) \\
%
	=& \;
	g_{4\beta/n}\left(
		(\vecz^{(i)}_k-\vecz^{(j)}_k) - (\vecz^{(i)}_{k+1}-\vecz^{(j)}_{k+1})
	\right)
	\;
	g_{ \beta/n}\left(\vecz^{(i)}_k-\vecz^{(i)}_{k+1}
		- \frac{1}{2}\left(\vecz^{(i)}_k-\vecz^{(j)}_k\right)
		+ \frac{1}{2}\left(\vecz^{(i)}_{k+1}-\vecz^{(j)}_{k+1}\right) \right) \\
%
	=& \;
	g_{4\beta/n}\left( \vecy^{(ij)}_k - \vecy^{(ij)}_{k+1} \right)
	\;
	g_{ \beta/n}\left(\vecz^{(i)}_k-\vecz^{(i)}_{k+1}
		- \frac{1}{2} \vecy^{(ij)}_k
		+ \frac{1}{2} \vecy^{(ij)}_{k+1} \right) \\
%
	=& \;
	g_{4\beta/n}\left( \vecy^{(ij)}_k - \vecy^{(ij)}_{k+1} \right)
	\;
	g_{ \beta/n}\left(\left(\vecz^{(i)}_k-\vecz^{(i)}_{k+1}\right)
		- \frac{1}{2} \left(\vecy^{(ij)}_k - \vecy^{(ij)}_{k+1} \right)
		\right).
\end{align*}
Then we have
\begin{equation*}
\begin{aligned}
	&
	\lim_{n\to\infty}
		\intRdntimes{d}{n-1}
		\intRdntimes{d}{n-1}
	\; d\vecz^{(i)}_1 \cdots d\vecz^{(i)}_{n-1}
	\; d\vecy^{(ij)}_1 \cdots d\vecy^{(ij)}_{n-1} \\
	&
	\prod_{k=0}^{n-1}
		g_{4\beta/n}\left( \vecy^{(ij)}_k - \vecy^{(ij)}_{k+1} \right)
	\prod_{k=0}^{n-1}
		g_{ \beta/n}\left(\left(\vecz^{(i)}_k-\vecz^{(i)}_{k+1}\right)
			- \frac{1}{2} \left(\vecy^{(ij)}_k - \vecy^{(ij)}_{k+1} \right)
			\right)
	\left[
		1 -
		\exp\left\{
			-\frac{\beta}{n}
			\sum_{k=0}^{n-1} U\left(\vecy^{(ij)}_k\right)
		\right\}
	\right].
\end{aligned}
\end{equation*}
Re-arranging terms (permissible by \plainidx{Tonelli's theorem} due to the
non-negativity of the integrand) yields
\begin{equation}
\label{eqn:2w_to_1w_split_up}
\begin{aligned}
	&
	\lim_{n\to\infty}
		\intRdntimes{d}{n-1}
	\; d\vecy^{(ij)}_1 \cdots d\vecy^{(ij)}_{n-1}
	\prod_{k=0}^{n-1}
		g_{4\beta/n}\left( \vecy^{(ij)}_k - \vecy^{(ij)}_{k+1} \right)
%
	\left[
		1 -
		\exp\left\{
			\frac{4\beta}{n}
			\left(-\frac{1}{4}\right)
			\sum_{k=0}^{n-1} U\left(\vecy^{(ij)}_k\right)
		\right\}
	\right] \\
%
	&
		\intRdntimes{d}{n-1}
	\; d\vecz^{(i)}_1 \cdots d\vecz^{(i)}_{n-1}
%
	\prod_{k=0}^{n-1}
		g_{ \beta/n}\left(\left(\vecz^{(i)}_k-\vecz^{(i)}_{k+1}\right)
			- \frac{1}{2} \left(\vecy^{(ij)}_k - \vecy^{(ij)}_{k+1} \right)
			\right).
\end{aligned}
\end{equation}

The second line of equation \ref{eqn:2w_to_1w_split_up} collapses by the
iterated Chapman-Kolmogorov\index{Chapman-Kolmogorov equation} proposition
(proposition \ref{prop:g_iter}):
\begin{align*}
	&
	\intRdntimes{d}{n-1}
	\; d\vecz^{(i)}_1 \cdots d\vecz^{(i)}_{n-1}
%
	\prod_{k=0}^{n-1}
		g_{ \beta/n}\left(\left(\vecz^{(i)}_k-\vecz^{(i)}_{k+1}\right)
			- \frac{1}{2} \left(\vecy^{(ij)}_k - \vecy^{(ij)}_{k+1} \right)
			\right) \\
%
	&
	\intRdntimes{d}{n-1}
	\; d\vecz^{(i)}_1 \cdots d\vecz^{(i)}_{n-1}
%
	\prod_{k=0}^{n-1}
		g_{ \beta/n}\left(
		\left( \vecz^{(i)}_k - \frac{\vecy^{(ij)}_k}{2} \right)
		-
		\left( \vecz^{(i)}_{k+1} - \frac{\vecy^{(ij)}_{k+1}}{2} \right)
		\right) \\
%
	&= \;
	g_{\beta}\left(
		\left( \vecz^{(i)}_0 - \frac{\vecy^{(ij)}_0}{2} \right)
		-
		\left( \vecz^{(i)}_n - \frac{\vecy^{(ij)}_n}{2} \right)
	\right)
%
	= \;
	g_{\beta}\left(
		\left( \vecx_i - \frac{\vecx_i-\vecx_j}{2} \right)
		-
		\left( \vecxpii - \frac{\vecxpii-\vecxpij}{2} \right)
	\right) \\
%
	&= \;
	g_{\beta}\left(
		\frac{\vecx_i+\vecx_j}{2} - \frac{\vecxpii+\vecxpij}{2}
	\right).
\end{align*}

With this scale factor pulled out of equation \ref{eqn:2w_to_1w_split_up}, and
recognizing the first line of equation \ref{eqn:2w_to_1w_split_up} from
definition \ref{defn:E_exp_U_bb}, we now have
\begin{equation}
\label{eqn:2w_to_1w_split_pre_norm}
\begin{aligned}
	g_{\beta}\left(
		\frac{\vecx_i+\vecx_j}{2} - \frac{\vecxpii+\vecxpij}{2}
	\right)
	&
	\int d\vecW_{0,4\beta}^{\vecx_i-\vecx_j,\vecxpii-\vecxpij}(\vecw^{(ij)})
	\left[
		1 -
		\exp\left\{
			-\frac{1}{4}
			\int_{0}^{4\beta}
			U\left(\vecw^{(ij)}_s \,ds \right)
		\right\}
	\right].
\end{aligned}
\end{equation}

It remains to restore the \plainidx{normalizations} from equation
\ref{eqn:2w_to_1w_normalization}.  From corollary \ref{cor:g2tg2t_g4t_gt_xy} we
have
\begin{align*}
	g_{2\beta}(\vecx_i-\vecxpii) \; g_{2\beta}(\vecx_j-\vecxpij)
	&=
	g_{4\beta}\left((\vecx_i-\vecx_j) - (\vecxpii-\vecxpij)\right)
	\;
	g_{ \beta}\left(\frac{(\vecx_i+\vecx_j) - (\vecxpii+\vecxpij)}{2}\right).
\end{align*}
The $g_{\beta}$ factor is cancelled by the $g_\beta$ factor in equation
\ref{eqn:2w_to_1w_split_pre_norm}; the $g_{4\beta}$ factor is precisely the
normalization factor to convert equation \ref{eqn:2w_to_1w_split_pre_norm}'s
$d\vecW$ to $d\vecWhat$.  (See definition \ref{defn:dWhat}.) Then equation
\ref{eqn:2w_to_1w_split_pre_norm} becomes
\begin{equation*}
\begin{aligned}
	\int d\vecWhat_{0,4\beta}^{\vecx_i-\vecx_j,\vecxpii-\vecxpij}(\vecw^{(ij)})
	\left[
		1 -
		\exp\left\{
			-\frac{1}{4}
			\int_{0}^{4\beta}
			U\left(\vecw^{(ij)}_s \,ds \right)
		\right\}
	\right]
\end{aligned}
\end{equation*}
which is what was to be shown.
\end{proof}
\end{prop}

%% ----------------------------------------------------------------
\subsection{Relation to spatial permutations}
\label{sec:connect}

Sections \ref{sec:U0} through \ref{sec:2w_to_1w} connect the permutation
Hamiltonian, obtained\index{spatial permutations} via the bosonic Feynman-Kac
formula, with the \plainidx{random-cycle model} as defined in section
\ref{sec:rcm}.  Specifically, we have, in correspondence with equation
\ref{eqn:rcm_H},
\begin{align}
\label{eqn:connect_H}
	H(\vecX,\pi) &= \sum_{i=1}^N \frac{1}{4\beta}\|\vecx_i-\vecxpii\|^2
		+ \sum_{1\le i<j \le N}
			V(\vecx_i, \vecxpii, \vecx_j, \vecxpij)
\end{align}
where
\begin{align}
\label{eqn:connect_V}
	V(\vecx_i,\vecxpii,\vecx_j, \vecxpij)
	=& \;
	\int d\vecWhat^{\vecx_i-\vecx_j,\vecxpii,\vecxpij}_{0,4\beta} (\vecwij)
	\left[
	1 - \exp\left\{
		-\frac{1}{4} \int_0^{4\beta} U\left(\vecw^{(ij)}_s\right)\,ds
	\right\}
	\right].
\end{align}

%% ================================================================
\newpage
\section{Model with full jump-pair interactions}
\label{sec:tji}

Equation \ref{eqn:connect_V} for the \plainidx{jump-pair interactions} was
\begin{align*}
	V(\vecx_i,\vecxpii,\vecx_j, \vecxpij)
	=& \;
	\int d\vecWhat^{\vecx_i-\vecx_j,\vecxpii,\vecxpij}_{0,4\beta} (\vecwij)
	\left[
	1 - \exp\left\{
		-\frac{1}{4} \int_0^{4\beta} U\left(\vecw^{(ij)}_s\right)\,ds
	\right\}
	\right].
\end{align*}
This may be interpreted as the probability that a \plainidx{Brownian bridge}
(figure \ref{fig:2w_to_1w} on page \pageref{fig:2w_to_1w}) running from
$\vecx_i-\vecx_j$ to $\vecxpii-\vecxpij$ in time $4\beta$ intersects a
\plainidx{ball of radius $a$} centered at the origin.  (Proposition
\ref{prop:2w_to_1w} showed that this is exactly equal to the
\plainidx{collision probability} for two bridges, one running from $\vecx_i$ to
$\vecxpii$ and another running from $\vecx_j$ to $\vecxpij$, in time $4\beta$.)

Ueltschi's 2007 paper \cite{bib:qmath} has little more to say about this full
jump-pair interaction.  There are (at least) three things which can be done
with it:

\begin{itemize}

\item Computing it directly using simulation methods is an interesting
statistical problem which will be discussed further in my dissertation.  (To
date, I have found that simulation of this equation is of prohibitive
computational expense.  Nonetheless, my dissertation will quantify this
expense.)

\item One may also hope that this equation could be written in terms of
\plainidx{special functions}.  Our research on this matter, and our contacts
with experts in Brownian bridges, has not produced a special-function
expression.

\item Although one may not simplify all interaction pairs, one may extract the
strongest interaction pairs (namely, two-cycles) and simplify those.  This is
discussed in section \ref{sec:n2}.

\end{itemize}

See also section \ref{sec:future} where future research is sketched.

%% ================================================================
\newpage
\section{Simple model with two-cycle jump-pair interactions}
\label{sec:n2}

For the simplified two-cycle-interaction model, unlike the fully interacting
model of section \ref{sec:tji}, one obtains expressions for the pressure,
critical density, and critical temperature for the weakly interacting Bose gas.
These appear as perturbations to the known expressions for the ideal gas.

%% ----------------------------------------------------------------
\subsection{Motivation}

Equation \ref{eqn:connect_V} is an expression for the jump-pair interaction.
In section \ref{sec:tji}, we reviewed the difficulties in either estimating or
simplifying this expression in general.

The pair-jump potential is the collision probability for two bridges, one
running from $\vecx_i$ to $\vecxpii$ and another running from $\vecx_j$ to
$\vecxpij$, in time $4\beta$; one expects this probability to be highest for
\plainidx{two-cycles}.  (In figure \ref{fig:2w_to_1w}, one expects increased
collision probability if $\vecxpii=\vecx_j$.)  Thus, one may choose to neglect
all $V_{ij}$ terms except those involving two-cycles.

\begin{notn}
The shorthand notation
$$
	i \twocyc j
$$
means that $i$ and $j$ participate in a two-cycle under $\pi$, i.e. $\pi(i)=j$,
$\pi(j)=i$, and $i < j$.
\end{notn}

Then equation
\ref{eqn:connect_H} becomes
\begin{equation}
\label{eqn:two_N2_hams}
\begin{aligned}
	H(\vecX,\pi) &= \sum_{i=1}^N \frac{1}{4\beta}\|\vecx_i-\vecxpii\|^2
		+ \sum_{1\le i<j \le N}
			V(\vecx_i, \vecxpii, \vecx_j, \vecxpij) \\
	\approx
	\Htilde_P(\vecX,\pi) &= \sum_{i=1}^N \frac{1}{4\beta}\|\vecx_i-\vecxpii\|^2
		+ \sum_{i \twocyc j} V(\vecx_i, \vecxpii, \vecx_j, \vecxpij).
\end{aligned}
\end{equation}

An unpublished computation of Ueltschi and Betz shows that, for two-cycles, the
jump-pair interaction (equation \ref{eqn:connect_V}) simplifies significantly
to
\begin{equation}
\label{eqn:Vxyyx_2a}
	V(\vecx_i,\vecxpii,\vecxpii,\vecx_i)
	= \frac{2a}{\|\vecx_i-\vecxpii\|} + O(a^2),
\end{equation}
where $a$ is the radius of the interparticle hard-core potential $U$.

The key point is that the Brownian bridges of equation \ref{eqn:connect_V} have
been simplified out completely for this two-cycle-interaction model.  To
estimate expectations of the random variable $\varrho_{m,n}$ (section
\ref{sec:model_defns}), one may simply use \plainidx{Markov-Chain Monte Carlo
techniques} to sample random permutations with weights derived from equation
\ref{eqn:Halpha}.  (The details of such MCMC simulations will be explained in
my dissertation.)

%% ----------------------------------------------------------------
\subsection{Hamiltonian with $N_2(\pi)$}
\label{sec:n2_ham}

Equation \ref{eqn:Vxyyx_2a} is quite adequate for use in Monte Carlo
simulations.  A further modification is made, though, which facilitates the
computation of thermodynamic properties of the two-cycle-interaction model in
the following section.

Letting $N_2(\pi)$ denote the number of \plainidx{two-cycles} in the
permutation $\pi$, one would like to have a simple Hamiltonian\index{$\alpha$}
of the form
\begin{align}
\label{eqn:Halpha}
	\Halpha_P(\vecX,\pi) &= \sum_{i=1}^N \frac{1}{4\beta}\|\vecx_i-\vecxpii\|^2
		+ \alpha N_2(\pi).
\end{align}
It remains to connect the old parameter $a$ with the new parameter $\alpha$.
Note that the distance dependence in equation \ref{eqn:Vxyyx_2a} needs to be
averaged out:  in equation \ref{eqn:Halpha}, all two-cycles acquire the same
weight $\alpha$ regardless of $\|\vecx_i-\vecxpii\|$.  

\begin{prop}
\label{prop:alpha_a}
The first-order approximation
\begin{align}
\label{eqn:alpha_a}
	\alpha = \left(\frac{8}{\pi\beta}\right)^{1/2} a + O(a^2)
\end{align}
is valid for computing expectations.
\begin{proof}
It may seem odd that the distances $\|\vecx_i-\vecx_j\|$ are replaced by a
common $\alpha$.  Remember, however, that (see section \ref{sec:model_defns})
we are averaging over particle positions $\vecx_i,\vecx_j\in\Lambda$.  Equating
the two Hamiltonians in \ref{eqn:two_N2_hams}, we have
\begin{align*}
	\int_{\Lambda^N}d\vecX
	&
		\exp\left\{
			-\frac{1}{4\beta} \sum_{i=1}^N \|\vecx_i-\vecxpii\|^2
		\right\}
		\exp\left\{
			-\sum_{i \twocyc j} V(\vecx_i, \vecxpii, \vecx_j, \vecxpij)
		\right\} \\
	&=
	\int_{\Lambda^N}d\vecX
		\exp\left\{
			-\frac{1}{4\beta} \sum_{i=1}^N \|\vecx_i-\vecxpii\|^2
		\right\}
		\exp\left\{
			-\sum_{i \twocyc j} \alpha
		\right\}
	\\
%
	\int_{\Lambda^N}d\vecX
	&
		\prod_{i=1}^N
		\exp\left\{
			-\frac{1}{4\beta} \|\vecx_i-\vecxpii\|^2
		\right\}
		\prod_{i \twocyc j}
		\exp\left\{
			-V(\vecx_i, \vecxpii, \vecx_j, \vecxpij)
		\right\} \\
	&=
	\int_{\Lambda^N}d\vecX
		\prod_{i=1}^N
		\exp\left\{
			-\frac{1}{4\beta} \|\vecx_i-\vecxpii\|^2
		\right\}
		\prod_{i \twocyc j}
		\exp\left\{
			\alpha
		\right\}.
\end{align*}
By Tonelli's theorem (given the positivity of the integrand) we may iterate the
integrals as we wish.  All $i$'s which do not participate in two-cycles factor
out, canceling on both sides.  Since $\pi(i)=j$ and $\pi(j)=i$, and
substituting expression \ref{eqn:Vxyyx_2a} for $V$, we are left with
\begin{align*}
	&
	\prod_{i \twocyc j}
	\int_{\Lambda} \int_{\Lambda} \,d\vecx_i \,d\vecx_j \,
	\exp\left\{ -\frac{1}{2\beta} \|\vecx_i-\vecx_j\|^2 \right\}
	\exp\left\{ -\frac{2a}{\|\vecx_i-\vecx_j\|} \right\} \\
	&=
	\prod_{i \twocyc j}
	\int_{\Lambda} \int_{\Lambda} \,d\vecx_i \,d\vecx_j \,
	\exp\left\{ -\frac{1}{2\beta} \|\vecx_i-\vecx_j\|^2 \right\}
	\exp\left\{ -\alpha\right\}.
\end{align*}
It suffices to apply the same technique to each $i \twocyc j$ pair.
Let
$$
	\vecr_{ij}=\vecx_i-\vecx_j.
$$
As in section \ref{sec:2w_to_1w}, the change-of-variables Jacobian has
determinant with absolute value 1, so without need for correction factors we
write
\begin{align*}
	\int_{\Lambda} \int_{\Lambda} \,d\vecx_i \,d\vecr_{ij} \,
	\exp\left\{ -\frac{1}{2\beta} \|\vecr_{ij}\|^2 \right\}
	\exp\left\{ -\frac{2a}{\|\vecr_{ij}\|} \right\}
	&=
	\int_{\Lambda} \int_{\Lambda} \,d\vecx_i \,d\vecr_{ij} \,
	\exp\left\{ -\frac{1}{2\beta} \|\vecr_{ij}\|^2 \right\}
	\exp\left\{ -\alpha\right\}.
\end{align*}
Taylor-expanding $\exp(-2a/\|\vecr_{ij}\|)$ and $\exp(-\alpha)$ to first order,
since we are interested in weak interactions, gives us
\begin{align*}
	\int_{\Lambda} \int_{\Lambda} \,d\vecx_i \,d\vecr_{ij} \,
	\exp\left\{ -\frac{1}{2\beta} \|\vecr_{ij}\|^2 \right\}
	\left( 1 -\frac{2a}{\|\vecr_{ij}\|} \right)
	&=
	\int_{\Lambda} \int_{\Lambda} \,d\vecx_i \,d\vecr_{ij} \,
	\exp\left\{ -\frac{1}{2\beta} \|\vecr_{ij}\|^2 \right\}
	\left(1 - \alpha\right).
\end{align*}
After distributing the exponentials over the differences, the first products
are Gaussians which cancel out on both sides.  What is left is
\begin{align*}
	2a
	\int_{\Lambda} \int_{\Lambda} \,d\vecx_i \,d\vecr_{ij} \,
	\frac{ \exp\left\{ -\frac{1}{2\beta} \|\vecr_{ij}\|^2 \right\} }
		{\|\vecr_{ij}\|}
	&=
	\alpha
	\int_{\Lambda} \int_{\Lambda} \,d\vecx_i \,d\vecr_{ij} \,
	\exp\left\{ -\frac{1}{2\beta} \|\vecr_{ij}\|^2 \right\}.
\end{align*}
The integral over $d\vecx_i$ gives the volume of $\Lambda$ on both sides, which
cancels; for $d\vecr_{ij}$, we pass from $\Lambda$ to $\R^3$, ignoring boundary
corrections which will disappear in the infinite-volume limit.  The right-hand
integral over $d\vecr_{ij}$ is $(2\pi\beta)^{3/2}$ (section
\ref{sec:gaussians}).  Then
\begin{align*}
	\alpha
	&= \frac{2a}{ (2\pi\beta)^{3/2}} \,
	\int_{\R^3} \,d\vecr \,
	\frac{ \exp\left\{ -\frac{1}{2\beta} \|\vecr\|^2 \right\} }
		{\|\vecr\|} \\
%
	&= \frac{2a}{ (2\pi\beta)^{3/2}} \,
	\int_{\theta=0}^{\theta=2\pi}
	\,d\theta
	\int_{\phi=0}^{\phi=\pi}
	\,d\phi
	\sin\phi
	\int_{r=0}^{r=\infty}
	\,dr
	\,\,r
	e^{-r^2/2\beta} \\
%
	&= \frac{8\pi a}{ (2\pi\beta)^{3/2}} \,
	\int_{r=0}^{r=\infty}
	\,dr \, r \,
	e^{-r^2/2\beta} \\
%
	&= \frac{8\pi \beta a}{ (2\pi\beta)^{3/2}}
	= \left( \frac{8}{\pi\beta} \right)^{1/2} \, a.
%
\end{align*}
(The $r$ integral really should start at $r=2a$ rather than $r=0$, since we
have a hard-core potential which prohibits $\vecx_i$ from being within $2a$ of
$\vecx_j$.  This would introduce a factor of $e^{-2a^2/\beta}$ into the above.
In the small-$a$ regime, though, we set this to 1 without remorse.)
\end{proof}
\end{prop}

%% ----------------------------------------------------------------
\subsection{Pressure}

Here we find the \plainidx{pressure} of the two-cycle-interaction model and
compare it to the pressure of the \plainidx{ideal gas}.  The latter is
\cite{bib:bu}:
\begin{equation}
\label{eqn:p_0}
	p^{(0)}(\beta,\mu) =
	-\frac{1}{\beta} \int_{\R^d}
	\log\left(1 - e^{-\beta \left(4\pi^2\|\veck\|^2 -\mu\right)}\right) d\veck.
\end{equation}
(The \bfidx{chemical potential} $\mu$ is defined in section \ref{sec:mu}.

\begin{prop}
The pressure of the two-cycle-interaction model is
\begin{equation}
\label{eqn:p_alpha}
	p^{(\alpha)}(\beta, \mu) =
	p^{(0)}(\beta,\mu)
	%-\frac{1}{\beta} \int_{\R^d}
	%\log\left(1 - e^{-\beta \left(4\pi^2 \|\veck\|^2-\mu\right)}\right) d\veck
	-
	\frac{e^{2\beta\mu}}{2^{11/2}\pi^{3/2}\beta^{5/2}}(1-e^{-\alpha}).
\end{equation}
\begin{proof}
% XXX prove me.
See section 7.1 of \cite{bib:bu}, or section 5 of \cite{bib:qmath}.  The proof
involves a lengthy decomposition into Fourier modes; I will exposit it as time
permits.
\end{proof}
\end{prop}

% XXX footnote about $\xi(\vecx)$ and $\veps(\veck)$.

%%% ----------------------------------------------------------------
%\subsection{Free energy}
%
%XXX \bfidx{free energy} (\plainidx{Legendre transform}):
%$$
%	f^{(\alpha)}(\beta,\rho) = \sup_\mu\left[
%		\rho \mu - p^{(\alpha)}(\beta,\mu)
%		\right].
%$$

%% ----------------------------------------------------------------
\subsection{Conjecture for the critical density}
\label{sec:N2_rho_c}

Recall from section \ref{sec:rho} (see also equation \ref{eqn:ideal_rho_c})
that the \plainidx{critical density} for the ideal Bose gas is
\begin{equation}
\label{eqn:rho_c_0}
	\rho_c^{(0)} = \int_{\R^d} \frac{d\veck}{e^{4\beta\pi^2|\veck|^2}-1}.
\end{equation}

As I continue to work on this random-cycle project, I certainly have more to
learn about thermodynamics.  In particular, one might inquire why the partial
derivative in the following proposition defines the critical density.
(\cite{bib:huang} approaches the subject from a different perspective than
Ueltschi does.) The best I can say at present is:

\begin{itemize}

\item The chemical potential is defined to be change in energy per additional
particle, with fixed volume and entropy, i.e. $\mu = (\partial E/\partial
N)_{S,V}$.

\item Particles in the ground state (condensed particles) contribute nothing to
the pressure.

\end{itemize}

\begin{prop}
\label{prop:rho_c_alpha}
The critical density for the two-cycle-interaction model is
\begin{equation}
\label{eqn:rho_c_alpha}
	\rho_c^{(\alpha)} = \frac{\partial p^{(\alpha)}}{\partial\mu}
		\Bigg|_{\mu=0-}
	= \rho_c^{(0)} - \frac{(1-e^{-\alpha})}{2^{9/2}\pi^{3/2}\beta^{3/2}}.
\end{equation}
\begin{proof}
Differentiating equation \ref{eqn:rho_c_alpha} through the integral sign, we
obtain
\begin{align*}
	%\frac{\partial p}{\partial \mu}\Big|_{\mu=0^-}
	\frac{\partial p}{\partial \mu}
&=
	-\frac{1}{\beta}
	\frac{\partial}{\partial \mu} \left[
	\int_{\R^d}
	\log\left(1 - e^{-4\pi^2\beta\|\veck\|^2} e^{\beta\mu}\right) d\veck
	\right]
	-
	\frac{\partial}{\partial \mu} \left[
	\frac{e^{2\beta\mu}}{2^{11/2}\pi^{3/2}\beta^{5/2}}(1-e^{-\alpha})
	\right] \\
&=
	-\frac{1}{\beta}
	\int_{\R^d}
	\frac
		{ -\beta e^{-4\pi^2\beta\|\veck\|^2} e^{\beta\mu} }
		{ 1 - e^{-4\pi^2\beta\|\veck\|^2} e^{\beta\mu} }
	d\veck
	-
	\frac{2\beta e^{2\beta\mu}}{2^{11/2}\pi^{3/2}\beta^{5/2}}(1-e^{-\alpha}) \\
%
	\frac{\partial p}{\partial \mu}\Big|_{\mu=0^-}
&=
	\int_{\R^d}
	\frac
		{ e^{-4\pi^2\beta\|\veck\|^2}  }
		{ 1 - e^{-4\pi^2\beta\|\veck\|^2}  }
	d\veck
	- \frac{(1-e^{-\alpha})}{2^{9/2}\pi^{3/2}\beta^{3/2}} \\
%
&=
	\int_{\R^d} \frac{d\veck}{e^{4\pi^2\beta\|\veck\|^2} - 1}
	- \frac{(1-e^{-\alpha})}{2^{9/2}\pi^{3/2}\beta^{3/2}}.
\end{align*}

As is familiar from real analysis, differentiation through the integral sign is
justified since
(i) the integrand is $\veck$-integrable for all $\mu<0$, and
(ii) the integrand has continuous partial derivative with respect to $\mu$.
\end{proof}
\end{prop}

\begin{conjecture}
Infinite cycles\index{infinite cycles} occur for densities above $\rho_c$.
\end{conjecture}
Note that the \plainidx{critical density} can be computed; what is not proved
is that it corresponds to \plainidx{infinite cycles}.  The following, weaker
result has been proved in \cite{bib:bu}.

\begin{thm}
For all $b<1$,
\begin{align*}
	\lim_{V\to\infty}
	\bbE_{\Lambda,N}\left(\varrho_{N^b,N}\right)
	\ge \rho - \frac{4}{(1+e^{-\alpha})^2} \, \rho_c^{(0)}
\end{align*}
where the limit is taken with fixed ratio $\rho=N/V$.
\end{thm}

Compare this to equation \ref{eqn:thm_1_rho} of section
\ref{sec:inf_cyc_existence}:  the $\rho_c$ in $\rho-\rho_c$ acquires the scale
factor $4/(1+e^{-\alpha})^2$, and we have only a lower bound.

\emph{Late note:}  A theorem of this type has been proved in the recent paper
\cite{bib:bu08}.

%% ----------------------------------------------------------------
\subsection{Shift in critical temperature}
\label{sec:delta_Tc}

\begin{prop}
\label{prop:delta_Tc}
For the weakly interacting two-cycle model,
$$
	\frac{T_c^{(a)}-T_c^{(0)}}{T_c^{(0)}}
	\approx 0.37 \rho^{1/3}a.
$$
\begin{proof}
Recall from section \ref{sec:critical_temperature}, and in particular figure
\ref{fig:critical_manifold}, that the critical region\index{critical
manifold}\index{critical temperature}\index{critical density} may be seen as a
manifold in $(\rho,\beta,a)$ space.  We will use the result
$$
	\frac{\D a}{\D\rho}
	\frac{\D \rho}{\D\beta}
	\frac{\D \beta}{\D a} = -1
$$
from proposition \ref{prop:triple_partial_product}.  (Since we are working on
the critical manifold, we take $\rho$ and $\beta$ to mean $\rho_c^{(a)}$ and
$\beta_c^{(a)}$, respectively.)  From equation \ref{eqn:ideal_rho_c_zeta}, the
\plainidx{critical line}, with $a=0$, has equation
\begin{align}
\label{eqn:critical_line}
	\rho_c^{(0)}
	{\beta_c^{(0)}}^{3/2}
	&= \frac{\zeta(3/2)}{(4\pi)^{3/2}}
\end{align}
where $\zeta(z)$ is the \plainidx{Riemann zeta function}.

First, we use proposition \ref{prop:rho_c_alpha} for the critical density, then
Taylor-expand in the small parameter $\alpha$ and use equation
\ref{eqn:alpha_a} to relate $\alpha$ and $a$:
\begin{align*}
	\rho_c^{(a)} &= \rho_c^{(0)} - \frac{1-e^{-\alpha}}{(8\pi\beta)^{3/2}} \\
%
	\frac{\rho_c^{(a)} - \rho_c^{(0)}}{\rho_c^{(0)}}
	&\approx
	\frac{-\alpha}{(8\pi\beta)^{3/2}}
	\frac {(4\pi\beta)^{3/2}} {\zeta(3/2)}
%
	=
	\frac{-8^{1/2} a}{\pi^{1/2} \beta^{1/2} (8\pi\beta)^{3/2}}
	\frac {(4\pi\beta)^{3/2}} {\zeta(3/2)}
%
	=
	\frac {- a} {\pi^{1/2} \beta^{1/2} \zeta(3/2)}.
\end{align*}
Letting
\begin{align}
\label{eqn:shift_prop_b}
	b &= \frac{1}{\zeta(3/2) \pi^{1/2}}
\end{align}
for brevity, we have
\begin{align*}
	\frac{\rho_c^{(a)}-\rho_c^{(0)}}{\rho_c^{(0)}}
	&= -\frac{ b a }{\beta^{1/2}}.
\end{align*}
Then
\begin{align*}
	a &= \frac{-\rho_c^{(a)} \beta^{1/2}}{\rho_c^{(0)} b}
		+ \frac{\beta^{1/2}}{b}
\end{align*}
and
\begin{align*}
	\frac{\D a}{\D\rho} &= \frac{-\beta^{1/2}}{\rho_c^{(0)} b}.
\end{align*}

Second, using equation \ref{eqn:critical_line},
\begin{align*}
	\rho_c^{(0)} &= \frac{\zeta(3/2)}{(4\pi\beta)^{3/2}}
\end{align*}
and so
\begin{align*}
	\frac{\D \rho_c^{(0)}}{\D\beta} &= \frac{-\zeta(3/2)}{(4\pi\beta)^{3/2}}.
\end{align*}
For small $a$ (see figure \ref{fig:critical_manifold} for intuition) we
use the approximation
\begin{align*}
	\frac{\D \rho_c^{(a)}}{\D\beta} &\approx \frac{\D \rho_c^{(0)}}{\D\beta}.
\end{align*}

Third, from section \ref{sec:critical_temperature} we expect
\begin{align*}
	\frac{T_c^{(a)}-T_c^{(0)}}{T_c^{(0)}} &= c \rho^{1/3} a.
\end{align*}
Putting $\beta=1/T$ gives
\begin{align*}
	\frac{\beta_c^{(a)}-\beta_c^{(0)}}{\beta_c^{(a)}}
	&= -c \rho^{1/3} a
\end{align*}
which we approximate by
\begin{align*}
	\frac{\beta_c^{(a)}-\beta_c^{(0)}}{\beta_c^{(0)}}
	&= -c \rho^{1/3} a.
\end{align*}
Solving for $\beta_c^{(a)}$ gives
\begin{align*}
	\beta_c^{(a)} &= \beta_c^{(0)} - \beta_c^{(0)} c \rho^{1/3} a
\end{align*}
and thus
\begin{align*}
	\frac{\D \beta}{\D a}
	&= -\beta_c^{(0)} c \rho^{1/3}.
\end{align*}

Combining the product of all three partial derivatives and using proposition
\ref{prop:triple_partial_product}, we have
\begin{align*}
	\left( \frac{\beta^{1/2}}{\rho_c^{(0)} b} \right)
	\left( \frac{\zeta(3/2)}{(4\pi\beta)^{3/2}} \right)
	\left( \beta_c^{(0)} c \rho^{1/3} \right)
	&= 1.
\end{align*}
Solving for $c$ gives
\begin{align*}
	c &=
	\frac{
		\rho_c^{(0)} \rho^{-1/3} \beta^{5/2}
	}{
		\beta_c^{(0)} \beta^{1/2}
	}
	\frac{ 2b \, (4\pi)^{3/2} }{ 3 \,\zeta(3/2) }.
\end{align*}
Approximating $\beta_c^{(0)}$ and $\rho_c^{(0)}$ by $\beta_c^{(a)}$ and
$\rho_c^{(a)}$, we have
\begin{align*}
	c &=
	\left( \rho^{2/3} \beta \right)
	\frac{ 2b \, (4\pi)^{3/2} }{ 3 \,\zeta(3/2) }
	=
	\left( \rho \beta^{3/2} \right)^{2/3}
	\frac{ 2b \, (4\pi)^{3/2} }{ 3 \,\zeta(3/2) }.
\end{align*}
From equation \ref{eqn:critical_line} we know that
\begin{align*}
	\rho_c^{(0)} \beta^{3/2} &= \frac{\zeta(3/2)}{(4\pi)^{3/2}}
\end{align*}
which we apply to $\rho_c^{(a)}$ as well.  Then
\begin{align*}
	c &=
	\left( \frac{\zeta(3/2)}{(4\pi)^{3/2}} \right)^{2/3}
	\frac{ 2b \, (4\pi)^{3/2} }{ 3 \,\zeta(3/2) }
	=
	\left( \frac{\zeta(3/2)^{2/3}}{(4\pi)} \right)
	\frac{ 2b \, (4\pi)^{3/2} }{ 3 \,\zeta(3/2) }
	=
	\frac{ 4b \, \pi^{1/2} }{ 3 \,\zeta(3/2)^{1/3} }.
\end{align*}
Since
\begin{align*}
	b &= \frac{1}{\zeta(3/2) \pi^{1/2}},
\end{align*}
\begin{align*}
	c &=
	\frac{4}{3} \zeta(3/2)^{-4/3} \approx 0.37.
\end{align*}
\end{proof}
\end{prop}

\begin{rem*}
This result applies for the two-cycle model, where the only permutation jumps
that interact are those which participate in two-cycles.  When longer cycles
are included, the shift in critical temperature is expected to be more
pronounced.  Thus, this result provides a rough lower bound on the true
constant $c$, which from other methods discussed in section
\ref{sec:critical_temperature} is believed to be approximately 1.3.  Further
work is needed (see the following section) before the random-cycle model can be
used to improve on the latter estimate.
\end{rem*}

%For these reasons, I have shown $\D\rho_c^{(a)}/\D a$ and $\D\beta_c^{(a)}/ \D
%a$ to be negative in figure \ref{fig:critical_manifold}.

%% ================================================================
\newpage
\section{Future work}
\label{sec:future}

%% ----------------------------------------------------------------
\subsection{Theoretical directions}
\label{sec:future_theory}

My research to date has shown that the full jump-pair interaction (see equation
\ref{eqn:connect_V} and section \ref{sec:tji}) is too computationally intensive
for practical use.  Ueltschi and Betz have simplified the Brownian bridges out
of the jump-pair interaction for the particular case of two-cycles (equation
\ref{eqn:Vxyyx_2a}); they also have some unpublished work on approximating the
full interaction (to first order in the scattering length $a$) as a Riemann
integral.  I have found that this, too, while better than the Brownian bridges,
is still computationally expensive.  If this integral could be represented as a
special function, then the full jump-pair-interacting case could be examined
experimentally.

Alternatively, notice that, for the two-cycle-interaction model, the distance
dependence in equation \ref{eqn:Vxyyx_2a} was averaged out in equation
\ref{eqn:Halpha}.  Since we are interested only in expected values of
$\rho_{m,n}$, the details of a particular interaction are of minor importance;
only the averaged behavior is of final interest.  Perhaps equation
\ref{eqn:connect_V} may be handled in a manner similar to proposition
\ref{prop:alpha_a}.

A very recent paper of Betz and Ueltschi \cite{bib:bu08} extends the two-cycle
Hamiltonian of section \ref{sec:n2} (equation \ref{eqn:Halpha}) to the form
\begin{align}
\label{eqn:Halphas}
	\HP(\vecX,\pi)
	&= \sum_{i=1}^N \frac{1}{4\beta}\|\vecx_i-\vecxpii\|^2
	+ \sum_{\ell=1}^N \alpha_\ell N_\ell(\pi)
\end{align}
where $N_\ell(\pi)$ counts the number of $\ell$-cycles in the permutation
$\pi$.  I have yet to evaluate how this recent development affects my own work.

%% ----------------------------------------------------------------
\subsection{Numerical directions}
\label{sec:future_numerics}

Simulations currently underway use the two-cycle-interaction model, with points
on a cubic unit lattice.  One would like to vary the positions of the points as
well, in order to simulate the point-process-configuration model.  (See the
last paragraph of section \ref{sec:model_defns} for a description of the
probability measures for these two models.)

Note that the choice between the lattice-configuration and the
point-process-configuration models is independent of the choice between the
two-cycle-interaction and full jump-pair-interaction models.

If the theoretical work outlined in section \ref{sec:future_theory} produces an
efficient expression for the full jump-pair interaction, then it will certainly
be examined experimentially.

%% ----------------------------------------------------------------
\subsection{Statistical directions}
\label{sec:future_stats}

Markov-chain Monte Carlo simulations compute sample means of $\varrho_{m,n}$
(see section \ref{sec:model_defns}) for the lattice and point-process
probability measures --- producing one mean for a given choice of lattice size
$N$, particle density $\rho$, inverse temperature $\beta$, and interaction
parameter $a$.

We will analyze the error of the resulting sample means.  For a large number of
trials, one expects a central-limit distribution for the estimated values of
$\varrho_{m,n}$; we also desire to have a practical estimator for the variance
of the sample mean.

To approach the infinite-volume limit in $N$, one needs to simulate larger and
larger lattice sizes, and then do \plainidx{finite-size scaling} on the
results.  This analysis, involving a regression of a carefully chosen form with
due sensitivity to sampling errors, will be a central part of my dissertation
work.

%% ================================================================
\newpage
\section{Acknowledgements}

Bill Faris' notes on the Feynman-Kac formula have proved invaluable, as has
Janek Wehr's spring 2008 course in stochastic differential equations.  Tom
Kennedy made it possible for me to participate in this research project, and
has been a patient and learned sounding board throughout.  Daniel Ueltschi has
shared with me his extensive knowledge and his enthusiasm, and has created an
environment of great mutual benefit.  Each of these gentlemen deserves my
gratitude; this work would not have been possible without any of them.

Additionally, much of the work here was done during the Summer School for
Current Topics in Mathematical Physics at the Erwin Schr\"odinger Institute in
July 2008.  The Institute provided a productive and fruitful atmosphere; the
University of Arizona Department of Mathematics helped make the visit possible.

%% ================================================================
\appendix
\newpage
\section{Quantum mechanics}

%For information on \plainidx{quantum mechanics}, see \cite{bib:griffiths},
%\cite{bib:nc}, and \cite{bib:sakurai}; see \cite{bib:huang}, \cite{bib:widom},
%and \cite{bib:pathria} for \plainidx{statistical mechanics}.  The topics
%presented here are a digestion of these sources, with unified vocabulary and
%notation, serving as a single point of reference for the main body of the
%paper.

For information on \plainidx{quantum mechanics}, see \cite{bib:griffiths},
\cite{bib:nc}, and \cite{bib:sakurai}; see \cite{bib:huang}, \cite{bib:widom},
and \cite{bib:pathria} for \plainidx{statistical mechanics}.  The topics
presented here are a digestion of these sources, with unified vocabulary and
notation, serving a single purpose:  so that I may explain to myself why
$\TrebH$ is of such central importance.

% \ref{sec:canon_ens_rho}

%% ----------------------------------------------------------------
\subsection{Postulates of quantum mechanics}
\label{sec:postulates}

I take elementary quantum mechanics as given; in particular, I will use bra-ket
notation without explanation.  (I should point out that I take the complex
inner product to be linear on the right, in accordance with
mathematical-physics convention.) However, I will recall the basic postulates
of quantum mechanics in order to smoothly elucidate the usefulness of the
otherwise mysterious density matrix.  In turn, density matrices are needed to
formulate quantum statistical mechanics in section \ref{sec:canon_ens_rho}.

\begin{itemize}

\item[(1)] To any isolated physical system is associated a separable complex
Hilbert space $\msH$, called a \emphidx{state space}, with inner product
$\braket{\phi}{\psi}$.  States (\plainidx{pure states}; see section
\ref{sec:pure_and_mixed}) of the physical system are described by unit vectors
(often called \emphidx{state vectors} or \emphidx{wave functions}) in $\msH$,
i.e. vectors $\psi$ such that $\braket{\psi}{\psi}=1$.  Vectors that are scalar
multiples of one another are considered to be equivalent.  (If $\psi_1 =
c\psi_2$, then, since $\psi_1$ and $\psi_2$ have norm 1, $|c|=1$.  That is,
$c=e^{i\theta}$ for some real $\theta$.  Thus arises the saying that state
vectors are distinct up to a \emphidx{phase factor}.) The state space of a
composite system is the \plainidx{tensor product} of the state spaces of the
component systems.

\item[(2)] State vectors evolve in time via \emphidx{unitary transformations}.
Specifically, the unitary operator is described by the \emphidx{Schr\"odinger
equation}:
$$
	i \hbar \frac{\D\psi(t,\vecx)}{\Dt}
	= \frac{-\hbar^2}{2m} \nabla^2 \psi(t,\vecx) + V(t,\vecx) \psi(t,\vecx)
	= H \psi(t,\vecx)
$$
i.e.
$$
	\frac{\D\psi}{\Dt}
	= \frac{-i H}{\hbar} \psi
$$
where $H$ is the Hermitian operator
$$
	H = \frac{-\hbar^2}{2m} \nabla^2 + V.
$$
The \plainidx{wave function} varies in position and time; the
\emphidx{potential function}, conventionally written with the letter $V$, may
also vary in position and time.  The potential represents external interactions
with the particle.

We can write down a solution
\begin{align}
\label{eqn:SE}
	\frac{d\psi}{dt} &= \frac{-i}{\hbar} \, H \psi \\
	\psi(t) &= e^{-iHt/\hbar} \psi(0).
\end{align}
One checks this solution by differentiating:
\begin{align*}
	\frac{d\psi}{dt} &= \frac{-iH}{\hbar} \, e^{-iHt/\hbar} \psi(0)
	= \frac{-iH}{\hbar} \, \psi(t)
	= \frac{-i}{\hbar} \, H \psi(t).
\end{align*}
As is well known, the imaginary exponential of a hermitian matrix is
\plainidx{unitary}.  The importance of this is that unitary matrices $U$ are
\plainidx{norm-preserving}.  That is, if $\|\psi\|=1$ then $\|U\psi\|$ is still
1; probability-density functions remain probability-density functions as they
evolve in time.  This fact is called the \emphidx{conservation of probability}.

%Real spectrum.  Existence of orthonormal basis of eigenvectors.
%Why they \emph{must} be orthonormal.

\item[(3)] An observable $\mcO$ corresponds to a densely defined Hermitian
operator $A:\msH\to \msH$.  Since $A$ is Hermitian, $\msH$ has an
\plainidx{orthonormal basis} $\{\phi_j\}$ of eigenvectors of $A$.  Since $\msH$
is \plainidx{separable}, the basis is a countable set.  In the case that $A$
has \plainidx{discrete spectrum}, the eigenstates of $A$ have respective
(necessarily real) eigenvalues $a_j$ and $A$ has a \emphidx{spectral
decomposition}
$$
	A = \sum a_j \ketbra{\phi_j}
$$
(The operators $P_j = \ketbra{\phi_j}$ are known as \emphidx{projection
operators}.) A (necessarily square-summable) wave function $\psi$ is a linear
combination of the eigenbasis:
$$
	\psi = \sum_{j} c_j \phi_j
	\qqwith
	c_j = \braket{\phi_j}{\psi}
	\qqand
	\sum_j |c_j|^2 = 1.
$$
A measurement of $\mcO$ has as its possible outcomes the numbers $a_j$ with
probabilities
$$
	\bbP(\mcO=a_j) = |\braket{\phi_j}{\psi}|^2.
$$
The expected value of $A$ is
\begin{align*}
	\expn{A} &= \sum_j a_j \bbP(\mcO=a_j)
	= \sum_j a_j |\braket{\phi_j}{\psi}|^2
	= \sum_j a_j \braket{\psi}{\phi_j} \braket{\phi_j}{\psi} \\
	&= \sum_j a_j
		\bra{\psi}\; \bigg(\ketbra{\phi_j}\bigg) \; \ket{\psi}
	= \sum_j \bra{\psi}\;
		\bigg(a_j \ketbra{\phi_j}\bigg)
		\; \ket{\psi} \\
	&= \braCket{\psi}{A}{\psi}.
\end{align*}
This is the \plainidx{average value} or \plainidx{mean} of all the
observations, averaged over repeated measurements on identical systems.

\item[(4)] Immediately after measurement of an observable $\mcO$, the state of
the system is described \emph{only} by an eigenstate $\phi_j$ of $A$.  This is
called the \emphidx{collapse of the wave function}.  The evolution of the
single-particle system thereafter is described by the Schr\"odinger equation,
with new initial conditions.

\end{itemize}

%% ----------------------------------------------------------------
\subsection{Pure and mixed states}
\label{sec:pure_and_mixed}

The terms \emphidx{pure state} and \emphidx{mixed state} are horribly
misleading.  Nonetheless, they are entrenched in the literature and cannot be
dispensed with.  A pure state is simply a state of a quantum-mechanical system:
the adjective \emph{pure} is superfluous.  A mixed state not a single state at
all:  it is a statistical description of a \emph{set} of many states.

\begin{defn}
An \bfidx{ensemble} is a list of $n$ distict state vectors $\{\psi_1, \ldots,
\psi_n\}$ along with respective probabilities $\{p_1, \ldots, p_n\}$ with $0
\le p_j \le 1$ and $\sum_{j=1}^n p_j = 1$.  (One may think of an ensemble as a
probability mass function.)
\end{defn}

\begin{defn}
A \bfidx{mixed state} is an ensemble.  A \bfidx{pure state} is (with abuse of
notation) either a single state, or an ensemble with $n=1$, or an ensemble with
$n \ge 1$ but only one $p_k=1$ and the remaining $p_j=0$ for $j\ne k$.
\end{defn}

\begin{rem*}
Pure and mixed states are not to be confused with the following:
\begin{itemize}
\item A quantum-mechanical state is a linear combination of eigenstates.
	Such states are pure, not mixed.
\item Entangled and separable states arise in the theory of multi-particle
	systems, e.g. the tensor product $\mcH_1 \tensor \mcH_2$ of Hilbert spaces
	$\mcH_1$ and $\mcH_2$.  An entangled state is an indecomposable tensor in
	$\mcH_1 \tensor \mcH_2$; a separable state is a decomposable tensor.
\end{itemize}
\end{rem*}

%% ----------------------------------------------------------------
\subsection{Density matrices}

Density matrices are a notational convenience for working with ensembles (mixed
states) of quantum systems.  I will develop\footnote{This development follows,
in part, that in \texttt{http://electron6.phys.utk.edu/QM1}, retrieved on June
12, 2008.} density matrices in three steps:  (1)  In this section, I will show
how to construct a density matrix from an ensemble.  (2) In section
\ref{sec:pure_dm}, I will show how density matrices serve to notationally
encapsulate the postulates of quantum mechanics as described in section
\ref{sec:postulates}.  (3) In section \ref{sec:mixed_dm} I will show how this
notation permits one to extend those postulates to the mixed-state case.

\begin{defn}
A \bfidx{density matrix} is a positive-definite matrix with trace 1.
\end{defn}

\begin{prop}
\label{prop:rho_from_ens}
A density matrix may be obtained from an ensemble by
$$
	\rho = \sum_{j=1}^n p_j \ketbra{\psi_j}.
$$
\begin{proof}
Let $\rho$ be such a matrix.
Fix an orthonormal basis $\{\phi_i\}$ for $\msH$.
I will use the characterization that a matrix $A$ is positive-definite iff it
is Hermitian and $\braCket\xi A \xi > 0$ for all non-zero $\xi$.  To see that
$\rho$ is Hermitian, write
\begin{align*}
	\rho^* &= \left(\sum_{j=1}^n p_j \ketbra{\psi_j}\right)^*
	= \sum_{j=1}^n p_j \left(\ketbra{\psi_j} \right)^*
	= \sum_{j=1}^n p_j \bra{\psi_j}^* \ket{\psi_j}^*
	= \sum_{j=1}^n p_j \ketbra{\psi_j}
	= \rho.
\end{align*}
Now let $0 \ne \xi \in \msH$.  Then
\begin{align*}
	\braCket \xi \rho \xi
	&= \BraCKet
			{\xi}
			{\sum_{j=1}^n p_j \ketbra{\psi_j}}
			{\xi}
	= \sum_{j=1}^n p_j \braket{\xi}{\psi_j} \braket{\psi_j}{\xi}
	= \sum_{j=1}^n p_j \braket{\xi}{\psi_j} \braket{\xi}{\psi_j}^*
	= \sum_{j=1}^n p_j | \braket{\xi}{\psi_j} |^2.
\end{align*}
Since this is a sum of non-negative real numbers, it is non-negative.  But by
positive-definiteness of the inner product, it is zero iff either $\xi$ is
zero, which it is not by hypothesis, or all of the $\psi_j$'s are zero, which
they are not since they all have norm 1.

To see that $\rho$ has trace 1, write
\begin{align}
	\Tr(\rho) &= \Tr\left(\sum_{k=1}^n p_k \ketbra{\psi_k}\right) \\
	&= \sum_{k=1}^n p_k \Tr\left(\ketbra{\psi_k}\right).
\end{align}
Write each $\psi_k$ in terms of the orthonormal basis as
$$
	\psi_k = \sum_\ell c_{k\ell} \phi_\ell.
$$
Note that $\ketbra{\psi_k}$ has $ij$th matrix element
$$
	\left( \ketbra{\psi_k} \right)_{ij} = c_{ki}\,c_{kj}^*.
$$
Then
$$
	\Tr\left( \ketbra{\psi_k} \right) =
		\sum_i c_{ki}\,c_{ki}^* =
		\sum_i | c_{ki} |^2 = 1,
$$
from which
$$
	\Tr(\rho) = \sum_{k=1}^n p_k = 1.
$$
\end{proof}
\end{prop}

Furthermore, we have a criterion for whether a density matrix corresponds
to a pure or mixed state.

\begin{prop}
If $\rho$ is a density matrix for an ensemble $\{\psi_1, \ldots, \psi_n\}$,
then $\Tr(\rho^2) \le 1$, with equality if and only if the ensemble is a pure
state.
\begin{proof}
First note that for a matrix $A$, we have
$$
	\Tr(A^2) =
		\sum_i
		\sum_j
		A_{ij}
		A_{ji}.
$$
Then
\begin{align*}
	\Tr(\rho^2) &=
		\sum_i
		\sum_j
		\rho_{ij}
		\rho_{ji} \\
	&=
		\sum_i
		\sum_j
		\left(
		\sum_{k=1}^n p_k
			c_{ki}
			c_{kj}^*
		\right)
		\left(
		\sum_{\ell=1}^n p_\ell
			c_{\ell j}
			c_{\ell i}^*
		\right)
		\\
	&=
		\sum_{k=1}^n p_k
		\sum_{\ell=1}^n p_\ell \;
		\sum_i
		\left(
			c_{\ell i}^*
			c_{ki}
		\right)
		\sum_j
		\left(
			c_{kj}^*
			c_{\ell j}
		\right)
		\\
	&=
		\sum_{k=1}^n p_k
		\sum_{\ell=1}^n p_\ell \;
		\braket{\psi_\ell}{\psi_k}
		\braket{\psi_k}{\psi_\ell}
		\\
	&=
		\sum_{k=1}^n p_k
		\sum_{\ell=1}^n p_\ell \;
		\braket{\psi_k}{\psi_\ell}^*
		\braket{\psi_k}{\psi_\ell}
		\\
	&=
		\sum_{k=1}^n p_k
		\sum_{\ell=1}^n p_\ell \;
		| \braket{\psi_k}{\psi_\ell} |^2.
\end{align*}
First suppose that one $p_{k_0}$ is 1 and the rest are zero.  (This is the pure
case.)  Then
\begin{align*}
	\Tr(\rho^2) &=
		| \braket{\psi_{k_0}}{\psi_{k_0}} |^2 = 1.
\end{align*}
Now suppose that at least two $p_k$'s are non-zero.  Decompose the sum as
\begin{align*}
	\Tr(\rho^2) &=
		\sum_{k=1}^n p_k^2
		| \braket{\psi_k}{\psi_k} |^2
	+
		\sum_{k=1}^n p_k
		\sum_{\ell \ne k} p_\ell \;
		| \braket{\psi_k}{\psi_\ell} |^2 \\
	&=
		\sum_{k=1}^n p_k^2
	+
		\sum_{k=1}^n
		\sum_{\ell \ne k}
		p_k
		p_\ell \;
		| \braket{\psi_k}{\psi_\ell} |^2.
\end{align*}
To finish the proof, note that $\sum_{k=1}^n p_k^2$ is the diagonal part of
$(\sum_{k=1}^n p_k)^2$.  We have
\begin{align*}
	1 = (1)^2 &= \left(\sum_{k=1}^n p_k\right)^2
	= \sum_{k=1}^n \sum_{\ell=1}^n
		p_k p_\ell
	= \sum_{k=1}^n \sum_{\ell \ne k}
		p_k p_\ell
		+
		\sum_{k=1}^n
		p_k^2
	\\
	\sum_{k=1}^n p_k^2 &= 1 - \sum_{k=1}^n \sum_{\ell \ne k} p_k p_\ell.
\end{align*}
Then
\begin{align*}
	\Tr(\rho^2) &=
	1 - \sum_{k=1}^n \sum_{\ell \ne k} p_k p_\ell
		(1 - \cos^2\theta_{k\ell})
\end{align*}
where
$$
	\braket{\psi_k}{\psi_\ell} = \|\psi_k\| \;\|\psi_\ell\|
		\;\cos\theta_{k\ell}.
$$
Note that all the terms in the sum are non-negative.  If two $p_k$'s are
non-zero --- without loss of generality, say $p_1$ and $p_2$ --- then the
product $p_1p_2$ is non-zero.  Furthermore, $\cos^2\theta_{1,2}$ is strictly
less than 1.  Otherwise, $\phi_1$ and $\phi_2$ (which are unit vectors)
would be identical up to global phase, contradicting the distinctness of the
members of the ensemble.  Thus, $\Tr(\rho^2) < 1$.
\end{proof}
\end{prop}

\begin{rem}
Two different ensembles can give the same density matrix.  In \cite{bib:nc} is
a theorem characterizing the conditions under which this can happen.
\end{rem}

% [XXX type this up and include an example.]

Lastly, we have two lemmas which will be used in section \ref{sec:pure_dm}.

\begin{lem}
\label{lem:rho_herm_idempotent}
The operator $\rho$ is Hermitian.  If $\rho=\ket\psi$,
then $\rho$ satisfies $\rho^2=\rho$.
\begin{proof}
For the first claim,
$$
	\rho^* =
	\left(\sum_k p_k \ketbra{\psi_k} \right)^* =
	\sum_k p_k \bra{\psi_k}^* \ket{\psi_k}^* =
	\sum_k p_k \ketbra{\psi_k} = \rho.
$$
For the second claim,
\begin{align*}
	\rho^2 = \ket\psi \braket\psi\psi \ket\psi = \ket\psi \bra\psi = \rho
\end{align*}
since $\braket \psi\psi = 1$.
\end{proof}
\end{lem}

\begin{lem}
\label{lem:expn_trace}
For a Hermitian operator $A$,
$$
	\braCket{\psi}{A}{\psi} = \Tr(A \ket\psi \bra\psi).
$$
\begin{proof}
The operator $\rho=\ketbra\psi$ is Hermitian by lemma
\ref{lem:rho_herm_idempotent}, so it makes sense to write $\braCket \psi A
\psi$.  Expand $\psi$ in terms of an orthonormal basis $\{\phi_j\}$ for $\msH$
as
$$
	\psi = \sum_j c_j \phi_j.
$$
Recall that
$$
	\braCket \psi A \psi = \sum_{ik} c_i^* A_{ik} c_k.
$$
On the other hand, for matrices $C$ and $D$,
\begin{align*}
	(CD)_{ij} &= \sum_k C_{ik} D_{kj}
	\qqand
	\Tr(CD)   = \sum_{ik} C_{ik} D_{ki},
\end{align*}
so
\begin{align*}
	\Tr(A \ketbra\psi) &= \sum_{ik} A_{ik} c_k c_i^*
	= \braCket \psi A \psi.
\end{align*}
\end{proof}
\end{lem}

%% ----------------------------------------------------------------
\subsection{Pure-state density matrices}
\label{sec:pure_dm}

For this section, fix a state $\psi$ and let
$$
	\rho = \ketbra\psi.
$$
Here, four key properties of the density matrix $\rho$ are described,
paralleling the four postulates of quantum mechanics listed in section
\ref{sec:postulates}.

\begin{itemize}

\item[(1)]  The state vector $\psi$ has norm 1.  The density matrix $\rho$
has trace 1, as proved in proposition \ref{prop:rho_from_ens}.

\item[(2)]  The state vector $\psi$ evolves in time according to
the \plainidx{Schr\"odinger equation}
$$
	\frac{d\psi}{dt} = \frac{-i}{\hbar} \, H \psi;
$$
the density matrix $\rho$ evolves in time according to
\begin{align}
\label{eqn:pure_state_drho_dt}
	\frac{d\rho}{dt} = \frac{-i}{\hbar} \, [H, \rho].
\end{align}

\begin{proof}
Using the product rule,
\begin{align*}
	\frac{d\rho}{dt}
		&= \frac{d\ket\psi}{dt} \bra\psi + \ket\psi \frac{d\bra\psi}{dt}.
\end{align*}
To find $d\bra\psi/dt$, take the adjoint of the Schr\"odinger equation and
recall that the Hamiltonian is Hermitian (so $H^*=H$):
\begin{align*}
	\left(\frac{d\ket\psi}{dt}\right)^*
		&= \left(\frac{-i}{\hbar} H \ket\psi \right)^* \\
	\frac{d\bra\psi}{dt}
		&= \frac{i}{\hbar} \bra\psi H^*
		= \frac{i}{\hbar} \bra\psi H.
\end{align*}
Then
\begin{align*}
	\frac{d\rho}{dt}
	&= \left(\frac{-i}{\hbar}\,H\ket\psi \right) \bra\psi
		+ \ket\psi \left( \frac{i}{\hbar}\,\bra\psi H \right)
	= \frac{-i}{\hbar}\,H \rho + \frac{i}{\hbar} \rho H
	= \frac{-i}{\hbar}\, [H,\rho].
\end{align*}
\end{proof}

In section \ref{sec:mixed_dm} we will find a solution for this equation.

\item[(3)]  Given an observable $\mcO$ with operator $A$, eigenstates
$\phi_j$, and corresponding eigenvalues $a_j$, we have
$$
	\bbP(\mcO = a_j) = |\braket{\phi_j}{\psi}|^2
	\qqand
	\expn{A} = \braCket{\psi}{A}{\psi};
$$
in terms of $\rho$, we have
$$
	\bbP(\mcO = a_j) = \Tr(\rho P_j)
	\qqand
	\expn{A} = \Tr(\rho A)
$$
where
$$
	P_j = \ketbra{\phi_j}.
$$

\begin{proof}
For the first claim,
\begin{align*}
	\bbP(\mcO = a_j) &= |\braket{\phi_j}{\psi}|^2
	= \braket{\psi}{\phi_j} \braket{\phi_j}{\psi}
	= \braCket{\psi}{P_j}{\psi}
	= \Tr(\rho P_j),
\end{align*}
where the last step follows from the \plainidx{centrality of trace}, along with
lemma \ref{lem:expn_trace}.  The second claim also follows from the lemma.
\end{proof}

\item[(4)]  After \plainidx{measurement} of the \plainidx{observable} $\mcO$
with operator $A$, the state vector $\psi$ is replaced by an eigenstate
$\phi_j$.  The density matrix $\ketbra\psi$ is replaced by
$\ket{\phi_j}\bra{\phi_j}$.

\end{itemize}

%% ----------------------------------------------------------------
\subsection{Mixed-state density matrices}
\label{sec:mixed_dm}

Now we suppose
$$
	\rho = \sum_k p_k \ketbra{\psi_k},
$$
and we enumerate the postulates of quantum mechanics in terms of the
\plainidx{density matrix} $\rho$.

\begin{itemize}

\item[(1)]  The density matrix $\rho$ has trace 1, as proved in proposition
\ref{prop:rho_from_ens}.

\item[(2)]  The density matrix $\rho$ evolves in time according to
\begin{align}
\label{eqn:mixed_state_drho_dt}
	\frac{d\rho}{dt} = \frac{-i}{\hbar} \, [H, \rho].
\end{align}

\begin{proof}
Using the product rule,
\begin{align*}
	\frac{d\rho}{dt}
		&= \sum_k p_k \frac{d\ket{\psi_k}}{dt} \bra{\psi_k}
		+ \sum_k p_k \ket{\psi_k} \frac{d\bra{\psi_k}}{dt}.
\end{align*}
As in section \ref{sec:pure_dm},
$$
	\left(\frac{d\ket{\psi_k}}{dt}\right)^* = \frac{i}{\hbar} \bra{\psi_k} H,
$$
so
\begin{align*}
	\frac{d\rho}{dt}
	&= \sum_k p_k \left(\frac{-i}{\hbar}\,H\ket{\psi_k} \right) \bra{\psi_k}
	+ \sum_k p_k \ket{\psi_k} \left( \frac{i}{\hbar}\,\bra{\psi_k} H \right) \\
%
	&= \frac{-i}{\hbar} \sum_k H p_k \ketbra{\psi_k}
	+ \frac{ i}{\hbar} \sum_k p_k \ketbra{\psi_k} H \\
%
	&= \frac{-i}{\hbar}\,H \rho + \frac{i}{\hbar} \rho H
	= \frac{-i}{\hbar}\, [H,\rho].
\end{align*}
\end{proof}

Furthermore, if \plainidx{wave functions} evolve by a unitary operator $U_t$
via
$$
	\psi_t = U_t \psi_0,
$$
then the density matrix evolves via
$$
	\rho_t = U_t \rho_0 U_t^*.
$$
In particular, if the Hamiltonian is \plainidx{time-independent}, we have
\begin{align}
\label{eqn:rho_evolution}
	\rho_t = e^{-iHt/\hbar} \rho_0 e^{iHt/\hbar}.
\end{align}
One verifies by differentiation that \ref{eqn:rho_evolution} solves
\ref{eqn:mixed_state_drho_dt}.

\item[(3)]  Given an observable $\mcO$ with operator $A$, eigenstates
$\phi_j$, and corresponding eigenvalues $a_j$,
$$
	\bbP(\mcO = a_j) = \Tr(\rho P_j)
	\qqand
	\expn{A} = \Tr(\rho A)
$$
where
$$
	P_j = \ketbra{\phi_j}.
$$

\begin{proof}
For the first claim, let $\mcO_k$ be the observation of $\psi_k$.  Using
\plainidx{conditional probability}, recalling result 3 from section
\ref{sec:pure_dm}, and using linearity of trace, we have
\begin{align*}
	\bbP(\mcO = a_j) &=
		\sum_k p_k \, \bbP(\mcO_k = a_j) \\
	&= \sum_k p_k \, \Tr\left( \ketbra{\psi_k} P_j \right) \\
	&= \Tr\left( \sum_k p_k \, \ketbra{\psi_k} P_j\right) \\
	&= \Tr\left( \left(\sum_k p_k \, \ketbra{\psi_k} \right)
		P_j\right) \\
	&= \Tr(\rho P_j).
\end{align*}
Similarly, for the second claim we have
\begin{align*}
	\expn{A} &= \sum_k p_k \expn{A_k} \\
	&= \sum_k p_k \Tr\left( \ketbra{\psi_k} A \right) \\
	&= \Tr\left( \sum_k p_k \ketbra{\psi_k} A \right) \\
	&= \Tr(\rho A).
\end{align*}
\end{proof}

\item[(4)]  After measurement of the observable $\mcO$ with operator $A$, the
density matrix $\rho$ is replaced by
$$
	\sum_j P_j \rho P_j.
$$
\begin{proof}
First,
\begin{align*}
	\sum_j P_j \rho P_j
	&= \sum_j
		\left(\ketbra{\phi_j}\right)
		\left(\sum_k p_k \ketbra{\psi_k}\right)
		\left(\ketbra{\phi_j}\right) \\
	&= \sum_{jk}
		p_k
		\ket{\phi_j}
		\braket{\phi_j}{\psi_k}
		\braket{\psi}{\phi_j}
		\bra{\phi_j} \\
	&= \sum_{jk}
		p_k
		\ket{\phi_j}
		| \braket{\phi_j}{\psi} |^2
		\bra{\phi_j} \\
	&= \sum_{jk}
		p_k
		|c_{kj}|^2
		\ketbra{\phi_j}.
\end{align*}
On the other hand, from result 4 in section \ref{sec:pure_dm} we know that each
state $\psi_k$ collapses to eigenstate $\phi_j$ with probability
$$
	|c_{kj}|^2 = |\braket{\phi_j}{\phi_k}|^2.
$$
Thus, again using conditional probability, the post-measurement ensemble should
be described by pure-state density matrices $\ketbra{\phi_j}$ with
probabilities $p_k|c_{kj}|^2$, i.e.
$$
	\sum_{jk} p_k |c_{kj}|^2 \ketbra{\phi_j} = \sum_j P_j \rho P_j.
$$
\end{proof}
\end{itemize}

%% ----------------------------------------------------------------
\subsection{Additional properties of density matrices}

\begin{lem}
\label{lem:sum_ketbra_phi_j_is_I}
If $\{\phi_j\}$ is an orthornomal basis for $\msH$, then
$$
	\sum_j \ketbra{\phi_j} = I.
$$
\begin{proof}
Take an arbitrary vector $\psi \in \msH$.  Expanded in terms of the basis,
this is
$$
	\psi = \sum_j c_j \phi_j
$$
where
$$
	c_j = \braket{\phi_j}{\psi}.
$$
Then
\begin{align*}
	\left( \sum_j \ketbra{\phi_j} \right) \ket\psi
	= \sum_j \ket{\phi_j} \braket{\phi_j}{\psi}
	= \sum_j \ket{\phi_j} c_j
	= \psi = I \psi.
\end{align*}
\end{proof}
\end{lem}

\begin{lem}
\label{lem:ebH_eig}
Eigenvectors of a self-adjoint operator $H$ are eigenvectors of $e^{-\beta H}$.
That is, if $H\phi_j = E_j\phi_j$ then
$$
	e^{-\beta H} \phi_j = e^{-\beta E_j} \phi_j.
$$
\begin{proof}
This may be justified formally by using the Taylor expansion for the
operator exponential:
\begin{align*}
	e^{-\beta H} \phi_j
	&= \sum_{k=0}^\infty \frac{(-\beta)^k H^k \phi_j}{k!}
	= \sum_{k=0}^\infty \frac{(-\beta)^k E_j^k \phi_j}{k!}
	= e^{-\beta E_j} \phi_j.
\end{align*}
\end{proof}
\end{lem}

%% ----------------------------------------------------------------
\subsection{Trace in coordinates}

With\index{trace in coordinates} respect to a basis $\{\phi_i\}$, $A$ has
matrix element $A_{ij}$
$$
	A_{ij} = \braCket{\phi_i}{A}{\phi_j}.
$$
Then
\begin{align}
\label{eqn:Tr_bracket}
	\Tr(A) &= \sum_i A_{ii}
	= \sum_i \braCket{\phi_i}{A}{\phi_i}.
\end{align}
If $\{\phi_i\}$ is an orthonormal eigenbasis for $A$, with respective
eigenvalues $\{\lambda_i\}$, then
\begin{align*}
	\Tr(A)
	= \sum_i \braCket{\phi_i}{\lambda_i}{\phi_i}
	= \sum_i \lambda_i \braket{\phi_i}{\phi_i}
	= \sum_i \lambda_i.
\end{align*}
Likewise, using lemma \ref{lem:ebH_eig},
\begin{align*}
	\TrebH
	= \sum_i \braCket{\phi_i}{e^{-\beta H}}{\phi_i}
	= \sum_i \braCket{\phi_i}{e^{-\beta E_i}}{\phi_i}
	= \sum_i e^{-\beta E_i} \braket{\phi_i}{\phi_i}
	= \sum_i e^{-\beta E_i}.
\end{align*}

%% ----------------------------------------------------------------
\subsection{Density matrix for the canonical ensemble}
\label{sec:canon_ens_rho}

Here we work out the density matrix for the \plainidx{canonical ensemble}.  The
$p_j$'s of proposition \ref{prop:rho_from_ens} are
$$
	p_j = \frac{e^{-\beta E_j}}{Z}.
$$
with $Z$ being a to-be-determined \plainidx{normalization factor} that makes
$\sum_j p_j=1$.  The $\phi_j$'s are eigenstates of the Hamiltonian with
respective eigenvalues $E_j$.  Then
\begin{align*}
	\rho
	&= \frac{1}{Z} \sum_j e^{-\beta E_j} \ketbra{\phi_j}
	= \frac{1}{Z} \sum_j e^{-\beta H} \ketbra{\phi_j}
	= \frac{1}{Z} e^{-\beta H} \sum_j \ketbra{\phi_j}
	= \frac{1}{Z} e^{-\beta H}.
\end{align*}
by lemmas
\ref{lem:sum_ketbra_phi_j_is_I}
and
\ref{lem:ebH_eig}.

It is now clear what the \plainidx{partition function} $Z$ is:  since trace is
linear, for a trace-class operator $A$ we have
$$
	\Tr\left(\frac{A}{\Tr(A)}\right) = \frac{\Tr(A)}{\Tr(A)} = 1.
$$
Thus, to normalize $\rho$ (see property 1 in section \ref{sec:mixed_dm}),
$$
	Z = \TrebH
$$
so
$$
	\rho = \frac{e^{-\beta H}}{\TrebH}.
$$
Then property (3) of section \ref{sec:mixed_dm} becomes
$$
	\expn{A} = \frac{\TrXebH{A})}{\TrebH}.
$$

%% ================================================================
\newpage
\section{Thermodynamics}

Here I collect certain facts about thermodynamics, gleaned from
\cite{bib:bu,bib:qmath,bib:lssy,bib:huang}.  In particular, I note here
questions which I would like to resolve to my own satisfaction during my
dissertation work.

%% ----------------------------------------------------------------
\subsection{Chemical potential}
\label{sec:mu}

The \bfidx{chemical potential} $\mu$ of a collection of particles is the
change in energy that results when a particle is added to the system:
$$
	dE = \mu\,dN.
$$

%% ----------------------------------------------------------------
\subsection{Grand-canonical partition function; Fourier space}

The \plainidx{partition function} $N$ in section \ref{sec:model_defns} is
what is known as a \bfidx{canonical partition function}:  it involves a fixed
particle number $N$.  In \plainidx{statistical mechanics}, one also considers a
\bfidx{grand-canonical partition function}, defined as
\begin{align*}
	\Xi(\beta,\Lambda,\mu)
	&= \sum_{N \ge 0} e^{\beta \mu N} Z(N).
\end{align*}

For the two-cycle-interaction model of section \ref{sec:n2}, we have
\begin{align*}
	\Xi(\beta,\Lambda,\mu)
%
	&= \sum_{N \ge 0} \frac{e^{\beta\mu N}}{N!}
		\sum_{\pi\in\SN}
		\int_{\Lambda^N} d\vecX
		\, e^{-H^{(\alpha)}(\vecX,\pi)} \\
%
	&= \sum_{N \ge 0} \frac{e^{\beta\mu N}}{N!}
		\sum_{\pi\in\SN}
		\int_{\Lambda^N} d\vecX
		\, \exp\left\{
			-\frac{1}{4\beta}\sum_{i=1}^N \|\vecx_i-\vecxpii\}^2
		\right\}
		\, \exp\left\{
			-\alpha N_2(\pi)
		\right\} \\
%
	&= \sum_{N \ge 0} \frac{e^{\beta\mu N}}{N!}
		\sum_{\pi\in\SN}
		\, \exp\left\{ -\alpha N_2(\pi) \right\}
		\int_{\Lambda^N} d\vecX
		\, \prod_{i=1}^N \exp\left\{
			-\frac{1}{4\beta} \|\vecx_i-\vecxpii\|^2
		\right\}.
\end{align*}
One may write down a grand-canonical partition function in momentum space
rather than in position space:  one transforms the $\vecx$ coordinates into
$\veck$ coordinates.  For the two-cycle-interaction model of section
\ref{sec:n2}, we have
\begin{align*}
	\Xihat(\beta,\Lambda,\mu)
	&= \sum_{N \ge 0} \frac{e^{\beta\mu N}}{N!}
		\sum_{\pi\in\SN}
		\, \exp\left\{ -\alpha N_2(\pi) \right\}
		\sum_{\veck_1, \ldots, \veck_N \in \Lambda^*}
		\, \prod_{i=1}^N \exp\left\{
			-4 \pi^2 \beta \|\veck_i\|^2
		\right\}
		\delta_{\veck_i, \veckpii}.
\end{align*}

An item to be merged into this discussion is the \bfidx{occupation number}.
Huang, p. 112 \cite{bib:huang}, derives the occupation number for the $\veck$th
state to be
\begin{align*}
	n_\veck
	= \frac{e^{\beta \mu}}{e^{\beta E_\veck} - e^{\beta \mu}}
	= \frac{1}{e^{\beta (E_\veck-\mu)} - 1}.
\end{align*}

%% ----------------------------------------------------------------
\subsection{Free energy}

The \bfidx{free energy} is defined in several ways.

In \cite{bib:lssy}, we have
$$
	F = -\frac{1}{\beta}\log Z.
$$

In \cite{bib:qmath}, the free energy is taken to be the
\plainidx{Legendre transform} of the pressure:
$$
	f^{(\alpha)}(\beta,\rho) = \sup_\mu\left[
		\rho \mu - p^{(\alpha)}(\beta,\mu)
		\right].
$$

%% ----------------------------------------------------------------
\subsection{Pressure}
\label{sec:p}

\cite{bib:huang} (p. 111) obtains the pressure of the Bose gas from first
principles.  He sums the momentum distribution for the Bose gas, with pressure
in terms of particle flux on the wall of a container.
$$
	p = \frac{2}{3} \int_{\R^3} \frac{d^3k}{8\pi^3}
		\frac{E_\veck}{e^{\beta(E_\veck-\mu)} - 1}.
$$

\cite{bib:lssy} define the pressure in terms of a partial derivative
involving the \plainidx{free energy}:
$$
	p = -\frac{\partial F}{\partial V} = \rho k_B T.
$$

\cite{bib:qmath} uses
$$
	p = \lim_{V \to \infty} \frac{1}{\beta V} \log Z(\beta, \Lambda, \mu).
$$

%% ----------------------------------------------------------------
\subsection{Density for the Bose gas}
\label{sec:rho}

\cite{bib:lssy} (p. 3) declare
$$
	\rho = \frac{1}{h^3} \int \frac{d\veck}
		{e^{\beta(H_\veck-\mu)}-1}.
$$

\cite{bib:huang} (p. 111) has
$$
	\rho = \frac{1}{2\pi^2} \int_0^\infty dk \frac{k^2}{e^{\beta(H_k-\mu)}-1}.
$$

\cite{bib:qmath} uses
$$
	\rho_c^{(\alpha)} = \frac{\partial p^{(\alpha)}}{\partial\mu}
		\Bigg|_{\mu=0-}
$$
which yields (equations \ref{eqn:rho_c_0} and \ref{eqn:rho_c_alpha})
\begin{align*}
	\rho_c^{(\alpha)}
	& = \frac{\partial p^{(\alpha)}}{\partial\mu}
		\Bigg|_{\mu=0-}
	= \rho_c^{(0)} - \frac{(1-e^{-\alpha})}{2^{9/2}\pi^{3/2}\beta^{3/2}} \\
%
	\rho_c^{(0)} &= \int_{\R^d} \frac{d\veck}{e^{4\beta\pi^2|\veck|^2}-1}.
\end{align*}

It appears that
$$
	\rho = \frac{\partial p}{\partial\mu}
	\qqand
	\rho_c = \frac{\partial p}{\partial\mu} \Bigg|_{\mu=0-}.
$$

%% ================================================================
\newpage
\section{Analysis results}

The results here are standard --- see, for example, \cite{bib:HN}.  They are
collected here for ready reference, using notation which is consistent with the
rest of the paper.  Most of the items here are elaborations on the sketches
provided in \cite{bib:faris}, or justifications of it-can-be-shown statements
in \cite{bib:qmath}.

Importantly, here I identify some otherwise-surprising factors of 2 which
appear in this paper as well as in \cite{bib:qmath} and \cite{bib:bu}:  the
\plainidx{generator of Brownian motion} is $\nabla^2/2$, but throughout the
present work, $H=-\nabla^2$.

%% ----------------------------------------------------------------
\subsection{Gaussians}
\label{sec:gaussians}

\begin{defn}
\label{defn:g_t}
For\index{$g_t$} $\vecx \in \R^d$, let
$$
	g_t(\vecx) = \frac{1}{(2\pi t)^{d/2}}\, e^{-|\vecx|^2/2t}.
$$
Then $g_t(\vecx-\veca)$ is a standard \plainidx{Gaussian}, with integral 1,
having \plainidx{mean} $\veca$ and \plainidx{variance} $t$ in each component.
\end{defn}

\begin{prop}
The family of $g_t$'s satisfies a \plainidx{Chapman-Kolmogorov equation}:
\label{prop:g_iter_2}
$$
	\int_{\R^d} g_s(\vecx-\vecz)g_t(\vecz-\vecy)\,d\vecz = g_{s+t}(\vecx-\vecy).
$$
\begin{proof}
This follows (albeit messily) by completing the square.
\end{proof}
\end{prop}

The following proposition is used for the definition of normalized Brownian
bridges in section \ref{sec:normalized_bridges}, as well as for the simplified
jump-pair interaction in section \ref{sec:2w_to_1w}.
\begin{prop}
\label{prop:g_iter}
\begin{align*}
	\intRdntimes{d}{n-1}
	&
	g_{t_1}(\vecx-\vecz_1)
	g_{t_2}(\vecz_1-\vecz_2)
	\cdots
	g_{t_{n-1}}(\vecz_{n-2}-\vecz_{n-1})
	g_{t_{n}}  (\vecz_{n-1}-\vecy)
	\\
%
	&
	\,d\vecz_1
	\,d\vecz_2
	\cdots
	\,d\vecz_{n-2}
	\,d\vecz_{n-1} \\
%
	&
	= g_{t_1 + \cdots + t_n}(\vecx-\vecy).
\end{align*}
\begin{proof}
Use induction, iterating the integrals and applying proposition
\ref{prop:g_iter_2}.
\end{proof}
\end{prop}

The following lemma is used in the proof of proposition
\ref{prop:2w_to_1w} in section \ref{sec:2w_to_1w}.
\begin{lem}
\label{lem:g2tg2t_g4t_gt_ab}
\begin{align*}
	g_{2t}(\veca_1-\veca_2)
	\;
	g_{2t}(\vecb_1-\vecb_2)
%
	&=
	g_{4t}\left((\veca_1-\vecb_1) - (\veca_2-\vecb_2)\right)
	\;
	g_{ t}\left(\veca_1-\veca_2
		- \frac{1}{2}\left(\veca_1-\vecb_1\right)
		+ \frac{1}{2}\left(\veca_2-\vecb_2\right) \right) \\
%
	&=
	g_{4t}\left((\veca_1-\vecb_1) - (\veca_2-\vecb_2)\right)
	\;
	g_{ t}\left(\frac{(\veca_1-\veca_2) + (\vecb_1-\vecb_2)}{2}\right).
\end{align*}
\begin{proof}
From definition \ref{defn:g_t}, the left-hand side is
\begin{align*}
	\frac{1}{(4\pi t)^{d}}\,
	\exp\left\{-\frac{1}{4t} \left(
		\|\veca_1-\veca_2\|^2 +
		\|\vecb_1-\vecb_2\|^2
	\right) \right\}.
\end{align*}
The right-hand side is
\begin{align*}
	&
	\frac{1}{(8\pi t)^{d/2}}\,
	\frac{1}{(2\pi t)^{d/2}}\,
	\exp\left\{-\frac{1}{8t}
		\|(\veca_1-\vecb_1) - (\veca_2-\vecb_2)\|^2
	\right\}
	\exp\left\{-\frac{1}{2t}
		\left\|
		\veca_1-\veca_2
		- \frac{ \left(\veca_1-\vecb_1\right) }{2}
		+ \frac{ \left(\veca_2-\vecb_2\right) }{2}
		\right\|^2
	\right\} \\
%
	&=
	\frac{1}{(4\pi t)^{d}}\,
	\exp\left\{-\frac{1}{8t}
		\|(\veca_1-\vecb_1) - (\veca_2-\vecb_2)\|^2
	\right\}
	\exp\left\{-\frac{1}{2t}
		\left\|
		\frac{\veca_1}{2}
		-\frac{\veca_2}{2}
		+\frac{\vecb_1}{2}
		-\frac{\vecb_2}{2}
		\right\|^2
	\right\} \\
%
	&=
	\frac{1}{(4\pi t)^{d}}\,
	\exp\left\{-\frac{1}{8t}
		\|(\veca_1-\veca_2) - (\vecb_1-\vecb_2)\|^2
	\right\}
	\exp\left\{-\frac{1}{8t}
		\|(\veca_1-\veca_2) + (\vecb_1-\vecb_2)\|^2
	\right\}.
%
\end{align*}
Since the scale factors are the same, it remains to show that
\begin{align*}
%
	2 \|\veca_1-\veca_2\|^2 + 2 \|\vecb_1-\vecb_2\|^2
%
	&=
%
	\|(\veca_1-\veca_2) - (\vecb_1-\vecb_2)\|^2
	+
	\|(\veca_1-\veca_2) + (\vecb_1-\vecb_2)\|^2.
%
\end{align*}
To see this, let
$$
	\vecu = \veca_1-\veca_2
	\qqand
	\vecv = \vecb_1-\vecb_2.
$$
Then we only need to show
\begin{align*}
	2 \|\vecu\|^2 + 2 \|\vecv\|^2
	&=
	\|\vecu - \vecv\|^2 + \|\vecu + \vecv\|^2.
\end{align*}
But this follows immediately since
\begin{align*}
	\|\vecu \pm \vecv\|^2 &= \|\vecu\|^2 \pm 2 \vecu\cdot\vecv + \|\vecv\|^2.
\end{align*}
\end{proof}
\end{lem}

The following corollary is also used in the proof of proposition
\ref{prop:2w_to_1w} in section \ref{sec:2w_to_1w}.  It is nothing more than a
relabeling of lemma \ref{lem:g2tg2t_g4t_gt_ab}, with
$$
	\veca_1 \mapsto \vecx_1, \quad
	\veca_2 \mapsto \vecy_1, \quad
	\vecb_1 \mapsto \vecx_2, \qand
	\vecb_2 \mapsto \vecy_2.
$$
\begin{cor}
\label{cor:g2tg2t_g4t_gt_xy}
\begin{align*}
	g_{2t}(\vecx_1-\vecy_1)
	\;
	g_{2t}(\vecx_2-\vecy_2)
%
	&=
	g_{4t}\left((\vecx_1-\vecx_2) - (\vecy_1-\vecy_2)\right)
	\;
	g_{ t}\left(\frac{(\vecx_1-\vecy_1) + (\vecx_2-\vecy_2)}{2}\right) \\
%
	&=
	g_{4t}\left((\vecx_1-\vecx_2) - (\vecy_1-\vecy_2)\right)
	\;
	g_{ t}\left(\frac{(\vecx_1+\vecx_2) - (\vecy_1+\vecy_2)}{2}\right).
\end{align*}
\end{cor}


%% ----------------------------------------------------------------
\subsection{$e^{-\beta H_0}$ as a convolution operator}
\label{sec:ebH0}

We seek to write $e^{-\beta H_0}$ as a \plainidx{convolution operator}, where
$H_0$ is a scaled \plainidx{Laplacian}.  Faris includes the arbitrary scale
factor $\sigma^2$ in
$
	H_0 = - \frac{\sigma^2}{2} \nabla^2;
$
we simply need to choose $\sigma^2=2$ to fit the present application,
then adapt the resulting scale factors to the standard Gaussian notation
found in definition \ref{defn:g_t}.

Additionally, we will use $\sigma^2/2=\beta/n$ in sections \ref{sec:fk1_exp}
and \ref{sec:fkboson_exp}.

In the following table, $\mcF$ denotes the \plainidx{Fourier transform};
$F_\beta$ and $\Fhat_\beta$ denote a Gaussian and its Fourier transform as
defined by Faris; $g$ and $\ghat$ are as in definition \ref{defn:g_t}.

$$
\begin{tabular}{|ccc|ccc|ccc|}
	\multicolumn{3}{c}{Arbitrary $\sigma^2$} &
	\multicolumn{3}{c}{$\sigma^2=1$} &
	\multicolumn{3}{c}{$\sigma^2=2$} \\
\hline
	&&&&&&&&\\
	$F_\beta(\vecx)$ & $=$ &
		$\displaystyle{
			\frac{e^{-|\vecx|^2/2\sigma^2 \beta}}{(2\pi\sigma^2 \beta)^{d/2}}
		}$ &
	$F_\beta(\vecx)$ & $=$ &
		$\displaystyle{
			\frac{e^{-|\vecx|^2/2 \beta}}{(2\pi \beta)^{d/2}}
		}$ &
	$F_\beta(\vecx)$ & $=$ &
		$\displaystyle{
			\frac{e^{-|\vecx|^2/4 \beta}}{(4\pi \beta)^{d/2}}
		}$ \\
	&&&&&&&&\\
%
	& $=$ & $g_{\sigma^2\beta}(\vecx)$ &
	& $=$ & $g_{\beta}(\vecx)$ &
	& $=$ & $g_{2\beta}(\vecx)$ \\
	&&&&&&&&\\
%
	$\Fhat_\beta(\veck)$ & $=$ &
		$\displaystyle{e^{-\sigma^2 \beta |\veck|^2/2}}$ &
	$\Fhat_\beta(\veck)$ & $=$ &
		$\displaystyle{e^{- \beta |\veck|^2/2}}$ &
	$\Fhat_\beta(\veck)$ & $=$ &
		$\displaystyle{e^{-\beta |\veck|^2}}$ \\
	&&&&&&&&\\
%
	$H_0$ &$=$& $-\frac{\sigma^2}{2} \nabla^2$ &
	$H_0$ &$=$& $-\frac{1}{2} \nabla^2$ &
	$H_0$ &$=$& $-\nabla^2$ \\
	&&&&&&&&\\
%
	$\mcF\left(-\frac{\sigma^2}{2} \nabla^2\right)$ & $=$ &
		$\frac{\sigma^2}{2}|\veck|^2$ &
	$\mcF\left( -\frac{1}{2} \nabla^2 \right)$ & $=$ &
		$\frac{1}{2}|\veck|^2$ &
	$\mcF\left( - \nabla^2 \right)$ & $=$ &
		$|\veck|^2$ \\
	&&&&&&&&\\
%
	$\mcF\left( e^{-\beta H_0} \right)$ & $=$ & $F_\beta(\veck)$ &
	$\mcF\left( e^{-\beta H_0} \right)$ & $=$ & $F_\beta(\veck)$ &
	$\mcF\left( e^{-\beta H_0} \right)$ & $=$ & $F_\beta(\veck)$ \\
	&&&&&&&&\\
%
	$(e^{-\beta H_0}f)(\vecx)$ & $=$ &
		$(e^{\frac{\sigma^2\beta}{2} \nabla^2}f)(\vecx)$ &
	$(e^{-\beta H_0}f)(\vecx)$ & $=$ &
		$(e^{\frac{\beta}{2} \nabla^2}f)(\vecx)$ &
	$(e^{-\beta H_0}f)(\vecx)$ & $=$ &
		$(e^{\beta \nabla^2}f)(\vecx)$ \\
	&&&&&&&&\\
%
	& $=$ & $F_\beta(\vecx) * f(\vecx)$ &
	& $=$ & $F_\beta(\vecx) * f(\vecx)$ &
	& $=$ & $F_\beta(\vecx) * f(\vecx)$ \\
	&&&&&&&&\\
%
	& $=$ & $g_{\sigma^2\beta}(\vecx) * f(\vecx)$ &
	& $=$ & $g_{\beta}(\vecx) * f(\vecx)$ &
	& $=$ & $g_{2\beta}(\vecx) * f(\vecx)$ \\
	&&&&&&&&\\
\hline
\end{tabular}
$$

We have used the $*$ notation for the \bfidx{convolution}:
$$
	g(\vecx) * f(\vecx) =
	(g*f)(\vecx) = \int_{\R^d} g(\vecx-\vecy) f(\vecy)\,d\vecy.
$$

We now prove the claims tabulated above.  The proof is split up into
a sequence of lemmas.

\begin{defn}
Following Faris, we define the \bfidx{Fourier transform} and
\bfidx{inverse Fourier transform} using the \plainidx{non-unitary}
angular-frequency\index{angular frequency} conventions for sign and scale.
Namely, for $f,g:\R^d\to\R$ which are in $L^1 \cap L^2$,
\begin{align}
	(\mcF f)(\veck) &= \fhat(\veck)
	= \int_{\R^d} f(\vecx) e^{-i\veck\cdot\vecx} \,d\vecx. \\
%
	(\mcF^{-1} g)(\vecx)
	&= \frac{1}{(2\pi)^d}
	\int_{\R^d} g(\veck) e^{+i\veck\cdot\vecx} \,d\veck.
\end{align}
\end{defn}

\begin{lem}
The Fourier transform\index{$\hat{g_t}$} of the one-dimensional $t$-variance
Gaussian (definition \ref{defn:g_t}) is
\label{lem:ghat1d}
\begin{align*}
	\ghat_t(k) &= e^{-tk^2/2}.
\end{align*}
\begin{proof}
Complete the square, use \plainidx{translation invariance} of the integral, and
recognize $1/(2\pi)^{1/2}$ as the integral of $e^{-x^2/2t}$:
\begin{align*}
	(\mcF g_t)
	&= \frac{1}{(2\pi t)^{1/2}} \int_\R e^{-x^2/2t} e^{-ikx}\,dx \\
	&= \frac{1}{(2\pi t)^{1/2}} \int_\R
		\exp\left\{-\frac{1}{2t} (x^2 + 2itkx)\right\}\,dx \\
	&= \frac{1}{(2\pi t)^{1/2}} \int_\R
		\exp\left\{-\frac{1}{2t} (x + itk)^2 - \frac{tk^2}{2} \right\}\,dx \\
	&= \frac{e^{-tk^2/2}}{(2\pi t)^{1/2}} \int_\R
		\exp\left\{-\frac{1}{2t} (x + itk)^2 \right\}\,dx \\
	&= \frac{e^{-tk^2/2}}{(2\pi t)^{1/2}} \int_\R
		\exp\left\{-\frac{1}{2t} x^2 \right\}\,dx \\
	&= e^{-tk^2/2}.
\end{align*}
\end{proof}
\end{lem}

\begin{prop}
The Fourier transform of the $d$-dimensional $t$-variance Gaussian (definition
\ref{defn:g_t}) is
\begin{align*}
	\ghat_t(\veck) &= e^{-t\|\veck\|^2/2}.
\end{align*}
\begin{proof}
First write
\begin{align*}
	(\mcF g_t)
	&= \frac{1}{(2\pi t)^{d/2}} \int_{\R^d}
		e^{-\|\vecx\|^2/2t} e^{-i\veck\cdot\vecx}\,d\vecx \\
	&= \left(\frac{1}{2\pi t}\right)^{d/2} \int_{\R^d}
		\left(\prod_{\ell=1}^d e^{-x_j^2/2t} e^{-ik_jx_j}\right)
		\,d\vecx.
\end{align*}
Since the integrand is non-negative, by Tonelli's theorem we may iterate the
integrals and apply lemma \ref{lem:ghat1d} to each.
\end{proof}
\end{prop}

\begin{lem}
Let $f:\R\to\R$ be such that $f(x),f'(x)$ and $f''(x)$ are in $L^1 \cap L^2$
and vanish at $\pm \infty$.  Then
$$
	\mcF\left(\frac{d^2f}{dx^2}\right) = -k^2 \fhat(k).
$$
\begin{proof}
Integrate by parts twice, applying the zero boundary conditions:
\begin{align*}
	\mcF\left(\frac{d^2f}{dx^2}\right)
	&= \int_\R \frac{d^2f}{dx^2} e^{-ikx}\,dx 
	= ik \int_\R \frac{df}{dx} e^{-ikx}\,dx
	= (ik)^2 \int_\R f(x) e^{-ikx}\,dx
	= -k^2 \fhat(k).
\end{align*}
\end{proof}
\end{lem}

\begin{lem}
\label{lem:FT_of_Lapl}
Let $f:\R^d\to\R$ be such that $f$ and all its first and second partial
derivatives are in $L^1 \cap L^2$ and vanish at $\pm \infty$.  Then
$$
	\mcF\left(\nabla^2 f\right) = -\|\veck\|^2 \fhat(\veck).
$$
\begin{proof}
Integrate by parts twice, applying the zero boundary conditions:
\begin{align*}
	\mcF\left(\nabla^2 f\right)
	&= \int_{\R^d} \nabla^2 f e^{-i\veck\cdot\vecx}\,d\vecx \\
	&= \int_{\R^d} \left(\sum_{j=1}^d \frac{d^2f}{dx_j^2}\right)
		\left(\prod_{\ell=1}^d e^{-ik_\ell x_\ell}\right) \,d\vecx \\
	&= \sum_{j=1}^d \int_{\R^d} \frac{d^2f}{dx_j^2}
		\left(\prod_{\ell=1}^d e^{-ik_\ell x_\ell}\right) \,d\vecx \\
	&= \sum_{j=1}^d (ik_j) \int_{\R^d} \frac{df}{dx_j}
		\left(\prod_{\ell=1}^d e^{-ik_\ell x_\ell}\right) \,d\vecx \\
	&= \sum_{j=1}^d (ik_j)^2 \int_{\R^d} f(\vecx)
		\left(\prod_{\ell=1}^d e^{-ik_\ell x_\ell}\right) \,d\vecx \\
	&= -\|\veck\|^2 \fhat(\veck).
\end{align*}
\end{proof}
\end{lem}

\begin{defn}
\label{def:operator_transform}
The \bfidx{Fourier transform} or \bfidx{operator transform} of an operator $A$
is
$$
	\Ahat = \mcF A \mcF^{-1}.
$$
Note this this is equivalent to
$$
	\Ahat \mcF = \mcF A.
$$
\end{defn}

That is, this is the frequency-domain operator which takes a $\veck$ function,
back-transforms to an $\vecx$ function, applies $A$, and forward-transforms:

\setlength{\unitlength}{1.0in}
\begin{picture}(6.5,1.2)
%\put(0.0,0.0){\framebox(6.5,1.2)}

\put(1.00,0.50){Spatial domain}
\put(5.00,0.50){Frequency domain}
\put(2.70,0.50){$A$}
\put(4.20,0.50){$\Ahat$}
\put(3.40,1.05){$\mcF^{-1}$}
\put(3.40,0.20){$\mcF$}

\put(3.00,0.90){\vector(0,-1){.75}}
\put(4.00,0.90){\vector(0,-1){.75}}

\put(3.90,0.95){\vector(-1,0){.80}}
\put(3.10,0.10){\vector(1,0){.80}}

\end{picture}

\begin{rem}
As a consequence,
\begin{align*}
	\Ahat &= \mcF A \mcF^{-1} \\
	\Ahat \mcF &= \mcF A \\
	\Ahat \mcF f &= (\mcF A) f \\
	\mcF (A f) &= \Ahat \fhat.
\end{align*}
\end{rem}

\begin{rem}
\label{rem:FT_scaled_Lapl}
We just saw in lemma \ref{lem:FT_of_Lapl}, with $A=\nabla^2$, that
\begin{align*}
	\mcF(\nabla^2 f) &= -\|\veck\|^2 \fhat.
\end{align*}
Therefore we conclude that the\index{Laplacian, Fourier transform of}
\plainidx{Fourier transform} of the scaled Laplacian is the multiplication
operator
\begin{align*}
	\mcF(c\nabla^2) &= -c\|\veck\|^2.
\end{align*}
\end{rem}

\begin{prop}
\label{prop:FT_exp_Lapl_1D}
\begin{align*}
	\mcF\left(e^{c\, d^2/dx^2}\right) &= e^{-ck^2}.
\end{align*}
\begin{proof}
Moving the transform through the infinite sum without justification, we have
\begin{align*}
	\mcF\left(e^{c\,d^2f/dx^2}\right)
	&= \mcF \left(\sum_{j=0}^\infty
		\frac{c^j}{j!} \frac{d^{2j}f}{dx^{2j}} \right) \\
	&= \sum_{j=0}^\infty
		\frac{c^j}{j!} \mcF\left( \frac{d^{2j}f}{dx^{2j}} \right) \\
	&= \sum_{j=0}^\infty \frac{c^j}{j!} (-k^{2j})\fhat(k) \\
	&= \left( \sum_{j=0}^\infty \frac{-(ck^2)^j}{j!} \right) \fhat(k) \\
	&= e^{-ck^2} \fhat(k).
\end{align*}
\end{proof}
\end{prop}

\begin{prop}
\label{prop:FT_exp_Lapl}
\begin{align*}
	\mcF\left(e^{c\nabla^2}\right) &= e^{-c\|\veck\|^2}.
\end{align*}
\begin{proof}
TBD.  The procedure is similar to \ref{prop:FT_exp_Lapl_1D}; there are some
combinatorial factors involving mixed partials which I need to work out on
paper.
\end{proof}
\end{prop}

\begin{prop}
The \plainidx{Fourier transform} of a \plainidx{convolution} is the product of
Fourier transforms:
\begin{align*}
	\mcF(f * g) &= \fhat \ghat.
\end{align*}
\begin{proof}
Start with
\begin{align*}
	\mcF\left( f*g \right)
	&= \mcF\left( \int_{\R^d}
		f(\vecy)g(\vecx-\vecy) \,d\vecy \right) \\
	&= \int_{\R^d} \int_{\R^d}
		f(\vecy)g(\vecx-\vecy) e^{-i\veck\cdot\vecx}d\vecy \,d\vecx.
\end{align*}
Now let $\vecz=\vecx-\vecy$.  Then, using Fubini's theorem since imaginary
exponentials have modulus one and we assume $f,g$ are $L^2$ so their product
is $L^1$ (their inner product exists),
\begin{align*}
	\mcF\left( f*g \right)
	&= \int_{\R^d} \int_{\R^d}
		f(\vecy)g(\vecz) e^{-i\veck\cdot(\vecy+\vecz)}d\vecy \,d\vecz \\
	&= \int_{\R^d} \int_{\R^d}
		f(\vecy) e^{-i\veck\cdot\vecy}
		g(\vecz) e^{-i\veck\cdot\vecz} \,d\vecy \,d\vecz \\
	&= \int_{\R^d} \left[ \int_{\R^d}
		f(\vecy) e^{-i\veck\cdot\vecy}
		d\vecy
		\right]
		g(\vecz) e^{-i\veck\cdot\vecz}
		d\vecz \\
	&= \fhat \ghat.
\end{align*}
\end{proof}
\end{prop}

\begin{prop}
The exponentiated scaled Laplacian may be written as an integral operator:
\begin{align*}
	\left(e^{c\nabla^2} f\right)(\vecx)
	&= \int_{\R^d} g_c(\vecx-\vecy) f(\vecy) \,d\vecy.
\end{align*}
\begin{proof}
Let $A=e^{-c\nabla^2}$.  Then we need to show
\begin{align*}
	Af &= g_c*f \\
	\mcF(Af) &= \mcF(g_c*f) \\
	\Ahat \fhat &= \ghat_c \fhat \\
	\Ahat &= \ghat_c.
\end{align*}
But this follows from proposition \ref{prop:FT_exp_Lapl}.
\end{proof}
\end{prop}

%% ----------------------------------------------------------------
\subsection{Operator trace}
\label{sec:operator_trace}

As\index{trace} described in section \ref{sec:canon_ens_rho}, we will want to
compute $\TrebH$.  Here is a partial result, assuming that $e^{-\beta H}$ can
be written as an integral operator.  The construction of that integral operator
will need to wait until section \ref{sec:fk1_ker}.

\begin{prop}
\label{prop:operator_trace}
If a trace-class operator $A$ on a separable Hilbert space has a $G(\vecx,
\vecy)$ such that
$$
	A\, f(\vecx) = \int G(\vecx,\vecy) f(\vecy) \,d\vecy,
$$
then
$$
	\Tr(A) = \int G(\vecx, \vecx) \,d\vecx.
$$
\begin{proof}
As shown in equation \ref{eqn:Tr_bracket}, with $\{\phi_j\}$ being a
(countable) basis for the (separable) Hilbert space $\msH$,
$$
	\Tr(A) = \sum_j \braCket{\phi_j}{A}{\phi_j}.
$$
Then
\begin{align*}
	\Tr(A)
	&= \sum_j \braCket{\phi_j}{A}{\phi_j} \\
	%&= \sum_j \braket{\phi_j} {\int G(\vecx, \vecy) \,\phi_j(\vecy)\,d\vecy} \\
	&= \sum_j \int \int
		\phi_j^*(\vecx) \,G(\vecx, \vecy) \,\phi_j(\vecy)\,d\vecy\,d\vecx \\
	&= \int \int \,G(\vecx, \vecy)
		\left(\sum_j \phi_j^*(\vecx) \,\phi_j(\vecy)\right) \,d\vecy\,d\vecx \\
	&= \int \int \,G(\vecx, \vecy)
		\delta(\vecx-\vecy) \,d\vecy\,d\vecx \\
	&= \int \,G(\vecx, \vecx) \,d\vecx.
\end{align*}

The appearance of the delta function is due to the following
\plainidx{completeness relation}.
\begin{lem}
For all $\vecx,\vecy$,
$$
	\left(\sum_j \phi_j^*(\vecx) \,\phi_j(\vecy)\right) = \delta(\vecx-\vecy).
$$
\begin{proof}
Let $\psi$ be an arbitrary wave function.  Then it is a linear combination
of the basis functions:
$
	\psi(\vecy) = \sum_j c_j \phi_j(\vecy)
$
where
$
	c_j = \int \phi_j^*(\vecx) \psi(\vecx) \,d\vecx.
$
Combining these two statements, we have
\begin{align*}
	\psi(\vecy) &= \sum_j \left(
		\int \phi_j^*(\vecx)\psi(\vecx)\,d\vecx \right) \phi_j(\vecy) \\
	&= \int \left(\sum_j \phi_j^*(\vecx)\phi_j(\vecy)\right)
		\psi(\vecx)\,d\vecx \\
	&= \int \delta(\vecx-\vecy) \psi(\vecx) d\vecx.
\end{align*}
\end{proof}
\end{lem}
\end{proof}
\end{prop}

\begin{rem}
Note that $G(\vecx,\vecy)$ only appears in an integral, and thus may be a
distribution rather than a function.  This will in fact be the case in section
\ref{sec:fk1}, where the $G(\vecx, \vecy)$ we construct will contain a
delta-function term.
\end{rem}

%%% ================================================================
%\newpage
%\section{Thermodynamics}
%
%intensive vs. extensive
%
%$E$
%$U$
%$A$
%
%$T$
%$V$
%$N$
%
%$S$
%$P$
%$\mu$
%
%pairs ...
%
%3 laws ...
%
%\begin{itemize}
%
%\item[(1)]
%\begin{align*}
%	dU &= T\,dS - P\,dV + \mu\,dN \\
%	dA &= -S\,dT - P\,dV + \mu\,dN \\
%	U  &= A - TS.
%\end{align*}
%
%Interp:
%	\begin{itemize}
%
%	\item Entropy describes how the free energy changes as the temperature
%	changes.  The product $T\,dS$ is \bfidx{heat} added to the system.
%
%	\item Pressure describes how the internal energy changes as volume changes.
%	The product $P\,dV$ is \bfidx{work} done on the system.
%
%	\item The chemical potential descries how the internal energy changes as
%	particles are added to the system.  The product $\mu\,dN$ is called [XXX?].
%
%	\end{itemize}
%
%\item[(2)]
%$$
%	\frac{\D S}{\D T} \ge 0
%$$
%and describe when.
%
%\item[(3)]
%$$
%	S = k \ln \Omega
%$$
%so $S$ is minimum when $\Omega=1$ which (XXX xref) happens at $T=0$.
%
%XXX This definition of $S$ is microscopic.  For this section, we need
%a macroscopic version \ldots .
%
%% How to reconcile this with degenerate ground state in the 1D Ising model?
%
%\end{itemize}

%% ----------------------------------------------------------------
\subsection{Triple product of partial derivatives}
\label{sec:triple_partial_product}

This is used for the proof of proposition \ref{prop:delta_Tc} in section
\ref{sec:delta_Tc}.

\begin{prop}
\label{prop:triple_partial_product}
Let $f:\R^3\to\R$ be continuously differentiable.  Let $(x_0,y_0,z_0)$ be a
point on the surface
$$
	f(x,y,z)=0
$$
where
$\Df/\Dx$, $\Df/\Dy$, and $\Df/\Dz$ are non-zero.  Then there is a neighborhood
of $(x_0,y_0,z_0)$ such that
$$
	\fracDxDy \fracDyDz \fracDzDx = -1.
$$
\begin{proof}
Since $\slashDfDx \ne 0$, by the implicit function theorem we can solve for $x$
and write
$$
	f(x(y,z),y,z) = 0.
$$
Differentiating with respect to $y$, we have
\begin{align*}
	\fracDfDx\fracDxDy + \fracDfDy &= 0 \\
	\fracDxDy &= - \frac{ \slashDfDy }{ \slashDfDx}.
\end{align*}
Likewise,
$$
	\fracDyDz = - \frac{ \slashDfDz }{ \slashDfDy}
	\qqand
	\fracDzDx = - \frac{ \slashDfDx }{ \slashDfDz}.
$$
Multiplying the three partials together, we obtain
$$
	\fracDxDy \fracDyDz \fracDzDx 
	= -
	\left(\frac{ \slashDfDy }{ \slashDfDx}\right)
	\left(\frac{ \slashDfDz }{ \slashDfDy}\right)
	\left(\frac{ \slashDfDx }{ \slashDfDz}\right)
	= -1.
$$
\end{proof}
\end{prop}

%%%% ================================================================
%%\newpage
%%\section{Statistical mechanics}
%%
%%XXX mention mechanics of $N$ classical or quantum-mechanical particles, when
%%$N$ can be quite large (cf. $\sim10^{23}$ for a typical laboratory experiment).
%%Average out the microscopic details and obtain simpler macroscopic quantities.
%%
%%N V E T and entropy ...
%%
%%CAREFUL (!) definitions of system, phase space, state, ensemble, distribution
%%
%%%% ----------------------------------------------------------------
%%\subsection{Microcanonical ensemble}
%%
%%$N V E$
%%
%%fixed-E subspaces as foliations (compact submanifolds) of phase space?
%%unif prob distr on cpct?
%%
%%key point is uniform probability.
%%
%%See section \ref{sec:1D_ising} for an example.
%%
%%%% ----------------------------------------------------------------
%%\subsection{Canonical ensemble}
%%
%%$N V T$
%%
%%XXX microstates and macrostates
%%
%%The \bfidx{Boltzmann distribution law}:  Probability of occurrence of
%%microstate with energy $E_i$ is proportional to
%%$$
%%	e^{-\beta E_i}.
%%$$
%%Probability of occurrence of macrostate with energy $E_i$ is proportional to
%%$$
%%	g(E_i) e^{-\beta E_i}
%%$$
%%where $g(E_i)$ is the number of microstates with common energy $E_i$.  One also
%%calls $g(E_i)$ the \emphidx{degeneracy} of the macrostate with energy $E_i$.
%%
%%See section \ref{sec:1D_ising} for an example.
%%
%%XXX motivate in two ways:  (1) Widom's independence argument; (2) Huang?
%%Pathria?  Find it and type it up.
%%
%%The \bfidx{partition function} or \bfidx{canonical partition function}
%%is obtained by \ldots
%%
%%\begin{align*}
%%	\bbP(E = E_i) &= \frac{e^{-\beta E_i}}{\sum_j e^{-\beta E_j}}
%%\end{align*}
%%so
%%\begin{align*}
%%	Z &= \sum_j e^{-\beta E_j}.
%%\end{align*}
%%
%%XXX $Z$ function of $T$, $V$, and $N$.
%%
%%XXX partials ... \bfidx{potentials} ... \cite{bib:widom}:
%%
%%Using conditional probability:
%%\begin{align*}
%%	U = \overline{E} &=
%%		\frac{\sum_i E_i e^{-\beta E_i}}{\sum_j E_j e^{-\beta E_j}}.
%%\end{align*}
%%Using the chain rule one immediately obtains
%%\begin{align*}
%%	U &= -\left( \frac{\D}{\D\beta} \ln Z \right).
%%\end{align*}
%%XXX:
%%\begin{align*}
%%	S &= k \ln \Omega \\
%%	S &= - \sum_i P_i \ln(P_i)
%%\end{align*}
%%XXX:
%%\begin{align*}
%%	U &= A + TS \\
%%	A &= -\ln Z/\beta \\
%%	S &= -\D A / \D T \\
%%	P &= -\D A / \D V \\
%%	\mu &= \D A / \D N.
%%\end{align*}
%%
%%XXX degeneracy of $E$ and $g(E)$ \ldots.
%%
%%%% ----------------------------------------------------------------
%%\subsection{1D Ising-model example}
%%\label{sec:1D_ising}
%%
%%XXX aligning nearest-neighbor energy, periodic boundary conditions:
%%$$
%%	E = -\sum_{i=1}^N \sigma_i \sigma_{i+1}
%%$$
%%with
%%$$
%%	\sigma_{N+1} \equiv \sigma_{1}.
%%$$
%%
%%XXX $N=10$, $\beta=1$:
%%
%%% E=10
%%% -+-+-+-+-+
%%% +-+-+-+-+-
%%%
%%% E=6
%%%
%%% ---+-+-+-+
%%% --+--+-+-+
%%% +++-+-+-+-
%%%
%%% E=2
%%%
%%% ----+-+--+
%%% -+-----+-+
%%% +++-+---+-
%%%
%%% E=-2
%%%
%%% -----++-+-
%%% +++--+++--
%%% ++++++-++-
%%%
%%% E=-6
%%%
%%% -----++---
%%% -+--------
%%% +++--+++++
%%%
%%% E=-10
%%%
%%% ----------
%%% ++++++++++
%%
%%$$
%%\begin{tabular}{|r|r|r|r|r|r|r|}
%%\hline
%%	\multicolumn{1}{|c|}{$i$} &
%%	\multicolumn{1}{|c|}{Representative spins} &
%%	\multicolumn{1}{|c|}{$g(E_i)$} &
%%	\multicolumn{1}{|c|}{$E_i$} &
%%	\multicolumn{1}{|c|}{$e^{-E_i}$} &
%%	\multicolumn{1}{|c|}{$\bbP(E=E_i)$}
%%	\\
%%\hline
%%
%%\hline 1 & $+-+-+-+-+-$ &   2 & $ 10$ &     0.0000454 &  $1.1\cdot 10^{-9}$ \\
%%\hline 2 & $+++-+-+-+-$ &  90 & $  6$ &     0.0024788 &  $2.7\cdot 10^{-6}$ \\
%%\hline 3 & $+++++-+-+-$ & 420 & $  2$ &     0.1353353 &  $6.8\cdot 10^{-4}$ \\
%%\hline 4 & $+++++++-+-$ & 420 & $- 2$ &     7.3890561 &  $0.037$ \\
%%\hline 5 & $+++++++++-$ &  90 & $- 6$ &   403.4287935 &  $0.434$ \\
%%\hline 6 & $++++++++++$ &   2 & $-10$ & 22026.4657948 &  $0.527$ \\
%%
%%\hline
%%\end{tabular}
%%$$
%%
%%$U = -7.95$
%%
%%XXX note that for high $\beta$, $U=-10$; for low $\beta$, $U=0$.
%%
%%%% ----------------------------------------------------------------
%%\subsection{Classical ideal gas}
%%
%%XXX two idealizations:  (1) classical/quantum; (2) ideal/real.
%%
%%XXX boltzmann statistics
%%
%%XXX here is where Widom gets into chemical complications I don't need ...
%%
%%XXX maxwell-boltzmann statistics
%%
%%XXX $N$-fold partition function decomposes as product (independence of
%%events ...)
%%
%%%% ----------------------------------------------------------------
%%\subsection{Classical real gas}
%%
%%XXX $N$-fold partition function no longer decomposes as product
%%(non-independence of events ...)
%%
%%%% ----------------------------------------------------------------
%%\subsection{Grand canonical ensemble}
%%
%%$$
%%	\Xi(T,V,\mu) = \sum_{N=0}^\infty e^{\mu\beta N} Z(T,V,N).
%%$$
%%
%%XXX grand partition function
%%
%%%% ----------------------------------------------------------------
%%\subsection{Order parameter}
%%
%%sth non-analytic at the phase xitn \ldots.



%%%% ----------------------------------------------------------------
%%\subsection{Fermi and Bose ideal gases}

%%%% ----------------------------------------------------------------
%%\subsection{The real Bose gas}



%%%% ================================================================
%%\newpage
%%\section{Bose-Einstein condensation in the ideal Bose gas}
%%\label{sec:ideal_bec}
%%
%%critical temperature
%%
%%critical density
%%
%%fraction of points in the condensation as function of $T$.

%% ================================================================
\newpage
\section{Brownian motion and Brownian bridges}

Here I recall some basic results, all in preparation for discussion of the
Feynman-Kac formulas in appendix \ref{sec:fk1} and section \ref{sec:fkboson}.
The presentation here is almost entirely an elaboration on \cite{bib:faris}.
For full information, see \cite{bib:gs}, \cite{bib:lawler}, or
\cite{bib:oksendal}.

%% ----------------------------------------------------------------
\subsection{Expectations and covariance}
\label{sec:expn_covar}
\newcommand{\fooa}{{\vecx,\vecx}}
\newcommand{\foob}{{0,T}}

\begin{notn}
\label{notn:cov}
The \plainidx{covariance} of $\R^d$-valued random variables
$$
	\vecX = \colvecthree{X_1}{\vdots}{X_d}
	\qand
	\vecY = \colvecthree{Y_1}{\vdots}{Y_d},
$$
having respective means\index{mean} $\vecmu$ and $\vecnu$, is written
$$
	\Cov(\vecX,\vecY) = \bbE[(\vecX-\vecmu)(\vecY-\vecnu)']
	= \bbE[\vecX\vecY'] - \vecmu\vecnu'.
$$
We write $\vecX'$ for the transpose of $\vecX$ from column to row vector.  Then
$\vecX\vecY'$ and $\vecmu\vecnu'$ denote the $d \times d$ \plainidx{outer
product} of $\vecX$ and $\vecY$ and $\vecmu$ and $\vecnu$, respectively, and
the expectations\index{expectation} are taken componentwise:
\begin{align*}
	\bbE[\vecX\vecY']
	&= \bbE\left[
		\colvecthree{X_1}{\vdots}{X_d}
		\rowvecthree{Y_1}{\vdots}{Y_d}
		\right]
	= \pcmatrix{
		\bbE[X_1 Y_1] & \cdots & \bbE[X_1 Y_d] \\
		\vdots        &        & \vdots        \\
		\bbE[X_d Y_1] & \cdots & \bbE[X_d Y_d] \\
		}.
\end{align*}
In case the covariance matrix has constant $c$ along the diagonal and zeroes
elsewhere, we write $\bbE[\vecX\vecY] = c I$.
\end{notn}

\begin{notn}
For
\index{$\mathbb{E}$, $\mathbb E_0^\vecx$, $\mathbb E_\fooa^\foob$}
\index{$\mathbb{P}$, $\mathbb P_0^\vecx$, $\mathbb P_\fooa^\foob$}
a \plainidx{stochastic process} $\vecX_t$, defined on, say, $t \in [0, \infty)$
or $t \in [0, T]$ for some positive real $T$, and for some Borel set $D \subset
\R^d$, we write the conditional probabilities\index{conditional probability}
$$
	\bbP_0^\veca(\vecX_t \in D) := \bbP(\vecX_t \in D \mid \vecX_0 = \veca)
$$
and
$$
	\bbP_{0,T}^{\veca,\vecb}(\vecX_t \in D)
	:= \bbP(\vecX_t \in D \mid \vecX_0 = \veca, \vecX_T = \vecb).
$$
Likewise, we write the conditional expectations\index{conditional expectation}
$$
	\bbE_0^\veca[F(\vecX_t)] := \bbE[F(\vecX_t) \mid \vecX_0 = \veca]
$$
and
$$
	\bbE_{0,T}^{\veca,\vecb}[F(\vecX_t)]
	:= \bbE[F(\vecX_t) \mid \vecX_0 = \veca, \vecX_T = \vecb].
$$
Similarly, conditional covariances\index{conditional covariance} are
$$
	\Cov_0^\veca(\vecX_t,\vecY_t)
	:= \bbE[(\vecX_t -\vecmu_t)(\vecY_t-\vecnu_t)
		\mid \vecX_0 = \veca, \vecY_0 = \veca]
$$
and
$$
	\Cov_{0,T}^{\veca,\vecb}(\vecX_t,\vecY_t)
	:= \bbE[(\vecX_t -\vecmu_t)(\vecY_t-\vecnu_t)
		\mid \vecX_0 = \veca, \vecY_0 = \veca, \vecX_T=\vecb, \vecY_T=\vecb].
$$
\end{notn}

\begin{notn}
\label{notn:expn_x1_xN}
For $N$ stochastic processes $\vecX^{(1)}, \ldots, \vecX^{(N)}$, we write the
conditional expectations\index{conditional expectation}
$$
	\bbE_{0}^{\veca^{(1)},\ldots,\veca^{(N)}}
	[F(\vecX^{(1)}_t, \ldots, \vecX^{(N)}_t)]
	:= \bbE\left[F(\vecX^{(1)}_t,\ldots,\vecX^{(N)})
		\mid
			\vecX^{(1)}_0 = \veca^{(1)}, \ldots, \vecX^{(N)}_0 = \veca^{(N)}
		\right]
$$
and
\begin{align*}
	\bbE_{0,T}^{\veca^{(1)},\vecb^{(1)}; \ldots;\veca^{(N)},\vecb^{(N)}}
	&[F(\vecX^{(1)}_t, \ldots, \vecX^{(N)}_t)] \\
	&:= \bbE\left[F(\vecX^{(1)}_t, \ldots, \vecX^{(N)}_t)
		\mid
			\vecX^{(1)}_0 = \veca^{(1)},
			\vecX^{(1)}_T = \vecb^{(1)},
			\ldots,
			\vecX^{(N)}_0 = \veca^{(N)},
			\vecX^{(N)}_T = \vecb^{(N)}
		\right].
\end{align*}
\end{notn}

%% ----------------------------------------------------------------
\subsection{Brownian motion}
\label{sec:bm}

In this section we review fundamental properties of \plainidx{Brownian motion},
highlighting its connections with the Gaussian function and laying notational
groundwork for the Feynman-Kac formulas in sections \ref{sec:fk1} and
\ref{sec:fkboson}.

Brownian motion (see \cite{bib:oksendal,bib:lawler,bib:gs} for more careful
treatments) is defined for $t \in [0, \infty) \mapsto \vecb_t \in \R^d$.  It
has the following properties:

\begin{itemize}

\item $\vecb_t$ is a \plainidx{stochastic process}, or random function of $t$,
which is almost surely continuous (with respect to the probability measure
defined below).

\item $\vecb_0 = 0$.

\item $\vecb_t$ has Gaussian distribution with \plainidx{mean} $0$ and
\plainidx{variance} $t$ in each of its $d$ components for all $t>0$, and it has
\plainidx{covariance}
$$
	\Cov(\vecb_s, \vecb_t) = \bbE[\vecb_s \vecb_t] = (s \wedge t) I
$$
where $s \wedge t$ denotes\index{wedge notation $s \wedge t$} $\min\{s, t\}$.
Specifically, the probability density function (\plainidx{PDF}) of $\vecb_t$ at
time $t$ is Gaussian with mean $0$ and component variances $t$.  For a box
$$
	D = [\ell^{(1)}, u^{(1)}] \times \cdots \times [\ell^{(d)}, u^{(d)}],
$$
we have
\begin{align}
\label{eqn:bmp1}
	\bbP(\vecb_t \in D) &= \int_D g_t(\vecy) \,d\vecy.
\end{align}
where $g_t$ is definition \ref{defn:g_t} for the Gaussian with variance $t$.

\item $\vecb_t$ has \emphidx{independent increments}:  For all $s<t<u<v$,
$\vecb_t-\vecb_s$ is independent of $\vecb_v-\vecb_u$.

\end{itemize}

It turns out (see also \cite{bib:faris}) that Brownian motion is uniquely
characterized by an extension of equation \ref{eqn:bmp1}.

\begin{prop}
\label{prop:bm_chzn}
\plainidx{Brownian motion} is uniquely characterized by the following:  for all
$n \ge 1$, all times $0 < t_1 < \ldots < t_n$, and all boxes
$$
	D_k = 
	[\ell_k^{(1)}, u_k^{(1)}]
	\times \cdots \times
	[\ell_k^{(d)}, u_k^{(d)}],
$$
\begin{equation}
\label{eqn:bmpn_bm}
\begin{aligned}
	\bbP
	\left(\bigcap_{k=1}^n \vecb_{t_k} \in D_k
	\right)
	=
		\int_{D_1}
	&
		\int_{D_2}
		\cdots
		\int_{D_n}
		\\
	&
		g_{t_1}(\vecz_1)
		g_{t_2-t_1}(\vecz_1-\vecz_2)
		\cdots
		g_{t_n-t_{n-1}}(\vecz_{n-1}-\vecz_{n})
		\, d\vecz_1 \,d\vecz_2 \cdots d\vecz_n.
\end{aligned}
\end{equation}
\end{prop}

\begin{rem}
That is, we set up some hoops and calculate the probability that a given
realization of the Brownian motion will jump through all of them (figure
\ref{fig:bm_hoops}).  The only family of functions which jump through such sets
of hoops with these normally distributed rates of success is the family of
realizations of Brownian motion.  This defines a \plainidx{probability measure}
on continuous functions $f:[0,+\infty) \to \R^d$ with $f(0)=0$.
\end{rem}

\begin{figure}[!htb]
\begin{center}
%\psfragscanon
\includegraphics{figures/comp_bm_hoops.eps}
\caption{Ten realizations of Brownian motion, moving past three boxes.
	\label{fig:bm_hoops}}
\end{center}
\end{figure}

\begin{rem*}
This is equivalent to saying that the \plainidx{joint PDF} of the random
variables $\vecb_{t_1}, \ldots, \vecb_{t_n}$ is
$$
	g_{t_1}(\vecz_1)
	g_{t_2-t_1}(\vecz_1-\vecz_2)
	\cdots
	g_{t_n-t_{n-1}}(\vecz_{n-1}-\vecz_n).
$$
\end{rem*}

\begin{rem}
\label{rem:ExFomega_bm}
Faris points out that one may generalize equation \ref{eqn:bmpn_bm}.
Rephrase it by writing
\begin{align*}
	f(\vecb_{t_1}, \ldots, \vecb_{t_n})
	&= \chi_{\cap_{k=1}^n \vecb_{t_k} \in D_k }
		(\vecb_{t_1}, \ldots, \vecb_{t_n}).
\end{align*}
Then equation \ref{eqn:bmpn_bm} becomes
\begin{equation}
\label{eqn:bmpnE_bm}
\begin{aligned}
	\bbE[f(\vecb_{t_1}, \ldots, \vecb_{t_n})] =
		\intRdntimes{d}{n}
		g_{t_1}(\vecz_1)
		&
		g_{t_2-t_1}(\vecz_1-\vecz_2)
		\cdots
		g_{t_n-t_{n-1}}(\vecz_{n-1}-\vecz_n)
		\\
	&
		f(\vecz_1, \ldots, \vecz_n)
		\; d\vecz_1 \,d\vecz_2 \cdots d\vecz_n.
\end{aligned}
\end{equation}
This holds not only for the characteristic function $\chi$ which appears in
equation \ref{eqn:bmpn_bm}, but for any Borel-measurable $f$.
\end{rem}

\begin{defn}
\label{defn:E_exp_U_bm}
We extend this to \plainidx{exponentiated integrals} by defining
\begin{equation}
\begin{aligned}
	\bbE\left[\exp\left\{\int_0^T f(\vecb_t)\,dt\right\} \right]
	:= \lim_{n\to\infty}
		\intRdntimes{d}{n}
		&
		g_{T/n}(\vecz_1)
		g_{T/n}(\vecz_1-\vecz_2)
		\cdots
		g_{T/n}(\vecz_{n-1}-\vecz_n)
		\\
	&
		\exp\left\{\frac{T}{n} \sum_{j=1}^n f(\vecz_j) \right\}
		\; d\vecz_1 \cdots d\vecz_n.
\end{aligned}
\end{equation}
\end{defn}

Now consider a realization of Brownian motion moving forward in time, having
already passed time $s$ with $\vecb_s = \vecy$ (figure \ref{fig:bm_cond}).  We
want to find the conditional density of $\vecb_t$.

\begin{figure}[!htb]
\begin{center}
%\psfragscanon
\includegraphics{figures/comp_bm_cond.eps}
\caption{Ten realizations of Brownian motion, conditioned on $\vecb_s=\vecy$.
	\label{fig:bm_cond}}
\end{center}
\end{figure}

\begin{lem}
The \plainidx{conditional density} of $\vecb_t$ given $\vecb_s = \vecy$, for $t
> s$, is
$$
	g_{t-s}(\vecz-\vecy).
$$
\begin{proof}
This is nothing more than $f_{\vecZ\mid \vecY}(\vecz\mid \vecy)$ from
elementary probability (see e.g. \cite{bib:gs}), where $\vecZ=\vecb_{t}$ and
$\vecY=\vecb_s$.  From equation \ref{eqn:bmpn_bm}, with $n=2$, the joint
density of $\vecY$ and $\vecZ$ is
$$
	g_{t-s}(\vecz-\vecy) g_s(\vecy-\vecx).
$$
We compute
\begin{align*}
	f_{\vecZ\mid \vecY}(\vecz \mid \vecy)
	&= \frac{f_{\vecZ,\vecY}(\vecz,\vecy)} {f_\vecY(\vecy)}
	= \frac{ g_{t-s}(\vecz-\vecy) g_s(\vecy-\vecx) } {g_s(\vecy-\vecx)}
	= g_{t-s}(\vecz-\vecy).
\end{align*}
\end{proof}
\end{lem}

%% ----------------------------------------------------------------
\subsection{Shifted Brownian motion}
\label{sec:shifted_bm}

Here we generalize Brownian motion (what \cite{bib:lawler} calls
\emphidx{standard Brownian motion}) by removing the requirement that the motion
start at 0.

\begin{notn}
For fixed $\vecx$, we write
\begin{equation}
\label{eqn:shifted_bm}
	\vecw^\vecx_t := \vecx + \vecb_t.
\end{equation}
If $\vecx$ is allowed to vary, we simply write $\vecw_t$.  However, this is
not well defined as a stochastic process unless we either set a probability
distribution for $\vecx$, or write $\vecw_t$ inside a conditional
expectation.  We will take the latter approach in this paper.  Also, we take
$\bbE^\vecx[\vecw_t]$ and $\bbE[\vecw^\vecx_t]$ to be equivalent.
\end{notn}

\begin{prop}
Shifted Brownian motion\index{shifted Brownian motion} has
\plainidx{conditional mean} and \plainidx{conditional covariance}
$$
	\bbE_0^\vecx[\vecw_t] = \vecx
	\qqand
	\Cov_0^\vecx[\vecw_s \vecw_t] = (s \wedge t) I.
$$
\begin{proof}
These follow immediately from equation \ref{eqn:shifted_bm}, along with the mean
and covariance for Brownian motion given in section \ref{sec:bm}.
\end{proof}
\end{prop}

Proposition \ref{prop:bm_chzn} becomes the following.

\begin{prop}
For all $n \ge 1$, all times $0 < t_1 < \ldots < t_n$, and all boxes
$$
	D_k = 
	[\ell_k^{(1)}, u_k^{(1)}]
	\times \cdots \times
	[\ell_k^{(d)}, u_k^{(d)}],
$$
we have
\begin{equation}
\label{eqn:bmpn}
\begin{aligned}
	\bbP_0^\vecx
	\left(\bigcap_{k=1}^n \vecw_{t_k} \in D_k
	\right)
	=
		\int_{D_1}
	&
		\int_{D_2}
		\cdots
		\int_{D_n}
		\\
	&
		g_{t_1}(\vecx-\vecz_1)
		g_{t_2-t_1}(\vecz_1-\vecz_2)
		\cdots
		g_{t_n-t_{n-1}}(\vecz_{n-1}-\vecz_n)
		\, d\vecz_1 \,d\vecz_2 \cdots d\vecz_n.
\end{aligned}
\end{equation}
\end{prop}

\begin{rem}
\label{rem:joint_PDF_shifted}
This is equivalent to saying that the \plainidx{joint PDF} of the random
variables $\vecw_{t_1}, \ldots, \vecw_{t_n}$, conditioned on $\vecw_0=\vecx$,
is
$$
	g_{t_1}(\vecz_1-\vecx)
	g_{t_2-t_1}(\vecz_1-\vecz_2)
	\cdots
	g_{t_n-t_{n-1}}(\vecz_{n-1}-\vecz_n).
$$
\end{rem}

\begin{rem}
\label{rem:ExFomega}
Equation \ref{eqn:bmpnE_bm} becomes
\begin{equation}
\label{eqn:bmpnE}
\begin{aligned}
	\bbE_0^\vecx[f(\vecw_{t_1}, \ldots, \vecw_{t_n})] =
		\intRdntimes{d}{n}
		g_{t_1}(\vecx-\vecz_1)
		&
		g_{t_2-t_1)}(\vecz_1-\vecz_2)
		\cdots
		g_{t_n-t_{n-1}}(\vecz_{n-1}-\vecz_n)
		\\
	&
		f(\vecz_1, \ldots, \vecz_n)
		\; d\vecz_1 \,d\vecz_2 \cdots d\vecz_n.
\end{aligned}
\end{equation}
\end{rem}

\begin{defn}
\label{defn:E_exp_U_wm}
Definition \ref{defn:E_exp_U_bm} becomes\index{exponentiated integrals}
\begin{equation}
\begin{aligned}
	\bbE_0^\vecx \left[\exp\left\{\int_0^T f(\vecw_t)\,dt\right\} \right]
	:= \lim_{n\to\infty}
		\intRdntimes{d}{n}
		&
		g_{T/n}(\vecz_1-\vecx)
		g_{T/n}(\vecz_1-\vecz_2)
		\cdots
		g_{T/n}(\vecz_{n-1}-\vecz_n)
		\\
	&
		\exp\left\{\frac{T}{n} \sum_{j=1}^n f(\vecz_j) \right\}
		\; d\vecz_1 \cdots d\vecz_n.
\end{aligned}
\end{equation}
\end{defn}

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

\cite{bib:faris} treats only bridges from $\vecx$ to $\vecx$; here, we treat
the more general case of bridges from $\vecx$ to $\vecy$.  This is needed for
\cite{bib:bu}, \cite{bib:qmath}, and in particular for section
\ref{sec:fkboson} of this paper.

\begin{defn}
The \bfidx{Brownian bridge} running from $\vecx$ at time 0 to $\vecy$ at time
$T$ is
\begin{align}
\label{eqn:bridge_process}
	\vecw_{0,t,T}^{\vecx,\vecy}
	&= \vecx + \vecb_t
		+ \frac{t}{T} \left(\vecy - \vecx - \vecb_T\right).
\end{align}
We may also write this as
\begin{align}
\label{eqn:bridge_process_aux}
	\vecw_{0,t,T}^{\vecx,\vecy}
	&= \frac{(T-t)\vecx + t\vecy}{T}
	+ \frac{T\vecb_t - t\vecb_T}{T}.
\end{align}
This looks ungainly, but it has the advantage that the \plainidx{mean} is
precisely the first term.
\end{defn}

\begin{notn}
As in section \ref{sec:shifted_bm}, if $\vecx$ and $\vecy$ are allowed to vary,
then we again have simply $\vecw_t$.  We will write $\vecw_t$
inside a conditional expectation in order to make this notation well defined as
a stochastic process.  As well, we take $\bbE^{\vecx,\vecy}_{0,T}[\vecw_t]$
and $\bbE[\vecw^{\vecx,\vecy}_{0,t,T}]$ to be equivalent.
\end{notn}

\begin{defn}
\label{defn:E_exp_U_bb}
As in definition \ref{defn:E_exp_U_wm}, we extend this to
\plainidx{exponentiated integrals} by defining
\begin{equation}
\begin{aligned}
	\bbE_{0,T}^{\vecx,\vecy}
		\left[\exp\left\{\int_0^T f(\vecw_t)\,dt\right\} \right]
	&:= \lim_{n\to\infty}
		\intRdntimes{d}{n-1} \\
%
	&
		g_{T/n}(\vecx-\vecz_1)
		g_{T/n}(\vecz_1-\vecz_2)
		\cdots
		g_{T/n}(\vecz_{n-2}-\vecz_{n-1})
		g_{T/n}(\vecz_{n-1}-\vecy)
		\\
%
	&
		\exp\left\{\frac{T}{n} \left(
			f(\vecx) + \sum_{j=1}^{n-1} f(\vecz_j) \right)
		\right\}
		\; d\vecz_1 \cdots d\vecz_{n-1} \\
%
	&:= \int \exp\left\{ \int_0^{T} f(\vecw_t)\,dt \right\}
		\,d\vecW_{0,T}^{\vecx,\vecy}(\vecw).
\end{aligned}
\end{equation}
\end{defn}

\begin{notn}
\label{notn:bridge_measure}
In the last line we have defined\index{$d\mathbf{W}$} the notation for the
bridge measure as used in \cite{bib:bu} and \cite{bib:qmath}.
\end{notn}

\begin{prop}
The bridge process $\vecw^{\vecx,\vecy}_{T,t}$ has mean\index{mean of Brownian
bridge}\index{Brownian bridge, mean and covariance of}
$$
	\vecx + \frac{t}{T} \left(\vecy-\vecx\right)
$$
and covariance\index{covariance of Brownian bridge}
$$
	I \left( s \wedge t - \frac{st}{T} \right).
$$
\begin{proof}
The mean follows immediately from taking the expectation of equation
\ref{eqn:bridge_process}.  For the covariance, we may take advantage of
equation \ref{eqn:bridge_process_aux} where the mean is readily subtractible,
and recall notation \ref{notn:cov}.  Then, since $\bbE[\vecb_s \vecb_t'] = (s
\wedge t)I$, we have
\begin{align*}
	\Cov^{\vecx,\vecy}_{0,T}[\vecw_{s},\vecw_t]
	%&= \bbE\left[ 
	%	\left(\frac{T\vecb_s - s\vecb_T}{T}\right)
	%	\left(\frac{T\vecb_t - t\vecb_T}{T}\right)'
	%	\right] \\
	&= \frac{1}{T^2} \, \bbE\left[ 
		\left(T\vecb_s - s\vecb_T\right)
		\left(T\vecb_t - t\vecb_T\right)'
		\right] \\
	&= \frac{1}{T^2} \, \left(
		T^2 \bbE[\vecb_s \vecb_t'] -
		s T \bbE[\vecb_T \vecb_t'] -
		t T \bbE[\vecb_s \vecb_T'] +
		s t \bbE[\vecb_S \vecb_T']
		\right) \\
	&= I \left( s \wedge t - \frac{st}{T} \right).
\end{align*}
\end{proof}
\end{prop}

\begin{rem}
If $s<t$ we may write this (perhaps more memorably) as
$$
	I \, \frac{s(T-t)}{T}.
$$
\end{rem}

\begin{prop}
\label{prop:two_bridges}
The \plainidx{difference of two Brownian bridges} is twice another Brownian
bridge.
\begin{proof}
Let
$$
	\vecx_1 + \vecb^{(1)}_t
		+ \frac{t}{T} \left(\vecy_1 - \vecx_1 - \vecb^{(1)}_T\right)
	\qand 
	\vecx_2 + \vecb^{(2)}_t
		+ \frac{t}{T} \left(\vecy_2 - \vecx_2 - \vecb^{(2)}_T\right)
$$
be two Brownian bridges, where $\vecb^{(1)}$ and $\vecb^{(2)}$ are independent
Brownian motions.  The difference of these bridges is
$$
	(\vecx_1 - \vecx_2) + (\vecb^{(1)}_t - \vecb^{(2)}_t)
	+ \frac{t}{T} \left((\vecy_1-\vecy_2) - (\vecx_1-\vecx_2)
	- (\vecb^{(1)}_T - \vecb^{(2)}_T)\right).
$$
This is a process running from $\vecx_1-\vecx_2$ to $\vecy_1-\vecy_2$ in time
$T$; the statistical properties of the difference of the two Brownian motions
remain to be found.

The mean of $\vecb^{(1)}_t-\vecb^{(2)}_t$ is zero by linearity of expecation.
The covariance is
\begin{align*}
	\bbE[(\vecb^{(1)}_s-\vecb^{(2)}_s) (\vecb^{(1)}_t-\vecb^{(2)}_t)]
%
	&=
	 \bbE[\vecb^{(1)}_s \vecb^{(1)}_t]
	-\bbE[\vecb^{(1)}_s \vecb^{(2)}_t]
	-\bbE[\vecb^{(2)}_s \vecb^{(1)}_t]
	+\bbE[\vecb^{(2)}_s \vecb^{(2)}_t] \\
%
	&= (s \wedge t) - 0 - 0 + (s \wedge t) = 2 (s \wedge t),
\end{align*}
where the cross terms are zero since $\vecb^{(1)}$ and $\vecb^{(2)}$ are
independent.  Reviewing the properties in section \ref{sec:bm} --- independence
of increments being trivial --- shows that the difference of the Brownian
motions is twice a Brownian motion.  Hence, the difference of the Brownian
bridges is twice a Brownian bridge.
\end{proof}
\end{prop}

%% ----------------------------------------------------------------
\subsection{Expectations over delta functions}

We continue to follow \cite{bib:faris} by computing \plainidx{expectations}
over \plainidx{delta functions}.  This may seem bizarre, but it is needed in
the proof of proposition \ref{prop:TrebH_bridge}.

\begin{lem}
\label{lem:gb0}
For the \plainidx{shifted Brownian motion} starting at $\vecx$, we have
$$
	\bbE_0^\vecx[\delta(\vecw_{T}-\vecy)] = g_{T}(\vecx-\vecy)
$$
where $g_t$ is definition \ref{defn:g_t} for the Gaussian with variance $t$.
\begin{proof}
Using remark \ref{rem:ExFomega} for the expectation, we have
\begin{align*}
	\bbE_0^\vecx[\delta(\vecw_{T}-\vecy)]
	&= \int_{\R^d}
		g_{T}(\vecx-\vecz)
		\delta(\vecz-\vecy) \,d\vecz
	= g_{T}(\vecx-\vecy).
\end{align*}
\end{proof}
\end{lem}

\begin{lem}
\label{lem:wmcovar}
The \plainidx{Wiener-measure covariance}\index{covariance} of the bridge
process $\vecw^{\vecx,\vecy}_{T,t}$ and the shifted Brownian motion $\vecw_{T}$
is zero.
\begin{proof}
\begin{align*}
	\Cov\left(\vecw^{\vecx,\vecy}_{T,t}, \vecw_{T}\right)
	&= \bbE\left[
		\left(\vecw^{\vecx,\vecy}_{T,t} - \vecx\right)\vecw_{T}' \right] \\
	&= \bbE\left[ \left(\frac{t}{T}(\vecy-\vecx) + \vecb_t
		- \frac{t}{T} \vecb_T\right) \vecb_T' \right] \\
	&= \frac{t}{T}(\vecy-\vecx) \,\bbE[\vecb_T]'
		+ \bbE[\vecb_t\vecb_T'] - \frac{t}{T} \,\bbE[\vecb_T\vecb_T'] \\
	&= \left(t - \frac{t}{T} T \right) = 0.
\end{align*}
\end{proof}
\end{lem}

The following proposition is needed for our treatment of the Feynman-Kac
formula in sections \ref{sec:fk1} and \ref{sec:fkboson_tr}.  There,
$F(\vecw)$ will be an exponentiated integral in the form of
definition \ref{defn:E_exp_U_bb}.

\begin{prop}
\label{prop:Ex_Exy}
Expectations over Brownian bridges and \plainidx{expectations} over shifted
Brownian motion are related by
$$
	\bbE_{0,T}^{\vecx,\vecy}[F(\vecw)]
	= \frac
		{\bbE_0^\vecx[F(\vecw) \delta(\vecw_{T}-\vecy)]}
		{g_{T}(\vecx-\vecy)}.
$$
\begin{proof}
We start with the right-hand side.  By construction, at time $T$ the Brownian
bridge coincides with the Brownian motion conditioned on taking value $\vecy$
at time $T$:  $\vecw_{T} = \vecy = \vecw^{\vecx,\vecy}_{T}$.  Thus
$$
	\frac
		{\bbE_0^\vecx[F(\vecw) \delta(\vecw_{T}-\vecy)]}
		{g_T(\vecx-\vecy)}
	=
	\frac
		{\bbE_0^\vecx[F(\vecw^{\vecx,\vecy}_{T}) \delta(\vecw_{T}-\vecy)]}
		{g_T(\vecx-\vecy)}.
$$
By lemma \ref{lem:wmcovar}, the Wiener-measure covariance of
$\vecw^{\vecx,\vecy}_{T,t}$ and $\vecw_{T}$ is zero.  Therefore the
expectation factors and the right-hand side becomes
$$
	\frac
		{
		\bbE_0^\vecx[F(\vecw^{\vecx,\vecy}_{T,\cdot})] \;
		\bbE_0^\vecx[\delta(\vecw_{T}-\vecy)]
		}
		{g_T(\vecx-\vecy)}.
$$
By lemma \ref{lem:gb0}, the denominator is equal to the second term in
the numerator.  Thus the right-hand side equals the left-hand side.
\end{proof}
\end{prop}

%% ----------------------------------------------------------------
\subsection{Normalized bridges}
\label{sec:normalized_bridges}

We modify definition \ref{defn:E_exp_U_bb} to obtain a \plainidx{normalized
measure}:  we set a scale factor, then prove that it is chosen correctly.

\begin{defn}
\label{defn:dWhat}
Let\index{$d\hat{\mathbf{W}}$}
\begin{equation}
	\,d\vecWhat_{0,T}^{\vecx,\vecy}(\vecw)
	=
	\frac{1}{g_{T}(\vecx-\vecy)}
	\,d\vecW_{0,T}^{\vecx,\vecy}(\vecw).
\end{equation}
\end{defn}

\begin{prop}
\label{prop:dW_normalization}
We have
$$
	\int \,d\vecWhat_{0,T}^{\vecx,\vecy}(\vecw) = 1.
$$
\begin{proof}
Starting with definition \ref{defn:E_exp_U_bb} for the non-normalized measure,
we have
\begin{align*}
	\int \,d\vecW_{0,T}^{\vecx,\vecy}(\vecw)
	&= \int \exp\left\{ \int_0^T 0 \, dt \right\}
		\,d\vecW_{0,T}^{\vecx,\vecy}(\vecw) \\
	&:= \lim_{n\to\infty}
		\intRdntimes{d}{n-1}
		g_{T/n}(\vecx-\vecz_1)
		\cdots
		g_{T/n}(\vecz_{n-1}-\vecy)
		\; d\vecz_1 \cdots d\vecz_{n-1} \\
	&= g_T(\vecx-\vecy)
\end{align*}
where the last step follows from proposition \ref{prop:g_iter}.
\end{proof}
\end{prop}

%% ----------------------------------------------------------------
\newpage
\section{Single-particle Feynman-Kac formulas}
\label{sec:fk1}

The material here is familiar from many references.  Here I review it for my
own understanding, as well as to fix notation and proof techniques which are
used for the derivation of the bosonic \plainidx{Feynman-Kac formula} in
section \ref{sec:fkboson}.

%% ----------------------------------------------------------------
\subsection{$e^{-\beta H}$ as expectation}
\label{sec:fk1_exp}

In section \ref{sec:ebH0}, for $H_0=-\nabla^2$, we saw how to interpret the
\plainidx{exponentiated operator}
$$
	e^{-\beta H_0} = e^{\beta \nabla^2}
$$
as an \plainidx{integral operator} with a \plainidx{Gaussian kernel}.  We now
ask, for $H = -\nabla^2+U(\vecx)$ where $U$ is a multiplication operator
corresponding to a potential energy, how to interpret
$$
	e^{-\beta H} = e^{-\beta (-\nabla^2 + U(\vecx))}
	= e^{\beta (\nabla^2 - U(\vecx))}.
$$

\begin{prop}
\label{prop:diffusion_FK}
For
$$
	H = -\nabla^2 + U(\vecx),
$$
we have
$$
	(e^{-\beta H} f)(\vecx)
	= \left(e^{\beta(\nabla^2 - U(\vecx))} f\right)(\vecx)
	= \bbE_0^\vecx \left[
		\exp\left\{
			-\frac{1}{2} \int_0^{2\beta} U(\vecw_s)\,ds
		\right\}
		f(\vecw_{2\beta})
	\right].
$$
\begin{proof}

The \plainidx{Trotter product formula} (see \cite{bib:simon} for a proof) says
that for self-adjoint operators $A$ and $B$,
\begin{align}
\label{eqn:trotter_AB}
	e^{\beta(A+B)}
	&= \lim_{n\to\infty} \left(e^{\beta A/n} e^{\beta B/n}\right)^n.
\end{align}
With $A=\nabla^2$ and $B=-U(\vecx)$, we have
\begin{align*}
	e^{\beta(\nabla^2-U(\vecx))}f(\vecx)
	&= \lim_{n\to\infty}
		\left(e^{\beta\nabla^2 / n} e^{-\beta U(\vecx)/n}\right)^n f(\vecx) \\
	&= \lim_{n\to\infty}
		e^{\beta\nabla^2 /n} e^{-\beta U(\vecx)/n}
		\left(e^{\beta\nabla^2 /n} e^{-\beta U(\vecx)/n}\right)^{n-1}
		f(\vecx).
\end{align*}
Recall that $e^{-\beta U(\vecx)/n}$ is simply a scalar.  Using the result of
section \ref{sec:ebH0} to write $e^{\beta \nabla^2 /n}$ as an integral
operator, we have
\begin{align*}
	(e^{-\beta H}f)(\vecx)
	&= \lim_{n\to\infty}
		\int_{\R^d} g_{2\beta/n}(\vecx-\vecz_1) 
		e^{-\beta U(\vecz_1)/n}
		\left(e^{\beta\nabla^2/n} e^{-\beta U(\vecz_1)/n}\right)^{n-1}
		f(\vecz_1)
		\,d\vecz_1.
\end{align*}
Repeating yields
\begin{align*}
	(e^{-\beta H}f)(\vecx)
	&= \lim_{n\to\infty}
		\int_{\R^d}
		\int_{\R^d}
		g_{2\beta/n}(\vecx-\vecz_1) 
		g_{2\beta/n}(\vecz_1-\vecz_2) 
	e^{-\beta U(\vecz_1)/n}
		e^{-\beta U(\vecz_2)/n}
		\left(e^{\beta\nabla^2/n} e^{-\beta U(\vecz_2)/n}\right)^{n-2}
		f(\vecz_2)
		\,d\vecz_1
		\,d\vecz_2,
\end{align*}
and finally
\begin{align*}
	(e^{-\beta H}f)(\vecx)
	= \lim_{n\to\infty}
		\intRdntimes{d}{n}
	&
		g_{2\beta/n}(\vecx-\vecz_1) 
		g_{2\beta/n}(\vecz_1-\vecz_2) 
		\cdots
		g_{2\beta/n}(\vecz_{n-1}-\vecz_n)  \\
	&
		e^{-\beta U(\vecz_1)/n}
		e^{-\beta U(\vecz_2)/n}
		\cdots
		e^{-\beta U(\vecz_n)/n}
		f(\vecz_n)
		\,d\vecz_1
		\,d\vecz_2 \dots
		\,d\vecz_n
	\\
%
	= \lim_{n\to\infty}
		\intRdntimes{d}{n}
	&
		g_{2\beta/n}(\vecx-\vecz_1) 
		g_{2\beta/n}(\vecz_1-\vecz_2) 
		\cdots
		g_{2\beta/n}(\vecz_{n-1}-\vecz_n)  \\
	&
		\exp\left\{-\frac{\beta}{n} \sum_{k=1}^n U(\vecz_k)\right\}
		f(\vecz_n)
		\,d\vecz_1
		\,d\vecz_2 \dots
		\,d\vecz_n.
\end{align*}

Now we recognize an integrand in the form of definition
\ref{defn:E_exp_U_wm},
and we can write
\begin{align*}
	(e^{-\beta H}f)(\vecx)
	&= \lim_{n \to \infty}
		\bbE_0^\vecx\left[
		\exp\left\{
			\frac{2\beta}{n}
			\left(-\frac{1}{2}\right)
			\sum_{k=1}^n U(\vecw_{2k\beta/n}) \right\}
		f(\vecw_{2\beta})
	\right].
\end{align*}
Interchanging limit and expectation by dominated convergence (since the
exponential is bounded above by 1) and recognizing the limit of the sum as a
Riemann integral, we obtain the desired result:
\begin{align*}
	(e^{-\beta H}f)(\vecx)
	&= \bbE_0^\vecx\left[
		\exp\left\{ -\frac{1}{2} \int_0^{2\beta} U(\vecw_s) \,ds \right\}
		f(\vecw_{2\beta})
	\right].
\end{align*}
\end{proof}
\end{prop}

\begin{rem}
Note that in the $U=0$ case, we have
$$
	(e^{-\beta H}f)(\vecx))
	= e^{\beta \nabla^2} f(\vecx)
	= \bbE_0^\vecx \left[ f(\vecw_{2\beta}) \right].
$$
Yet in section \ref{sec:ebH0} we also saw that
$$
	e^{\beta\nabla^2} f(\vecx)
	= \int_{\R^d} g_{2\beta}(\vecx-\vecy) f(\vecy) \,d\vecy.
$$
These two points of view are reconciled by using the \plainidx{Law of the
Unconscious Statistician}\index{Unconscious Statistician, Law of} from
elementary probability.  Namely, for a measurable function $h$ of a random
variable $Y$,
$$
	\bbE[h(\vecY)] = \int h(\vecy) f_{\vecY}(\vecy) \,d\vecy
$$
where $f_{\vecY}(\vecy)$ is the PDF for $\vecY$.  Here, for each fixed $\beta$
and $\vecx$, $\vecY$ is $\vecw_{2\beta}$, $h(\vecY)$ is $f(\vecw_{2\beta})$,
and the PDF of $\vecw_{2\beta}$ is $g_{2\beta}(\vecy-\vecx)$ as a function of
$\vecy$ as noted in remark \ref{rem:joint_PDF_shifted}.
\end{rem}

%% ----------------------------------------------------------------
\subsection{$e^{-\beta H}$ as an integral operator}
\label{sec:fk1_ker}

Section \ref{sec:operator_trace} showed how to compute $\TrebH$ if $e^{-\beta
H}$ can be represented as an \plainidx{integral operator}.  Proposition
\ref{prop:diffusion_FK} constructed the operator's kernel for the $U=0$ case.
Here we construct the \plainidx{kernel} for the general case, then apply it to
find an expression for the \plainidx{trace}.

\begin{prop}
\label{prop:ebHG}
If
\begin{equation*}
	H = -\nabla^2 + U(\vecx),
\end{equation*}
then
\begin{equation}
\label{eqn:temp1}
	e^{-\beta H} f(\vecx) = \int G_{2\beta,U}(\vecx, \vecy) f(\vecy) \,d\vecy
\end{equation}
where
\begin{equation}
\label{eqn:temp2}
	G_{2\beta,U}(\vecx,\vecy) = \bbE_0^\vecx\left[
		\exp\left\{
			-\frac{1}{2} \int_0^{2\beta} U(\vecw_s)\,ds
		\right\}
		\,\delta(\vecw_{2\beta}-\vecy)
	\right].
\end{equation}
\begin{proof}
Inserting equation \ref{eqn:temp2} into the right-hand side
of \ref{eqn:temp1} gives
\begin{align*}
	\int G_{2\beta,U}(\vecx, \vecy) f(\vecy) \,d\vecy
	&= \int
	\bbE_0^\vecx\left[
		\exp\left\{
			-\frac{1}{2} \int_0^{2\beta} U(\vecw_s)\,ds
		\right\}
		\,\delta(\vecw_{2\beta}-\vecy)
	\right]
	\,f(\vecy) \,d\vecy.
\end{align*}
Now interchange expectation and integral (by Tonelli's theorem, since the
integrand is non-negative) and integrate out the \plainidx{delta function}:
\begin{align*}
	\int G_{2\beta,U}(\vecx, \vecy) f(\vecy) \,d\vecy
	&= \bbE_0^\vecx\left[
		\int
		\exp\left\{
			-\frac{1}{2} \int_0^{2\beta} U(\vecw_s)\,ds
		\right\}
		\,\delta(\vecw_{2\beta}-\vecy)
	\,f(\vecy) \,d\vecy
	\right] \\
%
	&= \bbE_0^\vecx \left[
		\exp\left\{
			-\frac{1}{2} \int_0^{2\beta} U(\vecw_s)\,ds
		\right\}
	\,f(\vecw_{2\beta})
	\right] \\
%
	&= e^{-\beta H} f(\vecx).
\end{align*}
where the last step follows from proposition \ref{prop:diffusion_FK}.
\end{proof}
\end{prop}

%% ----------------------------------------------------------------
\subsection{$\Tr(e^{-\beta H})$ using Brownian bridges}
\label{sec:fk1_tr}

\begin{prop}
\label{prop:TrebH_bridge}
The \plainidx{trace} may be computed using \plainidx{Brownian bridges} as
follows:
\begin{align*}
	\TrebH
	&=
	g_{2\beta}(0) \int
	\bbE_{0,2\beta}^{\vecx,\vecx} \left[
			\exp\left\{
				-\frac{1}{2} \int_0^{2\beta} U(\vecw_s)\,ds
			\right\}
		\right]
	d\vecx.
\end{align*}
\begin{proof}
Using proposition \ref{prop:operator_trace}, we have
$$
	\TrebH
	=
	\int G_{2\beta,U}(\vecx,\vecx)\,d\vecx;
$$
proposition \ref{prop:ebHG} gives us an expression for $G(\vecx,\vecy)$.
Setting $\vecy=\vecx$ in equation \ref{eqn:temp2}, we have
\begin{align*}
	\TrebH
	&= \int
		\bbE_0^\vecx\left[
			\exp\left\{
				-\frac{1}{2} \int_0^{2\beta} U(\vecw_s)\,ds
			\right\}
			\,\delta(\vecw_{2\beta}-\vecx)
		\right]
	\,d\vecx.
\end{align*}
Using proposition \ref{prop:Ex_Exy}, we may convert this expectation over
\plainidx{Brownian motion} into an expectation over \plainidx{Brownian
bridges}:
\begin{align*}
	\TrebH
	&= g_{2\beta}(0)\, \int
		\bbE_{0,{2\beta}}^{\vecx,\vecx} \left[
			\exp\left\{
				-\frac{1}{2} \int_0^{2\beta} U(\vecw_s)\,ds
			\right\}
		\right]
	\,d\vecx.
\end{align*}
\end{proof}
\end{prop}

\begin{rem}
Using definition \ref{defn:E_exp_U_bb}, we may match the notation of
\cite{bib:bu} and \cite{bib:qmath}:
$$
	\TrebH =
		\int d\vecx
		\int d\vecW^{2\beta}_{\vecx\vecx}(\vecw)
		\,
		\exp \left\{- \frac{1}{2} \int_0^{2\beta} U(\vecw_s)\,ds \right\}.
$$
\end{rem}
