\chapter{Monte Carlo methods}


\subsection{Direct sampling}

Revisiting the problem in the previous chapter, we want to evaluate the expectation  
%
\begin{equation}
\label{spirit_2}
\langle \mathcal{O} \rangle = \frac{1}{Z} \sum_\mathcal{C} \mathcal{O}(\mathcal{C})  Z(\mathcal{C}) = \sum_\mathcal{C} \mathcal{O}(\mathcal{C})  p(\mathcal{C})
\end{equation}
%
where $\mathcal{O}(\mathcal{C})$ ($Z(\mathcal{C})$) is the measurement (weight) for configuration $\mathcal{C}$.
Here, $p(\mathcal{C}) \in [0,1]$ is interpreted as the probability of $\mathcal{C}$, obtained from the normalization of the configuration weight $Z(\mathcal{C})$. 
Suppose we have an {\it{erogic}} algorithm that is capable of stochastically generating any configuration $\mathcal{C}$ with equal likelihood, a direct sampling scheme can be employed to estimate $\langle \mathcal{O}\rangle$ as follows.
\\ \\
Direct sampling scheme:
\begin{enumerate}
\item Generate a random configuration $\mathcal{C}_i$.
\item Generate a random number $u_i \in [0,1]$. 
\item Perform a measurement
\begin{eqnarray}
\mathcal{O}_i = \left\{ \begin{array}{cc}  \mathcal{O}(\mathcal{C}_i)  & , \,\, u_i < p(\mathcal{C}_i)  \\ 0 & , \,\, \mathrm{otherwise}   \end{array}  \right.   \nonumber
\end{eqnarray}
\item Repeat 1-3 for $i=1,2,3,\cdots, M$
\item The {\it{law of large numbers}} in probability theory guarantees the Monte Carlo convergence of
\begin{eqnarray}
\langle \mathcal{O}\rangle_{\mathrm{MC}} = \frac{1}{M} \sum_{i=1}^M \mathcal{O}_i  \xrightarrow{ M \rightarrow \infty }  \langle \mathcal{O}\rangle \nonumber
\end{eqnarray}
towards the true expectation $\langle \mathcal{O}\rangle$ in the limit of large $M$. 
The Monte Carlo error is given by
\begin{eqnarray}
\Delta_{\mathrm{MC}} = \frac{\sigma}{\sqrt{M}} \nonumber
\end{eqnarray}
where $\sigma$ is the unbiased estimate of the standard deviation.
\end{enumerate}
%
Unfortunately, despite itself being so conceptually trivial, direct sampling schemes are unsuitable to our original problem for the following reasons.  
%
\begin{enumerate}
\item $Z$ is the partition function to be determined. Being itself unknown, $p(\mathcal{C})$ is unobtainable from the normalization of $Z(\mathcal{C})$. 
\item There may exist some configurations $\mathcal{C}$ that weigh much more than other configurations, say $\mathcal{C}'$, i.e. $Z(\mathcal{C}) \gg Z(\mathcal{C}')$. By randomly generating configurations with equal likelihood, the Monte Carlo convergence may become extremely slow.
\end{enumerate}
%
To resolve this, we should employ Markov-chain sampling schemes, as illustrated in the coming section.

\subsubsection{Law of Large Numbers}

In probability theory, Monte Carlo convergence is guaranteed by the (weak) law of large numbers.
Due to its fundamental importance, this subsection has been allocated for a quick proof.
With a little bit of algebra,  
%
\begin{eqnarray}
var(\langle \mathcal{O}\rangle_{\mathrm{MC}}) &=& \mathbf{E} \left[ \left(  \langle \mathcal{O}\rangle_{\mathrm{MC}} -   \langle \mathcal{O}\rangle \right)^2  \right]   \nonumber   \\
\frac{\sigma^2}{M} &=& \int  \left(  x -   \langle \mathcal{O}\rangle \right)^2  f_{\langle \mathcal{O}\rangle_{\mathrm{MC}}} (x)  \,\, dx   \nonumber \\
 &\geq&  \left( \int_{-\infty}^{\langle \mathcal{O}\rangle-\epsilon}  + \int_{\langle \mathcal{O}\rangle+\epsilon}^{\infty}  \right) \left(  x -   \langle \mathcal{O}\rangle \right)^2  f_{\langle \mathcal{O}\rangle_{\mathrm{MC}}} (x)  \,\, dx     \,\, , \,\,\,\,\,   \forall \epsilon > 0 \nonumber \\
 &\geq&  \epsilon^2 P( \left|  \langle \mathcal{O}\rangle_{\mathrm{MC}} -   \langle \mathcal{O} \rangle \right| > \epsilon)  \nonumber 
 \end{eqnarray}
%
where $f_{\langle \mathcal{O}\rangle_{\mathrm{MC}}} (x)$ is the probability density function of $\langle \mathcal{O}\rangle_{\mathrm{MC}}$,
we arrive at the Chebyshev's inequality or the weak law of large numbers
%
\begin{equation}
P( \left|  \langle \mathcal{O}\rangle_{\mathrm{MC}} -   \langle \mathcal{O} \rangle \right| > \epsilon) \leq \frac{\sigma^2}{M \epsilon^2}
\end{equation}
%
which guarantees the Monte Carlo convergence of $\langle \mathcal{O}\rangle_{\mathrm{MC}}$ to $\langle \mathcal{O}\rangle$ at large $M$.


\subsection{Markov chain sampling}

A Markov chain is a sequence of events, in this case configurations, i.e. $\{ \mathcal{C}_1 , \mathcal{C}_2 , \cdots , \mathcal{C}_M \}$, such that any configuration $\mathcal{C}_{i+1}$ in the sequence is in fact derived from its parent $\mathcal{C}_i$, thus preserving memory.
\\ \\
Define transition probability $\Gamma(\mathcal{C}_i \rightarrow \mathcal{C}_{i+1})$ such that
%
\begin{eqnarray}
\label{transition_probability_1}
Z(\mathcal{C}_{i+1}) &=& \sum_{\mathcal{C}_i} \Gamma(\mathcal{C}_i \rightarrow \mathcal{C}_{i+1}) Z(\mathcal{C}_i) \\
\label{transition_probability_2}
1 &=& \sum_{\mathcal{C}_i} \Gamma(\mathcal{C}_i \rightarrow \mathcal{C}_{i+1}) 
\end{eqnarray}
%
Note that the choice for $\Gamma(\mathcal{C}_i \rightarrow \mathcal{C}_{i+1})$ is not unique, so long if the above relations are satisfied.
Usually, $\Gamma(\mathcal{C}_i \rightarrow \mathcal{C}_{i+1})$ is expressed in the form of a $\Gamma$-matrix with row entry $\mathcal{C}_{i+1}$ and column entry $\mathcal{C}_i$.
Here, $\Gamma$ characterizes the entire Markov process, and the sum of each row in $\Gamma$ must be 1 in accordance to equation \ref{transition_probability_2}.


\subsubsection{Principle of ergodicity}

For any Markov process to be physical, it has to satisfy the {\it{principle of ergodicity}}, i.e. any configuration $\mathcal{C}_i$ can be derived from any other configuration $\mathcal{C}_{i-\alpha}$ in some finite $\alpha$ steps.
%
\begin{enumerate}
\item If there exists any diagonal component of $\Gamma$ to be 1, then the Markov process is not ergodic. We say that the Markov process is {\it{absorbing}}.
\item If all components of $\Gamma^\alpha$ is non-zero for some $\alpha \in \mathbf{Z}^+$, then the Markov process must be ergodic. We say that the Markov process is {\it{regular}}. In fact, all regular Markov processes are ergodic, but the reverse is not true.
\item Even if there exist any zero component in $\Gamma$, the Markov process may still be ergodic.
\end{enumerate}


\subsubsection{Principle of detailed balance}

For any Markov process to be {\it{reversible}}, it has to satisfy the {\it{principle of detailed balance}}, i.e.
%
\begin{equation}
\Gamma(\mathcal{C}_i \rightarrow \mathcal{C}_{i+1}) Z(\mathcal{C}_i) = \Gamma(\mathcal{C}_{i+1} \rightarrow \mathcal{C}_i) Z(\mathcal{C}_{i+1})
\end{equation}
%
where it automatically satisfies equation \ref{transition_probability_1}.
Here, we stress again that the choice for $\Gamma$ is not unique, and this flexibility is exactly what makes Markov chain Monte Carlo (MCMC) so attractive. 
There are 2 common choices for $\Gamma$, namely
%
\begin{enumerate}
\item heatbath transition matrix:
%
\begin{equation}
\Gamma(\mathcal{C}_i \rightarrow \mathcal{C}_{i+1}) = \frac{Z(\mathcal{C}_{i+1})}{Z}
\end{equation}
%
\item Metropolis transition matrix:
%
\begin{equation}
\Gamma(\mathcal{C}_i \rightarrow \mathcal{C}_{i+1}) = \min \left\{ 1 , \frac{Z(\mathcal{C}_{i+1})}{Z(\mathcal{C}_i)}
\right\}
\end{equation}
%
\end{enumerate}
where {\it{importance sampling}} has been intrinsically built in, thus speeding up Monte Carlo convergence. 
For our original problem, the Metropolis scheme is more beneficial because it does not require the partition function $Z$.


\subsubsection{Measurements and autocorrelations}

\begin{enumerate}
\item The {\it{law of large numbers}} in probability theory guarantees the Monte Carlo convergence of
\begin{eqnarray}
\langle \mathcal{O}\rangle_{\mathrm{MC}} = \frac{1}{M} \sum_{i=1}^M \mathcal{O}_i  \xrightarrow{ M \rightarrow \infty }  \langle \mathcal{O}\rangle \nonumber
\end{eqnarray}
towards the true expectation $\langle \mathcal{O}\rangle$ in the limit of large $M$. 
\item The Markov process possesses memory, and its history is characterized by the autocorrelation function
\begin{eqnarray}
\mathrm{Cor}_i (\mathcal{O}_i, \mathcal{O}_{i+\tau}) = \frac{ \langle \mathcal{O}_i \mathcal{O}_{i+\tau} \rangle_i  - \langle \mathcal{O}_i \rangle_i^2    }{ \langle \mathcal{O}_i^2 \rangle_i  - \langle \mathcal{O}_i \rangle_i^2 } \nonumber
\end{eqnarray} 
such that $\mathrm{Cor}_i (\mathcal{O}_i, \mathcal{O}_{i+\tau}) \in [-1,1]$. 
Here, an autocorrelation of +1 (-1) means perfect correlation (anti-correlation).
Summing over all $\tau$, we have an estimate of the (integrated) autocorrelation time
%
\begin{eqnarray}
\tau_\mathcal{O}^{\mathrm{int}} = \sum_{\tau} \mathrm{Cor}_i (\mathcal{O}_i, \mathcal{O}_{i+\tau}) \nonumber
\end{eqnarray}
%
\item The Monte Carlo error has to be corrected by
%
\begin{eqnarray}
\Delta_\mathrm{MC} = \frac{\sigma}{\sqrt{M}} \sqrt{ 1 + 2 \tau_\mathcal{O}^\mathrm{int} }  \nonumber
\end{eqnarray}
%
\item Measurements are usually expensive in Monte Carlo simulations. Therefore, it is wiser to measure only once per certain number of MC steps in order to reduce autocorrelation, or memory.
\end{enumerate}


\subsubsection{Metropolis algorithm}

\begin{enumerate}
\item From the present configuration $\mathcal{C}_i$, propose configuration $\mathcal{C}$ randomly.
\item Generate a random number $u \in [0,1]$. 
\item Accept or reject the proposal, i.e.
\begin{eqnarray}
\mathcal{C}_{i+1}  = \left\{ \begin{array}{ccc} \mathcal{C} & , &  u < \min \left\{ 1 , \frac{Z(\mathcal{C})}{Z(\mathcal{C}_i)} \right\}  \\ \mathcal{C}_i  & , & \mathrm{otherwise} \end{array} \right.
\end{eqnarray}
\item Perform a measurement every certain number of MC steps.
\end{enumerate}


\section{Quantum Monte Carlo}

\subsection{Path integral - continuous time}

Quantum Monte Carlo (QMC) is the use of Monte Carlo methods to solve complex quantum many-body problems.
Here, we restrict ourselves to continuous-time path integral QMC.
Revisiting our problem in the previous chapter, we want to random walk over all configurations $\mathcal{C} = \{ m , \, i_1 \cdots i_m \, ,  \, \tau_1\cdots \tau_m \,  | 0 < \tau_1 < \cdots < \tau_m < \beta  \}$ in a Markov chain Monte Carlo simulation, with importance sampling in accordance to the configuration weight
%
\begin{equation}
Z(\mathcal{C}) = e^{-\beta \epsilon_1} \left( e^{\tau_1 \epsilon_1} \, V_{i_1i_2}\,  e^{-\tau_1 \epsilon_2} \right) \cdots \left( e^{\tau_m \epsilon_m} \, V_{i_mi_1}\,  e^{-\tau_m \epsilon_1} \right) \,\, .
\end{equation}
%
Specifically speaking, one needs to come up with (detailed) balanced moves from one configuration to another, such that any configuration can be reached in some finite number of steps in the Markov process.
Here, we stress again that the design of such moves are not unique.
Using the same example in the previous chapter, we shall illustrate some possible moves in the following.

\subsubsection{Insertion/Annihilation of Vertex}

\begin{figure}
%\label{worldline_1d_1}
\includegraphics[width=7cm]{./figures/worldlines-1d.jpg}
\includegraphics[width=7cm]{./figures/worldlines-1d-2.jpg}
%\caption{Example of a worldline diagram of 1D boson Hubbard model with 10 sites. 
%Left: Configuration $\mathcal{C}_1$.
%Right: Configuration $\mathcal{C}'_1$.
%Horizontal axis: Site number $0,\cdots,9$. Vertical axis: Imaginary time $\tau/\beta \in [0,1)$ periodic. Dotted line: 0 particle; Solid line: 1 particle; Double solid line: 2 particles; Triple solid line: 3 particles. In this example, we set $t=1$, $U=4$, $\mu_0 = \cdots \mu_9 = 0$, $\beta = 2$.}
\end{figure}
%
See figure \ref{worldline_1d_1}. 
In this example, we propose to move from configuration $\mathcal{C}_1$ to configuration $\mathcal{C}'_1$ by annihilating a pair of vertices. 
Here, we have for configuration $\mathcal{C}_1$,
%
\begin{enumerate}
\item $\mathcal{C}_1 = \{  m = 8 \, ; \, i_1 \cdots i_8 \,;\, \tau_1 \cdots \tau_8 \,|\, 0 < \tau_1 < \cdots < \tau_8 < \beta  \}$
\item $\forall \, i,\tau \in \mathcal{C}_1$, we have
\begin{eqnarray}
\tau_1/\beta = 0.10  & :  |i_1\rangle = | 0 , 1 , 0 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_1 = 16 \nonumber \\
\tau_2/\beta = 0.30  & :  |i_2\rangle = | 0 , 0 , 1 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_2 = 16 \nonumber \\
\tau_3/\beta = 0.35  & :  |i_3\rangle = | 0 , 0 , 1 , 1 , 3 , 0 , 0 , 1 , 2 , 0 \rangle & ;  \epsilon_3 = 16 \nonumber \\
\tau_4/\beta = 0.50  & :  |i_4\rangle = | 0 , 0 , 2 , 0 , 3 , 0 , 0 , 1 , 2 , 0 \rangle & ;  \epsilon_4 = 20 \nonumber \\
\tau_5/\beta = 0.70  & :  |i_5\rangle = | 0 , 1 , 1 , 0 , 3 , 0 , 0 , 1 , 2 , 0 \rangle & ;  \epsilon_5 = 16 \nonumber \\
\tau_6/\beta = 0.75  & :  |i_6\rangle = | 0 , 1 , 1 , 0 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_6 = 16 \nonumber \\
\tau_7/\beta = 0.80  & :  |i_7\rangle = | 0 , 1 , 0 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_7 = 16 \nonumber \\
\tau_8/\beta = 0.90  & :  |i_8\rangle = | 1 , 0 , 0 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_8 = 16 \nonumber 
\end{eqnarray}
\item $\forall \, i \in \mathcal{C}_1$, we have 
\begin{eqnarray}
V_{i_1 i_2} = 1 \,\, , \,\, V_{i_2 i_3} = 1 \,\, , \,\, V_{i_3 i_4} = 1 \,\, , \,\, V_{i_4 i_5} = \sqrt{2}  \nonumber \\
V_{i_5 i_6} = \sqrt{2} \,\, , \,\, V_{i_6 i_7} = 1 \,\, , \,\, V_{i_7 i_8} = 1 \,\, , \,\, V_{i_8 i_1} = 1 \nonumber
\end{eqnarray}
\end{enumerate}
%
and for configuration $\mathcal{C}'_1$,
%
\begin{enumerate}
\item $\mathcal{C}'_1 = \{  m = 6 \, ; \, i'_1 \cdots i'_6 \,;\, \tau'_1 \cdots \tau'_6 \,|\, 0 < \tau'_1 < \cdots < \tau'_6 < \beta  \}$
\item $\forall \, i',\tau' \in \mathcal{C}'_1$, we have
\begin{eqnarray}
\tau'_1/\beta = 0.10  & :  |i'_1\rangle = | 0 , 1 , 0 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon'_1 = 16 \nonumber \\
\tau'_2/\beta = 0.30  & :  |i'_2\rangle = | 0 , 0 , 1 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon'_2 = 16 \nonumber \\
\tau'_3/\beta = 0.50  & :  |i'_3\rangle = | 0 , 0 , 2 , 0 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon'_3 = 20 \nonumber \\
\tau'_4/\beta = 0.70  & :  |i'_4\rangle = | 0 , 1 , 1 , 0 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon'_4 = 16 \nonumber \\
\tau'_5/\beta = 0.80  & :  |i'_5\rangle = | 0 , 1 , 0 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon'_5 = 16 \nonumber \\
\tau'_6/\beta = 0.90  & :  |i'_6\rangle = | 1 , 0 , 0 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon'_6 = 16 \nonumber 
\end{eqnarray}
\item $\forall \, i' \in \mathcal{C}'_1$, we have 
\begin{eqnarray}
V_{i'_1 i'_2} = 1 \,\, , \,\, V_{i'_2 i'_3} = 1 \,\, , \,\, V_{i'_3 i'_4} = \sqrt{2}  \nonumber \\
V_{i'_4 i'_5} = \sqrt{2} \,\, , \,\, V_{i'_5 i'_6} = 1 \,\, , \,\, V_{i'_6 i'_1} = 1 \nonumber
\end{eqnarray}
\end{enumerate}
%
Suppose we use Metropolis algorithm to balance the reverse move (creating the same pair of vertices), we could derive the following acceptance probabilities for this particular example.
%
\begin{eqnarray}
p(\mathcal{C}_1 \rightarrow \mathcal{C}'_1) &=&  1     \nonumber \\
p(\mathcal{C}'_1 \rightarrow \mathcal{C}_1) &=&  1     \nonumber 
\end{eqnarray}


\subsubsection{Move Vertex}

\begin{figure}
%\label{worldline_1d_2}
\includegraphics[width=7cm]{./figures/worldlines-1d-2.jpg}
\includegraphics[width=7cm]{./figures/worldlines-1d-3.jpg}
%\caption{Example of a worldline diagram of 1D boson Hubbard model with 10 sites. 
%Left: Configuration $\mathcal{C}_2$.
%Right: Configuration $\mathcal{C}'_2$.
%Horizontal axis: Site number $0,\cdots,9$. Vertical axis: Imaginary time $\tau/\beta \in [0,1)$ periodic. Dotted line: 0 particle; Solid line: 1 particle; Double solid line: 2 particles; Triple solid line: 3 particles. In this example, we set $t=1$, $U=4$, $\mu_0 = \cdots \mu_9 = 0$, $\beta = 2$.}
\end{figure}
%
See figure \ref{worldline_1d_2}. 
In this example, we propose to move the vertex over imaginary time from configuration $\mathcal{C}_2$ to $\mathcal{C}'_2$.
Using Metropolis algorithm to balance the reverse move, we can easily obtain the acceptance probabilities for this particular example.
%
\begin{eqnarray}
p(\mathcal{C}_2 \rightarrow \mathcal{C}'_2) &=&  \min \left\{ 1, \frac{e^{(0.35)(-4)}}{e^{{(0.5)(-4)}}}  \right\}  = 1     \nonumber \\
p(\mathcal{C}'_2 \rightarrow \mathcal{C}_2) &=&  \min \left\{ 1, \frac{e^{(0.5)(-4)}}{e^{{(0.35)(-4)}}}  \right\}  = e^{-0.6} \nonumber 
\end{eqnarray}


\subsubsection{Insertion/Removal of Particles}

See figure \ref{worldline_1d_3}.
In this example, we propose to remove a particle from configuration $\mathcal{C}_3$ to configuration $\mathcal{C}'_3$.
Here, we have for configuration $\mathcal{C}_3$
%
\begin{enumerate}
\item $\mathcal{C}_3 = \{  m = 6 \, ; \, i_1 \cdots i_6 \,;\, \tau_1 \cdots \tau_6 \,|\, 0 < \tau_1 < \cdots < \tau_6 < \beta  \}$
\item $\forall \, i,\tau  \in \mathcal{C}_3$, we have
\begin{eqnarray}
\tau_1/\beta = 0.10  & :  |i_1\rangle = | 0 , 1 , 0 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_1 = 16 \nonumber \\
\tau_2/\beta = 0.30  & :  |i_2\rangle = | 0 , 0 , 1 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_2 = 16 \nonumber \\
\tau_3/\beta = 0.35  & :  |i_3\rangle = | 0 , 0 , 2 , 0 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_3 = 20 \nonumber \\
\tau_4/\beta = 0.70  & :  |i_4\rangle = | 0 , 1 , 1 , 0 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_4 = 16 \nonumber \\
\tau_5/\beta = 0.80  & :  |i_5\rangle = | 0 , 1 , 0 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_5 = 16 \nonumber \\
\tau_6/\beta = 0.90  & :  |i_6\rangle = | 1 , 0 , 0 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_6 = 16 \nonumber 
\end{eqnarray}
\item $\forall \, i \in \mathcal{C}_3$, we have 
\begin{eqnarray}
V_{i_1 i_2} = 1 \,\, , \,\, V_{i_2 i_3} = 1 \,\, , \,\, V_{i_3 i_4} = \sqrt{2}  \nonumber \\
V_{i_4 i_5} = \sqrt{2} \,\, , \,\, V_{i_5 i_6} = 1 \,\, , \,\, V_{i_6 i_1} = 1 \nonumber
\end{eqnarray}
\end{enumerate}
%
and for configuration $\mathcal{C}'_3$,
%
\begin{enumerate}
\item $\mathcal{C}'_1 = \{  m = 6 \, ; \, i'_1 \cdots i'_6 \,;\, \tau'_1 \cdots \tau'_6 \,|\, 0 < \tau'_1 < \cdots < \tau'_6 < \beta  \}$
\item $\forall \, i',\tau' \in \mathcal{C}'_3$, we have
\begin{eqnarray}
\tau'_1/\beta = 0.10  & :  |i'_1\rangle = | 0 , 1 , 0 , 1 , 2 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon'_1 = 8 \nonumber \\
\tau'_2/\beta = 0.30  & :  |i'_2\rangle = | 0 , 0 , 1 , 1 , 2 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon'_2 = 8 \nonumber \\
\tau'_3/\beta = 0.35  & :  |i'_3\rangle = | 0 , 0 , 2 , 0 , 2 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon'_3 = 12 \nonumber \\
\tau'_4/\beta = 0.70  & :  |i'_4\rangle = | 0 , 1 , 1 , 0 , 2 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon'_4 = 8 \nonumber \\
\tau'_5/\beta = 0.80  & :  |i'_5\rangle = | 0 , 1 , 0 , 1 , 2 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon'_5 = 8 \nonumber \\
\tau'_6/\beta = 0.90  & :  |i'_6\rangle = | 1 , 0 , 0 , 1 , 2 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon'_6 = 8 \nonumber 
\end{eqnarray}
\item $\forall \, i' \in \mathcal{C}'_3$, we have 
\begin{eqnarray}
V_{i'_1 i'_2} = 1 \,\, , \,\, V_{i'_2 i'_3} = 1 \,\, , \,\, V_{i'_3 i'_4} = \sqrt{2}  \nonumber \\
V_{i'_4 i'_5} = \sqrt{2} \,\, , \,\, V_{i'_5 i'_6} = 1 \,\, , \,\, V_{i'_6 i'_1} = 1 \nonumber
\end{eqnarray}
\end{enumerate}
%
Using Metropolis algorithm to balance the reverse move of inserting the same particle, we can easily derive the acceptance probabilities for this particular example.
%
\begin{eqnarray}
p(\mathcal{C}_3 \rightarrow \mathcal{C}'_3) &=&  \min \left\{ 1, \frac{e^{-\beta \epsilon'_1}}{e^{-\beta \epsilon_1}}  \right\}  = 1     \nonumber \\
p(\mathcal{C}'_3 \rightarrow \mathcal{C}_3) &=&  \min \left\{ 1, \frac{e^{-\beta \epsilon_1}}{e^{-\beta \epsilon'_1}}  \right\}  = e^{-16} \nonumber 
\end{eqnarray}
%
\begin{figure}
%\label{worldline_1d_3}
\includegraphics[width=7cm]{./figures/worldlines-1d-3.jpg}
\includegraphics[width=7cm]{./figures/worldlines-1d-4.jpg}
%\caption{Example of a worldline diagram of 1D boson Hubbard model with 10 sites. 
%Left: Configuration $\mathcal{C}_3$.
%Right: Configuration $\mathcal{C}'_3$.
%Horizontal axis: Site number $0,\cdots,9$. Vertical axis: Imaginary time $\tau/\beta \in [0,1)$ periodic. Dotted line: 0 particle; Solid line: 1 particle; Double solid line: 2 particles; Triple solid line: 3 particles. In this example, we set $t=1$, $U=4$, $\mu_0 = \cdots \mu_9 = 0$, $\beta = 2$.}
\end{figure}
%


\subsubsection{Non-ergodicity issue}

With only the moves mentioned above, the Markov process is solely capable of random-walking over a class of configurations, which is unfortunately not universal.
For example, the Markov process is not capable of creating a particle that winds over a periodic system (i.e. horizontally from one boundary to another) starting from either configuration $\mathcal{C}_1$, $\mathcal{C}_2$ or $\mathcal{C}_3$.
Unless one can "cut" some closed lines and move these segments randomly before "joining" them eventually, the Markov process will never be ergodic with only the moves mentioned above.
\\ \\
Having realised this, Prokof'ev and Svistunov proposed the {\it{worm algorithm}} in 1998, which is capable of sampling all configurations not only ergodically, but also efficiently.
This will be the subject matter of the next section.

\subsection{Worm algorithm}

Starting from any worldline configuration, the worm algorithm begins by either inserting or deleting a line segment over an infinitesimal time at any random site and random time.
The two ends of the line segment are called worms.
By fixing one of the worms (wormtail) and stochastically moving the other (wormhead)
%
\begin{enumerate}
\item over time to create/annihilate particles; or,
\item across sites to insert/remove vertices, 
\end{enumerate}
%
until eventually removing the wormpair upon touching the wormhead with the wormtail, the worm algorithm generates an ergodic Markov chain of configurations.  
Before we write down the detailed balanced equations for the stochastic movement of the wormhead based on Pollet's implementation of directed worm algorithm\cite{Pollet}, we need a proper theoretical formalism of worms which is the subject matter of the coming section.

\subsubsection{Worms, Extended configurations, and Green's functions}

\begin{figure}
\centering
%\label{worldline_1d_5}
\includegraphics[width=9cm]{./figures/worldlines-1d-5.jpg}
\caption{Example of an extended worldline diagram of 1D boson Hubbard model with 10 sites. 
Horizontal axis: Site number $0,\cdots,9$. Vertical axis: Imaginary time $\tau/\beta \in [0,1)$ periodic. Dotted line: 0 particle; Solid line: 1 particle; Double solid line: 2 particles; Triple solid line: 3 particles. In this example, we set $t=1$, $U=4$, $\mu_0 = \cdots \mu_9 = 0$, $\beta = 2$.}
\end{figure}


A sole purpose of worms is to either insert or delete line segments over imaginary times. Therefore, worms are simply creation ($\hat{b}^\dag$) or annihilation ($\hat{b}$) operators.
During the stochastic process, the worms can be anywhere in the worldline diagram, "penetrating" an existing configuration $\mathcal{C}$.
We call such configuration in the presence of worms an extended configuration $\mathcal{E}$.
Without loss of generality, consider an extended configuration $\mathcal{E}_{pq} = \{ m \, , \, i_1 \cdots i_p \cdots i_q \cdots i_m \, , \, \tau_1 \cdots \tau_p \cdots \tau_q \cdots \tau_m \, | \,  0 < \tau_1 < \cdots < \tau_p < \cdots < \tau_q < \cdots < \tau_m < \beta \} $ of $m$ vertices at imaginary times $\tau_1 \cdots \tau_m$, and the annihilation (creation) operator at time $\tau_p$ ($\tau_q$).
The extended configuration weight reads
%
\begin{equation}
\label{extended_configuration_weight}
Z(\mathcal{E}_{pq}) = e^{-\beta \epsilon_1} \left( e^{-\tau_1 \epsilon_1} \, V_{i_1i_2}\,  e^{\tau_1 \epsilon_2} \right) \cdots \left( e^{-\tau_p \epsilon_p} \, b_p\,  e^{\tau_p \epsilon_p'} \right) \cdots \left( e^{-\tau_q \epsilon_q} \, b_q^\dag\,  e^{\tau_q \epsilon_q'} \right) \cdots  \left( e^{-\tau_m \epsilon_m} \, V_{i_mi_1}\,  e^{\tau_m \epsilon_1} \right)   
\end{equation}
%
which sums up to the Green's function after normalisation
%
\begin{eqnarray}
G (x_p , \tau_p \,;\, x_q , \tau_q) &=&  \frac{1}{Z}  \sum_{m=0}^\infty \, \sum_{i_1 \cdots i_m } \sum_{i_p,i_q} e^{-\beta \epsilon_1} \int_0^\beta d\tau_1 \cdots \int_0^{\tau_{m-1}} d\tau_m \, \nonumber \\ 
 & & \times \left( e^{-\tau_1 \epsilon_1} \, V_{i_1i_2}\,  e^{\tau_1 \epsilon_2} \right) \cdots  \nonumber \\
 & & \times \cdots \left( e^{-\tau_p \epsilon_p} \, b_p\,  e^{\tau_p \epsilon_p'} \right) \cdots \left( e^{-\tau_q \epsilon_q} \, b_q^\dag\,  e^{\tau_q \epsilon_q'} \right) \cdots \nonumber \\
 & & \times  \cdots \left( e^{-\tau_m \epsilon_m} \, V_{i_mi_1}\,  e^{\tau_m \epsilon_1} \right)  
\end{eqnarray}
%
Figure \ref{worldline_1d_5} illustrates an example of an extended worldline diagram for the case of 1D systems.
In this example, the extended worldline diagram reads the following information:
%
\begin{enumerate}
\item $\mathcal{E}_{pq} = \{  m = 7 \, ; \, i_1 \cdots i_p \cdots i_q \cdots i_7 \,;\, \tau_1 \cdots \tau_p \cdots \tau_q \cdots \tau_7  \}$
\item $\forall \, i,\tau \in \mathcal{E}_{pq}$, we have
\begin{eqnarray}
\tau_1/\beta = 0.10  & :  |i_1\rangle = | 0 , 1 , 0 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_1 = 16 \nonumber \\
\tau_q/\beta = 0.20  & :  |i_q\rangle = | 0 , 1 , 0 , 1 , 3 , 1 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_q = 16 \nonumber \\
\tau_2/\beta = 0.30  & :  |i_2\rangle = | 0 , 0 , 1 , 1 , 3 , 1 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_2 = 16 \nonumber \\
\tau_3/\beta = 0.40  & :  |i_3\rangle = | 0 , 0 , 1 , 1 , 3 , 0 , 2 , 0 , 2 , 0 \rangle & ;  \epsilon_3 = 20 \nonumber \\
\tau_4/\beta = 0.50  & :  |i_4\rangle = | 0 , 0 , 2 , 0 , 3 , 0 , 2 , 0 , 2 , 0 \rangle & ;  \epsilon_4 = 24 \nonumber \\
\tau_p/\beta = 0.60  & :  |i_p\rangle = | 0 , 0 , 2 , 0 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_p = 20 \nonumber \\
\tau_5/\beta = 0.70  & :  |i_5\rangle = | 0 , 1 , 1 , 0 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_5 = 16 \nonumber \\
\tau_6/\beta = 0.80  & :  |i_6\rangle = | 0 , 1 , 0 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_6 = 16 \nonumber \\
\tau_7/\beta = 0.90  & :  |i_7\rangle = | 1 , 0 , 0 , 1 , 3 , 0 , 1 , 0 , 2 , 0 \rangle & ;  \epsilon_7 = 16 \nonumber 
\end{eqnarray}
\item $\forall \, i \in \mathcal{E}_{pq}$, we have 
\begin{eqnarray}
V_{i_1 i_q} = 1 \,\, , \,\, b_q = 1 \,\, , \,\, V_{i_2 i_3} = 1 \nonumber \\ 
V_{i_3 i_4} = \sqrt{2}  \,\, , \,\, V_{i_4 i_p} = \sqrt{2} \,\, , \,\, b_p^\dag = \sqrt{2} \nonumber \\
V_{i_5 i_6} = \sqrt{2} \,\, , \,\, V_{i_6 i_7} = 1 \,\, , \,\, V_{i_7 i_1} = 1 \nonumber
\end{eqnarray}
\end{enumerate}
%

\subsubsection{Density Matrices}

Density matrices are Green's functions at equal time, ie.
%
\begin{equation}
\langle b_q^\dag b_p \rangle = G(x_p , \tau \,;\, x_q , \tau)
\end{equation}
%
for any $\tau$.
They are useful theoretical quantities for the computation of experimentally relevant physical observables, such as lattice momentum density 
%
\begin{equation}
\langle n_{\bf{k}} \rangle = \langle \hat{b}^\dag_{\bf{k}} \hat{b}_{\bf{k}} \rangle = \frac{1}{N} \sum_{pq} \langle \hat{b}_q^\dag \hat{b}_p \rangle e^{i{\bf{k}}\cdot ({\bf{x}}_q - {\bf{x}}_p )} \,\, ,
\end{equation}
%
where N is the total number of bosons, and the local density
%
\begin{equation}
\langle n_p \rangle = \langle b_p^\dag b_p \rangle
\end{equation}
%
coincides with the diagonal density matrix element.
In practice, decoherence effects can be theoretically observed in the rapid vanishing of the offdiagonal density matrix elements.
\\ \\
The implementation of density matrices will be left to a later section.


