\chapter{Quantum Monte Carlo}
\label{chapter_QMC}

This chapter introduces an efficient Quantum Monte Carlo method \cite{prokofev1}, implemented via the directed worm algorithm \cite{pollet1}, that is used to solve the boson Hubbard model  \ref{boson_hubbard_model} exact numerically at finite temperature.


\section{Quantum statistical mechanics}

\subsection{Feynman path-integral formalism}

Rewriting the boson Hubbard hamiltonian \ref{boson_hubbard_model} as $\hat{H} - \mu \hat{N} = \hat{H}_0 - \hat{V}$ such that 
%
\begin{eqnarray}
\hat{H_0} &=& \frac{U}{2} \sum_i n_i(n_i-1) - \sum_i \mu_i n_i  \\
\hat{V}      &=& t \sum_{\langle i,j \rangle} \hat{b}^\dag_i  \hat{b}_j \,\, ,
\end{eqnarray}
%
the equation of motion for the evolution operator $\hat{U}(\tau,\tau_0)$ becomes
%
\begin{equation}
\label{eom}
\frac{\partial}{\partial \tau} \hat{U}(\tau,\tau_0) = \hat{V}(\tau) \hat{U}(\tau,\tau_0)
\end{equation}
%
in the interaction picture, taking $\hbar=1$ \cite{fetter}.
With
%
\begin{eqnarray}
\label{evolution_operator}
\hat{U}(\tau,\tau_0) &=& e^{\hat{H}_0 \tau} \, e^{-\hat{H} (\tau-\tau_0)} \, e^{-\hat{H}_0 \tau_0} \\
\hat{V}(\tau) &=& e^{\hat{H}_0 \tau} \, \hat{V} \, e^{-\hat{H}_0 \tau} \,\, ,
\end{eqnarray}
%
equation (\ref{eom}) becomes
%
\begin{equation}
\hat{U}(\tau,\tau_0) - \hat{U}(\tau_0,\tau_0) = \int_{\tau_0}^{\tau} d\tau_1  \hat{V}(\tau_1) \hat{U}(\tau_1,\tau_0)
\end{equation}
%
or ($\tau=\beta$, $\tau_0=0$),
%
\begin{eqnarray}
\hat{U}(\beta) &=& 1 + \int_0^\beta d\tau_1  \hat{V}(\tau_1) \hat{U}(\tau_1) \nonumber \\
    &=& 1 +  \int_0^\beta d\tau_1 \,  \hat{V}(\tau_1) + \int_0^\beta d\tau_1 \int_0^{\tau_1} d\tau_2 \, \hat{V}(\tau_1) \hat{V}(\tau_2) + \cdots
\end{eqnarray}
%
According to equation (\ref{evolution_operator}), we have $\hat{U}(\beta) = e^{\beta\hat{H}_0} \, e^{-\beta\hat{H}}$, thus the partition function
%
\begin{equation}
Z = \mathrm{Tr} \, e^{-\beta \hat{H}} = \mathrm{Tr} \, e^{-\beta \hat{H}_0} \, \hat{U}(\beta)
\end{equation}
%
In Fock basis $\{ | i \rangle \}$, the partition function can be written as
%
\begin{equation}
\label{partition_function_integral}
Z = \sum_{m=0}^\infty \, \sum_{i_1 \cdots i_m } e^{-\beta \epsilon_1} \int_0^\beta d\tau_1 \cdots \int_0^{\tau_{m-1}} d\tau_m \, \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}
%
where $\hat{H}_0 |i\rangle = \epsilon_i |i\rangle$ and $V_{ij} = \langle i | \hat{V} | j \rangle$. 


\subsection{Configuration and worldlines diagram}

The integrand of the partition function Z in equation \ref{partition_function_integral} is uniquely determined if the {\it{configuration}} $\mathcal{C} = \{ m , \, i_1 \cdots i_m \, ,  \, \tau_1\cdots \tau_m \,  | 0 < \tau_1 < \cdots < \tau_m < \beta  \}$ is known.
Defining the configuration weight
%
\begin{equation}
\label{configuration_weight}
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}
%
we have now a simpler notation
%
\begin{eqnarray}
\label{spirit1}
Z &=& \sum_\mathcal{C} Z(\mathcal{C})   \\
\label{spirit2}
\langle \mathcal{O} \rangle &=& \frac{1}{Z} \sum_\mathcal{C} \mathcal{O}(\mathcal{C})  Z(\mathcal{C}) 
\end{eqnarray}
%
for diagonal observables $\langle \mathcal{O} \rangle = \frac{1}{Z} \, \mathrm{Tr} \, \hat{\mathcal{O}} e^{-\beta \hat{H}}$ \cite{prokofev1}.
Here, $\mathcal{O}(\mathcal{C})$ is called the measurement of $\hat{\mathcal{O}}$ for configuration $\mathcal{C}$.
In practice, configurations are visualized as worldlines diagrams. 
Figure \ref{worldline_1d} illustates an example of a worldlines diagram that corresponds to the configuration $\mathcal{C} = \{  m = 8 \, ; \, i_1 \cdots i_8 \,;\, \tau_1 \cdots \tau_8 \,|\, 0 < \tau_1 < \cdots < \tau_8 < \beta  \}$.
In this example, 
%
\begin{figure}[h]
\label{worldline_1d}
\centering
\includegraphics[width=10cm]{./figures/worldlines-1d.jpg}
\caption{Example of a worldlines 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}
%
\begin{enumerate}
\item $\forall \, i,\tau \in \mathcal{C}$, 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}$, 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}
\item For some diagonal observables $\mathcal{O}$, we have
\begin{eqnarray}
N  &:& N(\mathcal{C}) = 8 \nonumber \\
E_0 &:& E_0 (\mathcal{C}) = 16 \nonumber 
\end{eqnarray}
\end{enumerate}
%


\section{Markov chain sampling}

The whole spirit of Quantum Monte Carlo (QMC) simulation is to obtain statistics for $\langle \mathcal{O} \rangle$ in equation \ref{spirit2} through random walking in the configuration space efficiently according to some stochastic algorithm in mind.
A Markov chain random walk is a sequence of 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 \cite{markov_chain}.
With importance sampling in accordance to the configuration weight \ref{configuration_weight}, the Monte Carlo estimate
%
\begin{equation}
\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{equation}
%
converges to $\langle \mathcal{O} \rangle$ with increasing statistics\footnote{This convergence is guaranteed by the law of large numbers \cite{probability_theory}.}.
The stochastic algorithm which governs the Markov chain random walk does not need to be unique, so long it satisfies both the principle of ergodicity, as well as the principle of detailed balance \cite{markov_chain}
%
\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}
%
Here, the transition matrix $\Gamma$ has been chosen to follow Metropolis, i.e.
%
\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}
%
which satisfies detailed balance automatically \cite{metropolis}.


\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\footnote{The frequency of measurements should be adjusted according to the autocorrelation time in the measurments.}.
\end{enumerate}

\newpage


\section{Directed worm algorithm}

\subsection{Local updates and non-ergodicity}

\begin{figure}[h]
\includegraphics[width=6cm]{./figures/worldlines-1d.jpg}
\includegraphics[width=6cm]{./figures/worldlines-1d-2.jpg}
\caption{Left $\rightarrow$ ($\leftarrow$) Right : insertion (deletion) of vertices}
\label{worldline_1d_1}
\end{figure}
%
\begin{figure}[h]
\includegraphics[width=6cm]{./figures/worldlines-1d-2.jpg}
\includegraphics[width=6cm]{./figures/worldlines-1d-3.jpg}
\caption{Left $\rightarrow$ ($\leftarrow$) Right : move vertex over time}
\label{worldline_1d_2}
\end{figure}
%
\begin{figure}[h]
\includegraphics[width=6cm]{./figures/worldlines-1d-3.jpg}
\includegraphics[width=6cm]{./figures/worldlines-1d-4.jpg}
\caption{Left $\rightarrow$ ($\leftarrow$) Right : delete (insert) 1 particle}
\label{worldline_1d_3}
\end{figure}
%
\noindent \\ 
\noindent Figure \ref{worldline_1d_1}, \ref{worldline_1d_2}, and \ref{worldline_1d_3} illustrate some possible moves to update configurations in the Markov chain random walk.
These moves result in local updates of closed worldlines that is unfortunately not universal.
For instance, 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 any configuration. 
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 these moves.
Having realised this, Prokof'ev and Svistunov proposed the {\it{worm algorithm}} in 1998 \cite{prokofev1}, which is capable of sampling all configurations ergodically and efficiently.
After nearly a decade, Pollet proposed the {\it{directed worm algorithm}} \cite{pollet1} which is a faster implementation compared to the original version.



\subsection{Worms and extended configurations}

\begin{figure}[h]
\centering
\includegraphics[width=9cm]{./figures/worldlines-1d-5.jpg}
\caption{Example of an extended worldlines 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$.}
\label{worldline_1d_5}
\end{figure}

\noindent An open worldline has 2 end points, called {\it{worms}} \cite{prokofev1}, as illustrated by a blue and a red dot in figure \ref{worldline_1d_5}.
The blue dot, called the {\it{wormtail}}, is fixed at a specific location in space-time, whereas the red dot, called the {\it{wormhead}}, moves around.
In the directed worm algorithm, the wormhead moves in a specific direction until bouncing occurs.
In fact, worms are either $\hat{b}(\tau_p)$ or $\hat{b}^\dag(\tau_q)$ operators, that extend a configuration to an {\it{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 is naturally defined as
%
\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 \cite{fetter} after normalization
%
\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}
%
Density matrices are Green's functions at equal time, i.e.
%
\begin{equation}
\langle b_q^\dag b_p \rangle = G(x_p , \tau \,;\, x_q , \tau) \,\, .
\end{equation}
%

\subsection{Updates in extended configuaration space}

\begin{figure}[h]
\centering
\label{extended_worldline_1}
\includegraphics[width=7cm]{./figures/extended-worldlines-1.jpg}
\end{figure}

\noindent We begin studying detailed balance by looking at the movement of wormhead from extended configuration X to Y as illustrated above.
From X to Y, the $\hat{b}^\dag$ wormhead 
%
\begin{enumerate}
\item moves forward in time from $\tau_p$ to $\tau_v$ unhalted,
\item jumps from site i to site j, thus creating a $\hat{b}_i^\dag \hat{b}_j$ vertex at time $\tau_v$, 
\item moves forward in time from $\tau_v^+$ to $\tau_{p'}$ halted. 
\end{enumerate}
%
The reverse move, from Y to X, invokes the $\hat{b}^\dag$ wormhead to
%
\begin{enumerate}
\item move backward in time from $\tau_{p'}$ to $\tau_v^+$ halted,
\item jump from site j to site i, thus annihilating the $\hat{b}_i^\dag \hat{b}_j$ vertex at time $\tau_v$,
\item move backward in time from $\tau_v^-$ to $\tau_p$ halted.
\end{enumerate}
%
According to equation \ref{extended_configuration_weight}, we have
%
\begin{equation}
\frac{Z(Y)}{Z(X)} = \frac{ \left( e^{-\tau_{p'} \epsilon_v} \, \langle i'_p | \hat{b}_j^\dag | i_v \rangle   \,  e^{\tau_{p'} \epsilon_p'} \right)   \left( e^{-\tau_v \epsilon_p} \, \langle  i'_p | t \, \hat{b}_i^\dag \hat{b}_j  |i_v  \rangle \,  e^{\tau_v \epsilon_v} \right)  } { \left( e^{-\tau_p \epsilon_p} \, \langle  i'_p | \hat{b}_i^\dag | i_p \rangle   \,  e^{\tau_p \epsilon_p'} \right) } \,\, .
\end{equation}
%
Since $Z(X) P(X \rightarrow Y) = Z(Y) P(Y \rightarrow X)$, i.e.
%
\begin{equation}
\frac{P(X \rightarrow Y)}{P(Y \rightarrow X)} = \left( \frac{ \epsilon_p \, e^{-\epsilon_p \tau_{vp}} \, e^{-\epsilon_v \tau_{p'v}}  }{  e^{-\epsilon'_p \tau_{p'v} }  \, e^{-\epsilon'_p \tau_{vp} }   }  \right) \times \left( \frac{ t \langle i'_p  | \hat{b}_j^\dag | i_v  \rangle \langle i'_p |  \hat{b}_i^\dag \hat{b}_j  | i_v  \rangle }{ \epsilon_p \langle  i'_p| \hat{b}_i^\dag | i_p \rangle }   \right) 
\end{equation}
%
(eg. $\tau_{vp} = \tau_v - \tau_p$), what remains is simply to allocate the balanced probabilities for the sub-moves.



\subsubsection{Unhalted move in time}

\begin{figure}[h]
\centering
\label{extended_worldline_2}
\includegraphics[width=4cm]{./figures/extended-worldlines-2.jpg}
\end{figure}

\noindent The first sub-move involves the movement of wormhead over time, as illustrated in the above figure.
Here, the transition probability from extended configuration $X_1$ to $X_2$ is set as
%
\begin{equation}
\label{exponential_distribution}
P(X_1 \rightarrow X_2) = \epsilon_p e^{-\epsilon_p \tau_{p'p} }   \,\, .
\end{equation}
%
\begin{enumerate}
\item By unhalted move over time itself, it does not satisfy detailed balance. To balance it, the unhalted move must be immediately followed by an attempt to create a vertex at that new time $\tau_{p'}$. (See later section.)
\item In probability theory, the event of an attempt to create a vertex is naturally regarded as a Poisson process, a direct consequence of definition \ref{exponential_distribution}, ie.  the time interval until the random new time $\tau_{p'}$ follows an exponential distribution. 
\item For the transition probability to make sense, the mean occurrence of such a Poisson event is $\epsilon_p$ must be strictly positive, ie. $\epsilon_p > 0$. 
In implementation, the new time $\tau_{p'}$ is stochastically generated with acceptance probability 1 such that 
%
\begin{equation}
\label{new_unhalted_time}
\tau_{p'} - \tau_p = \tau_{p'p} = - \frac{ \log (1-u) }{ \epsilon_p }
\end{equation}
%
where $u \in [0,1)$ is an uniformly distributed random number.
If $\tau_{p'} > \beta$, then $\tau_{p'}$ will be "modulus"-ed by $\beta$, i.e. $\tau_{p'} \rightarrow \tau_{p'} \, \% \, \beta$.
\item In practice, the value of $\epsilon_p$ can turn out to be enormous, when the system size gets large for instance. When this happens, the algorithm becomes very inefficient. See later section for a clever fix.
\end{enumerate}


\subsubsection{Halted move in time}

\begin{figure}[h]
\centering
\label{extended_worldline_3}
\includegraphics[width=4cm]{./figures/extended-worldlines-3.jpg}
\end{figure}

\noindent Implementing equation \ref{new_unhalted_time}, a new time $\tau_{p'}$ is stochastically generated with chances that it exceeds $\tau_v$ such that $\tau_{p'} > \tau_v$ $\beta$-periodically.
When that happens, the algorithm halts the wormhead just before the vertex at time $\tau_v^-$.
For the halted movement of wormhead over time, the transition probability from extended configuration $X_1$ to $X_2$ must therefore be
%
\begin{equation}
P(X_1 \rightarrow X_2) = \int_{\tau_{vp}}^\infty  \epsilon_p \, e^{-\epsilon_p \tau_{p'p} }  \,\, d  \tau_{p'p}
\end{equation}
%
or
%
\begin{equation}
\label{integrated_exponential_distribution}
P(X_1 \rightarrow X_2) = e^{-\epsilon_p \tau_{vp} } 
\end{equation}
%
\begin{enumerate}
\item By halted move over time itself, it does not satisfy detailed balance. To balance it, the halted move must be immediately followed by an attempt to delete or relink the vertex at time $\tau_v$. (See later section.)
\item Halting of this kind only happens when the wormhead collides with a vertex of conjugate type, ie. either a $\hat{b}^\dag$ wormhead colliding with a $\hat{b}$-vertex (in this example), or a $\hat{b}$ wormhead colliding with a $\hat{b}^\dag$-vertex.
\item There are occasions whereby the wormhead collides with a vertex of the same type. See the following figure.
\begin{figure}[h]
\centering
\label{extended_worldline_4}
\includegraphics[width=4cm]{./figures/extended-worldlines-4.jpg}
\end{figure}
When this happens, the algorithm first halts the wormhead just before the vertex at time $\tau_v^-$. Next, the wormhead crosses the vertex with new time $\tau_v^+$ with probability 1, which is a consequence of the commutation relations $[\hat{b},\hat{b}] = [\hat{b}^\dag,\hat{b}^\dag] = 0$. Finally, a new time is stochastically generated again according to equation \ref{new_unhalted_time}. Its validity will be justified in a later section.
\item An useful property of Poisson events is its {\it{memorylessness}}, ie. the event whereby the wormhead moves unhaltedly to a new time $\tau_{p'}$ is stochastically equivalent to the event whereby the wormhead first moves haltedly to an intermediate time $\tau_v$, then moves unhaltedly to $\tau_{p'}$. One can check easily from equations \ref{exponential_distribution} and \ref{integrated_exponential_distribution}. 
\end{enumerate}
%

\subsubsection{Inserting, deleting, and relinking vertex}  

The figure below is a class of extended configurations that illustrate vertex insertion, deletion, and relinking for 1D systems.
(2D and 3D are natural generalisations.)  
Meanwhile, note that the time interval shown on the figure from $\tau_v$ to $\tau_v^+$ is indeed over-exaggerated only for the purpose of clarity.
In fact, being at time $\tau_v^+$, the wormhead is just infinitesimally above the vertex at time $\tau_v$ for extended configurations $Y_1$, $Y_2$.
%
\begin{figure}[h]
\centering
\label{extended_worldline_5}
\includegraphics[width=12cm]{./figures/extended-worldlines-5.jpg}
\end{figure}
%
The following are the only 2 sub-moves in the algorithm that involves an attempt to jump wormhead across sites.
%
\begin{enumerate}
\item Starting from $X^{(f)}$,
\begin{enumerate}
\item the wormhead could jump to either $Y_1^{(f)}$ or $Y_2^{(f)}$, and inserts a vertex along the way;
\item the wormhead could bounce, or turn around, to $X^{(b)}$.
\end{enumerate}
\item Starting from $Y_1^{(b)}$,
\begin{enumerate}
\item the wormhead could jump to $X^{(b)}$, and deletes the vertex along the way;
\item the wormhead could jump to $Y_2^{(f)}$,  and relinks the vertex along the way;
\item the wormhead could bounce, or turn around, to $Y_1^{(f)}$.
\end{enumerate}
\end{enumerate}
%
The algorithm fixes the following transition probabilities.
%
\begin{eqnarray}
P( X , Y_1 , Y_2 \rightarrow X ) &=& \epsilon_{v^-} \,  \langle i_{v^+}  | \hat{b}_j^\dag |  i_{v^-}  \rangle  \\
P( X , Y_1 , Y_2 \rightarrow Y_1 ) &=&  t \, \langle i_{v^+} | \hat{b}^\dag_i | i_v \rangle  \langle i_v | \hat{b}_i \hat{b}_j^\dag | i_{v^-} \rangle       \\
P( X , Y_1 , Y_2 \rightarrow Y_2 ) &=&  t \, \langle i_{v^+} | \hat{b}^\dag_k | i_v \rangle  \langle i_v | \hat{b}_k \hat{b}_j^\dag | i_{v^-} \rangle                    
\end{eqnarray}
%
To choose which of the sub-moves, one has either the choice of heatbath algorithm or Metropolis algorithm.
It is the preference of the author to choose the latter.
This completes the discussion for all sub-moves of the algorithm in the extended configuration space.


\subsubsection{Crossing vertex}

\begin{figure}[h]
\centering
\label{extended_worldline_6}
\includegraphics[width=4cm]{./figures/extended-worldlines-6.jpg}
\end{figure}

\noindent For the above extended configurations, the relevant weights are
%
\begin{eqnarray}
Z(X) &=& e^{-\epsilon_p \tau_p} \langle i_v | \hat{b}^\dag | i_p \rangle e^{\epsilon_v \tau_p }  \times e^{-\epsilon_v \tau_v} \langle i'_v | \hat{b}^\dag | i_v \rangle e^{\epsilon'_v \tau_v }  \nonumber \\
Z(Y) &=& e^{-\epsilon_p \tau_v} \langle i_v | \hat{b}^\dag | i_p \rangle e^{\epsilon_v \tau_v }  \times e^{-\epsilon_v \tau_v^+} \langle i'_v | \hat{b}^\dag | i_v \rangle e^{\epsilon'_v \tau_v^+ }  \nonumber
\end{eqnarray}
%
or
%
\begin{equation}
\frac{Z(Y)}{Z(X)} = \frac{ e^{- \epsilon_p \tau_{vp} } }{ e^{-\epsilon_v \tau_{vp}} } \nonumber
\end{equation}
%
Balancing with $Z(X) P( X \rightarrow Y ) = Z(Y) P( Y \rightarrow X )$, the algorithm sets
%
\begin{eqnarray}
P(X \rightarrow Y) &=& e^{- \epsilon_p \tau_{vp} }    \nonumber \\
P(Y \rightarrow X) &=& e^{- \epsilon_v \tau_{vp} }    \nonumber 
\end{eqnarray}
%
The above transition probability suggests that whenever the wormhead collides with a vertex of the same type at time $\tau_v$, it first halts at time $\tau_{v}^-$ and then crosses the vertex with probability 1 to a new time $\tau_v^+$. This global move is balanced.


\subsubsection{Summary of moves in extended configuration space}

The following is the flow which summarises all balanced moves in the extended configuration space. 
Without loss of generality, suppose that the wormhead is $\hat{b}^\dag$ and moves forward at time $\tau_p$.

\begin{enumerate}
\item Propose a new time $\tau_{p'}$ for the wormhead according to equation \ref{new_unhalted_time}. 
\item If the wormhead is not halted by any upcoming vertex during the move, then
\begin{enumerate}
\item either the wormhead jumps to a neighbouring site at time $\tau_v^+$ (moving in the same direction) and creates a vertex along the way,
\item or the wormhead turns around in direction (i.e.. bounces).
\end{enumerate}
\item If the wormhead is halted by an upcoming vertex at time $\tau_v$ during the move, then it is halted just before the vertex at time $\tau_v^-$. Either one of the following will occur.
\begin{enumerate}
\item If the upcoming vertex is the wormtail, the wormpair (both the wormhead and wormtail) will be removed. (See next section.) 
\item If the upcoming vertex is of the same type ($\hat{b}^\dag$), the wormhead crosses over the vertex to a new time $\tau_v^+$ with probability 1.
\item If the upcoming vertex is of conjugate type ($\hat{b}$), the wormhead 
\begin{enumerate}
\item either the wormhead jumps to the neighbouring site at time $\tau_v^+$ (moving in the same direction), and deletes the vertex along the way,
\item or the wormhead jumps to the neighbouring site of its neighbouring site at time $\tau_v^-$ (moving in opposite direction), and relinks the vertex along the way,
\item or the wormhead turns around in direction (i.e. bounces).
\end{enumerate}
\end{enumerate}
\item If the wormpair still exists, repeat step 1. Note that if the wormhead now moves backward in time, then the time signs all change in the above.
\end{enumerate}

\subsubsection{Rescaling the diagonal energy}

\begin{figure}[h]
\centering
\label{extended_worldline_7}
\includegraphics[width=1.5cm]{./figures/extended-worldlines-7.jpg}
\end{figure}

\noindent For the algorithm to be more efficient, the diagonal energy should be rescaled. 
A smart choice is
%
\begin{eqnarray}
\epsilon_< &\leftarrow& \epsilon_< - \min \{  \epsilon_< , \epsilon_> \}  + \epsilon_{\mathrm{offset}}  \\
\epsilon_> &\leftarrow& \epsilon_> - \min \{  \epsilon_< , \epsilon_> \}  + \epsilon_{\mathrm{offset}}
\end{eqnarray}
%
where the offset energy $\epsilon_{\mathrm{offset}}$ can be any choice of positive number, i.e. $\epsilon_{\mathrm{offset}}>0$.
Advanced users tune $\epsilon_{\mathrm{offset}}$ for better algorithmic efficiency.
This choice ensures that $\epsilon_< , \epsilon_> > 0$.


\subsubsection{Updating density matrix}

Whenever the wormhead touches the imaginary plane crossing the time of the wormtail, the worms are at equal time. 
When this happens, the corresponding density matrix element is updated. 



\subsection{Updates in configuration space}

The objective of the directed worm algorithm is to generate new configurations via the extended configuration space stochastically.
A configuration becomes an extended configuration upon the creation of a wormpair, and the vice versa upon its annihilation.

\subsubsection{Wormpair creation/ annihilation}

\begin{figure}[h]
\centering
\label{extended_worldline_8}
\includegraphics[width=3cm]{./figures/extended-worldlines-8.jpg}
\end{figure}

\noindent To insert randomly a wormpair, at any time $\tau \in [0,\beta)$, and at any of the L sites, there are $4\beta L$ ways of doing so to any existing configuration, say X.
The factor 4 comes from the fact that, first, either a $\hat{b}-\hat{b}^\dag$ or $\hat{b}^\dag-\hat{b}$ wormpair is created, and second, either the $\hat{b}$ or $\hat{b}^\dag$ worm is the wormhead, thus the other the wormtail.
Therefore, any configuration must be $4\beta L$-degenerate in the extended configuration space, and the relevant weights become
%
\begin{eqnarray}
Z(X) &=& 4\beta L \\
Z(Y) &=& n
\end{eqnarray}
%
where n is the number of particles enclosed within the wormpair.
Since $Z(X) P(X \rightarrow Y) = Z(Y) P(Y \rightarrow X)$, the algorithm defines the transition probabilities as
%
\begin{eqnarray}
P(X \rightarrow Y) &=&  \frac{1}{4\beta L} \\
P(Y \rightarrow X) &=&  \frac{1}{n} \,\, .
\end{eqnarray}
%
\begin{enumerate}
\item In the wormpair-creation process, the probability of locating a random time and a random site is $\frac{1}{4\beta L}$. Therefore, the wormpair is always created with probability 1. 
\item In the wormpair-annihilation process, the probability of choosing any 1 of the n particles to be annihilated is $\frac{1}{n}$. Therefore, the wormpair is always removed with probability 1.
\end{enumerate}
%

\subsubsection{Summary of moves in configuration space}

\begin{enumerate}
\item Starting from any configuration, randomly pick a time and a site.
\item At that location, either insert a $\hat{b} - \hat{b}^\dag$ or a  $\hat{b}^\dag - \hat{b}$ wormpair with 50-50 chance. The following exceptions may arise.
\begin{enumerate}
\item If there is no particle at that location, then insert a $\hat{b}^\dag - \hat{b}$ wormpair with 50\% chance. The other 50\% chance goes to not inserting any wormpair.
\item If the number of particles at that location has reached a maximum, then insert a $\hat{b} - \hat{b}^\dag$ wormpair with 50\% chance. The other 50\% chance goes to not inserting any wormpair.
\end{enumerate}
\item Pick with 50-50 chance either $\hat{b}$ or $\hat{b}^\dag$ as the wormhead, and thus the other being the wormtail.
\item The wormhead random-walks in the extended configuration space, until the wormhead collides into the wormtail.
\item Remove the wormpair.
\item Repeat step 1.
\end{enumerate}


\subsubsection{Updating density matrix}

\begin{enumerate}
\item Just before the annihilation of the wormpair, the diagonal density matrix element is updated.
\item Should any insertion of wormpair fail, the diagonal density matrix element is updated.
\end{enumerate}



