There are many different QMC methods; common feature: Monte Carlo evaluation of
high-dimensional integrals relevant for quantum systems.

\section{Classification of QMC methods}
\begin{enumerate}
 \item Systems: \begin{itemize}
 \item continuum systems (in real space)
 \item lattice systems (in particular: electrons)
 \item single impurity model (relevant for DMFT)
 \item quantum Heisenberg model \ldots
\end{itemize}
\item regime: \begin{itemize}
 \item finite temperatures (canonical ensemble)
 \item ground state (projective QMC)
 \item variational MC (upper bound on $E_{GS}$
\end{itemize}
\item approach: \begin{itemize}
 \item wave function
\item density matrix
 \item auxiliary field
\end{itemize}
\item statistics: \begin{itemize}
 \item "`boltzmannons"' (distinguishable particles)
 \item bosons
 \item fermions
\end{itemize}
\end{enumerate}


\section{Path Integral Monte Carlo (PIMC)}

\begin{swort}{Characterization}
\begin{multicols}{2}
\includegraphics[width=60mm]{bilder/3_1.png}\\
 \begin{itemize}
 \item continuum system
 \item finite temperature
 \item based on density matrix
 \item works best for boltzmannons, good for bosons, fermions require
fixed-node approximation
\item employs analogy between path integrals and polymers
\item beads connected by springs
 \item interaction only between beads of different ring polymers at "equal
time"' (equal color)
\end{itemize}
\end{multicols}
\end{swort}
\begin{swort}{Foundation}
 quantum statistical mechanics
\end{swort}
In terms of \hl{eigen}states $\ket{\Phi_i}$, thermal expectation values of
observables may be written as
\begin{equation}
 \langle\hat{O}\rangle=Z^{-1}\sum_i\langle\Phi_i\vert\hat{O}\vert\Phi_i\rangle
e^{-\beta E_i}
\end{equation}
where $\beta=\frac{1}{k_BT}$, the partition function is 
\begin{equation}
 Z=\sum_i e^{-\beta E_i}
\end{equation}
and
\begin{equation}
 E_i=\langle\Phi_i\vert\hat{\hami}\vert\Phi_i\rangle
\end{equation}
\\
Compare classical case:
\begin{equation}
\begin{split}
 Z_{class}&\propto\int d^{3N}r\int d^{3N}v e^{-\beta E\big(\lbrace r\rbrace,
\lbrace v\rbrace\big)}\\
&E=E_{pot}+E_{kin};\quad E_{kin}=\sum_{i=1}^{3N}\frac{1}{2}m_iv_i^2;\quad
E_{pot}=E_{pot}(\lbrace r\rbrace)\\
Z_{class}&\propto\bigg(\int d^{3N}r e^{-\beta E_{pot}(\lbrace r\rbrace)}\bigg)
\bigg(\int dv e^{-\beta\tfrac{1}{2}mv^2}\bigg)^{3N}\\
&=Z_{pot}\big(\frac{2\pi}{\beta m}\big)^{\frac{3N}{2}}
\end{split}
\end{equation}
Kinetic part has trivial solution $\rightarrow$ MC for static problem\\
\\
For an arbitrary basis (in which the Hamiltonian is in general not diagonal):
\begin{equation}
\begin{split}
 \langle\hat{O}\rangle&=\frac{1}{Z}\sum_{\alpha}\langle\alpha\vert\hat{O}
e^{-\beta\hat{\hami}}\vert\alpha\rangle\color{darkgreen}=\frac{1}{Z}
\operatorname{tr}(\hat{O}e^{-\beta\hat{\hami}})\\
&=\frac{1}{Z}\sum_{alpha}\sum_{\alpha'}\langle\alpha\vert\hat{O}\vert\alpha'
\rangle\langle\alpha'\vert e^{-\beta\hat{\hami}}\vert\alpha\rangle
\end{split}
\end{equation}
with $\displaystyle Z=\sum_{\alpha}\langle\alpha\vert e^{-\beta\hat{\hami}}
\vert\alpha\rangle=\operatorname{tr}(e^{-\beta\hat{\hami}})$\\
Central operator: {\color{red}\fbox{\hl{density matrix} $\color{black}\displaystyle e^{-\beta\hat{\hami}}$}}
\\ \\
Now choose \hl{position basis with labeled particles}.\\
Advantage: in this basis, all matrix elements are \hl{real} and \hl{positive}
(proof: later)\\
\hl{Position-space density matrix:}
\redbox{
 \rho(R,R',\beta)=\langle R\vert e^{-\beta\hat{\hami}}\vert R'\rangle=
\sum_i\Phi_i*(R)\Phi_i(R')e^{-\beta\hat{\hami}}
}
with position representations of the eigenstates
\begin{equation}
 \Phi_i(R')=\langle i\vert R'\rangle;\quad \Phi_i*(R)=\langle R\vert i=
\langle i\vert R\rangle*
\end{equation}
and position vectors $R=\lbrace\vec{r}_1,\ldots,\vec{r}_N\rbrace$.\\
In $d=3$ dimensions $\rho$ is a function of $6N+1$ variables. In position basis
, the expectation values become:
\begin{equation}
 \langle\hat{O}\rangle=Z^{-1}\int dRdR'\:\rho(R,R',\beta)\langle R'\vert
\hat{O}\vert R\rangle
\end{equation}
with $Z=\int dR\:\rho(R,R',\beta)$
\\ \\
\begin{swort}{Note}
 the product of two density matrices is a density matrix (all for the same
Hamiltonian!):
\begin{equation}
 e^{-\beta_1\hat{\hami}}e^{-\beta_2\hat{\hami}}=e^{-(\beta_1+\beta_2)
\hat{\hami}}
\end{equation}
\begin{itemize}
 \item in this case, operators can be treated like numbers since $\lbrack
\hat{\hami},\hat{\hami}\rbrack=0$
\item the product matrix is associated with a \hl{lower temperature}
\end{itemize}
\end{swort}
Written for  positions $\rightarrow$ \hl{convolution}
\redbox{
 \rho(R_1,R_3;\beta_1+\beta_2)=\int dR_2\rho(R_1,R_2;\beta_1)\rho(R_2,R_3;
\beta_2)
}


\subsection{Discrete path integrals}

Now consider $\hat{\hami}=\hat{T}+\hat{V}$ with the usual kinetic contribution
(nonrelativistic, no magnetic field):
\begin{equation}
 \hat{T}=\sum_{n=1}^N-\frac{1}{2m}\hat{p}^2_n\qquad(\hat{p}_n=-i\hbar
\vec{\nabla}_n\;\textrm{in position basis})
\end{equation}
and the potential contribution:
\begin{equation}
 \hat{V}=\int dR\:\underbrace{V(R)}_{\color{darkgreen}=V(\hat{R})}\ket{R}
\bra{R}
\end{equation}
The fundamental problem in quantum mechanics is that we know eigenbasis and
density matrix only for either the kinetic \hl{or} the potential part of the
Hamiltonian, not both:
\begin{equation}
\sideset{}{^{-\beta\hat{\hami}}}{
\operatorname*{e}}_{\substack{\textrm{needed,}\\ \textrm{unknown}}} \neq 
\quad\sideset{}{^{-\beta\hat{T}}}{
\operatorname*{e}}_{\substack{\textrm{known}\\ \textrm{(see below)}}}
\qquad\sideset{}{^{-\beta\hat{V}}}{
\operatorname*{e}}_{\substack{\textrm{known}\\ \textrm{(see below)}}}
\qquad(\textrm{for general }\beta)\\
\end{equation}
\begin{swort}{But}
 asymptotic factorization at high temperatures (classical limit),i.e. small
$\beta$. More generally: good high-T approximations for density matrix!
\end{swort}\\ \\
\begin{swort}{Notation}
 imaginary time step $\Delta\tau=\frac{\beta}{M}\quad(M\in\mathbb{N})$
\end{swort}
Exact discrete path integral ($M$ equal steps):
\redbox{
 \rho(R_0,R_M,\beta)&=&\int dR_1\int dR_2\ldots\int 
dR_{M-1}\rho(R_0,R_1;\Delta\tau)\\ &&\rho(R_1,R_2;\Delta\tau)
\ldots\rho(R_{M-1},R_M;\Delta\tau)
}
We will now derive an explicit approximation for this path integral.
\\
Campbell-Baker-Haussdorf formula:
\begin{equation}
 \exp\big\lbrack-\Delta\tau(\hat{T}+\hat{V})+\underbrace{\tfrac{1}{2}\Delta
\tau^2\lbrack\hat{T},\hat{V}\rbrack+\mathcal{O}(\Delta\tau^3)}_{
\textrm{negligible in limit of small }\Delta\tau}\big\rbrack=e^{-\Delta\tau
\hat{T}}e^{-\Delta\tau\hat{V}}
\end{equation}
$\longrightarrow$ \hl{Trotter formula:}
\redbox{
e^{-\beta(\hat{T}+\hat{V})}&=&\lim_{M\rightarrow\infty}\big\lbrack
e^{-\Delta\tau\hat{T}}e^{-\Delta\tau\hat{V}}\big\rbrack^M;\quad\Delta\tau=
\tfrac{\beta}{M}\\
&\approx&\big\lbrack e^{-\Delta\tau\hat{T}}e^{-\Delta\tau\hat{V}}\big\rbrack^M
\;\;\textrm{for }\Delta\tau\textrm{ small}
}
\hl{Primitive approximation} in position space for single time step:
\begin{equation}
 \rho(R_0,R_2;\Delta\tau)=\int dR_1\:\langle R_0\vert e^{-\Delta\tau\hat{T}}
\vert R_1\rangle\langle R_1\vert e^{-\Delta\tau\hat{V}}\vert R_2\rangle
\end{equation}
Potential operator is diagonal, therefore:
\begin{equation}
 \langle R_{1}\vert e^{-\Delta\tau\hat{V}}\vert R_2\rangle=e^{-\Delta\tau
V(R_1)}\delta(R_1-R_2)\tag{{\color{red}\textbf{*}}}
\end{equation}
For kinetic energy operator, introduce eigenfunction expansion in a cube with
side length $L$ and pbc, denoted by integer vector $\vec{n}\in\mathbb{R}^{3N}$:
\begin{equation}
 \Phi_{\vec{n}}(\vec{R})=\langle\vec{n}\vert\vec{R}\rangle=L^{-\tfrac{3N}{2}}
e^{-i\vec{k}_{\vec{n}}\cdot\vec{R}};\quad\vec{k}_{\vec{n}}=\frac{2\pi\vec{n}}
{L}\quad\textrm{plane waves}
\end{equation}
with eigen values $\displaystyle\frac{\hbar^2}{2m}\vec{k}_{\vec{n}}^2\equiv
\lambda\vec{k}_{\vec{n}}^2$
\begin{equation}
\begin{split}
 \langle R_0\vert e^{-\Delta\tau\hat{T}}\vert R_1\rangle&=\sum_{\vec{n}}
\langle R_0\vert\vec{n}\rangle\langle\vec{n}\vert e^{-\Delta\tau\hat{T}}\vert
\vec{n}\rangle\langle\vec{n}\vert R_1\rangle\\
&=\sum_{\vec{n}}L^{-\tfrac{3N}{2}}e^{-i\vec{k}_{\vec{n}}\cdot\vec{R}_0}
e^{-\Delta\tau\lambda\vec{k}_{\vec{n}}^2}L^{-\tfrac{3N}{2}}
e^{+i\vec{k}_{n}\cdot\vec{R}_1}\\
&=\sum_{\vec{n}}L^{-3N}\exp\big\lbrack-\Delta\tau\lambda\vec{k}_{\vec{n}}^2
-i\vec{k}_{\vec{n}}\cdot(\vec{R}_0-\vec{R}_1)\big\rbrack\\
&\approx(4\pi\lambda\Delta\tau)^{-\tfrac{3N}{2}}\exp\bigg\lbrack
-\frac{(\vec{R}_0-\vec{R}_1)^2}{4\lambda\Delta\tau}\bigg\rbrack\\
&\uparrow\\
&\color{darkgreen}\sum_{\vec{n}}\rightarrow\int d^{3n}n,\;\;\textrm{good for }
\lambda\Delta\tau\ll L^2\textrm{ then Gauss integral}
\end{split}\tag{{\color{red}$\Delta$}}
\end{equation}
(exact kinetic density matrix: $\Theta_3(z,q)\ldots$)
\\
Insert (lowest order) Trotter approximation, using $\color{red}\mathbf{*}$ and
$\color{red}\Delta$, into general discrete path integral:
\redbox{
\rho(\vec{R}_0,\vec{R}_M;\beta)&=&\int d\vec{R}_1\:\int d\vec{R}_2\:\ldots\:
\int d\vec{R}_{M-1}(4\pi\lambda\Delta\tau)^{-\frac{3NM}{2}}\\
&&\exp\Big\lbrack-\sum_{m=1}^M\big(\frac{(\vec{R}_{m-1}-\vec{R}_m)^2}{4\lambda
\Delta\tau}+\Delta\tau V(\vec{R}_{m})\big)\bigg\rbrack\color{darkgreen}\geq0\\
&&\color{red}\textrm{primitive approximation for density matrix}
}

\begin{multicols}{2}
 \includegraphics[width=60mm]{bilder/3_2.png}\columnbreak\\
$Z=\int d\vec{R}\:\rho(\vec{R},\vec{R};\beta)$ corresponds to classical
partition function of ring polymers with ideal springs within each polymer and
interpolymer interaction only at equal times ($\tilde{=}$ color).
\end{multicols}

More explicitly: position-dependent part of $Z$ for pair interactions
\begin{equation}
\begin{split}
 &V(\vec{R}_m)=\sum_{i=1}^N\sum_{j=i+1}^Nv\big(\vec{r}_i(\tau_m)-\vec{r}_j
(\tau_m)\big)\\
\textrm{is}
&\int d\vec{R}_1\:\int d\vec{R}_2\:\ldots\:\int d\vec{R}_Me^{-\beta\tilde{V}
(\lbrace\vec{R}\rbrace)}
\end{split}
\end{equation}
with
\redbox{
\tilde{V}(\lbrace\vec{R}\rbrace)&=&\sum_{i=1}^N\sum_{j=i+1}^N\frac{1}{M}
\sum_{m=1}^Mv\big(\vec{r}_i(\tau_m)-\vec{r}_j(\tau_m)\big)\\
&+&\sum_{i=1}^N{\color{darkgreen}\underbrace{\color{black}
\frac{M}{4\lambda\beta^2}}_{\tfrac{k}{2}}}\sum_{m=1}^N\big(
\vec{r}_i(\tau_{m-1})-\vec{r}_j(\tau_m)\big)^2
}

\subsection{Fractal nature}

In the case of a smooth imaginary-time path $\vec{r}(\tau)$ we would expect
$\vert\vec{r}(\tau_i)-\vec{r}(\tau_j)\vert^2\leq v_{max}^2\vert\tau_i-\tau_j
\vert^2$ however, in limit $V\rightarrow0$ we get
\begin{equation}
 \langle\vert\vec{r}_i(\tau_m)-\vec{r}_i(\tau_{m-1})\vert^2\rangle=
3(2\lambda\Delta\tau)
\end{equation}
$\rightarrow$ in the limit $\Delta\tau\rightarrow0$, each polymer acquires
infinite length! (smoothing + shading)

\subsection{Role of temperature}

In the polymer model, the finite size of each polymer is an entropic effect:
for fixed spring constant $k$, the size would vanish for $T\rightarrow0$.
\begin{tabbing}
 However: \=for fixed $M$; $k\propto T^2\Rightarrow\vert\vec{r}_i(\tau_m)-
\vec{r}_i(\tau_{m-1})\vert^2\propto T^{-1}\rightarrow\infty$\\
\>for fixed $\tau$: $k\propto T\Rightarrow\vert\vec{r}_i(\tau_m)-\vec{r}_i(
\tau_{m-1})\vert^2\propto T^0\rightarrow\textrm{const}$ and $M\rightarrow
\infty$
\end{tabbing}
in both cases: polymer diameter $\propto\sqrt{\lambda\beta}\rightarrow\infty$

\subsection{Connection to Feynman-Kac path integral}

For simplicity, write only for $Z$ (analogous for $\rho$):
\begin{equation}
\begin{split}
 Z&=\bigg\langle(4\pi\lambda\Delta\tau)^{-\tfrac{3NM}{2}}\prod_{m=1}^M
e^{-\frac{\vert\vec{R}_{m-1}-\vec{R}_m\vert^2}{4\lambda\Delta\tau}}e^{
-\Delta\tau V(\vec{R}_m)}\bigg\rangle_{\textrm{uniform}}\\
&=Z_{\textrm{free}}\bigg\langle\prod_{m=1}^M e^{-\Delta\tau V(\vec{R}_m)}\bigg
\rangle_{\textrm{free}}\quad{\color{darkgreen}\lbrace R_m\rbrace
\textrm{ generated according to $Z$ for }V=0}\\
&=Z_{\textrm{free}}\Big\langle\exp\big\lbrack-\sum_{m=1}^M\Delta\tau V(\vec{R}
_m)\big\rbrack\Big\rangle_{\textrm{free}}\\
&\xrightarrow{M\rightarrow\infty}Z_{\textrm{free}}\langle\exp\lbrack-\int_0
^\beta V\big(\vec{R}(\tau)\big)\:d\tau\rbrack\rangle_{\textrm{free}}
\end{split}
\end{equation}
\redbox{
e^{-\beta(F-F_{\textrm{free}})}=\frac{Z}{Z_\textrm{free}}=
\langle\exp\lbrack-\int_0^\beta V\big(\vec{R}(\tau)\big)\:d\tau\rbrack\rangle
_{\textrm{free}}
}

\subsection{Monte Carlo sampling of the partition function}

Not possible: naive simple sampling (take all $\frac{r_i}{\tau_m}$ as
independent random variables in $L^3$)
\begin{enumerate}
 \item single particle displacement moves ($\sim$ classical atoms)
\begin{enumerate}
 \item select particle $i$
 \item draw displacement $\vec{\delta}$from uniform distribution inside
cube $\lbrack-\frac{\Delta}{2},\frac{\Delta}{2}\rbrack^3$
 \item propose move $\vec{r}_i(\tau_m)\longrightarrow\vec{r}_i(\tau_m)
+\vec{\delta}$ for all m
 \item accept with probability $p=\min\lbrace1,e^{-\beta\Delta V}\rbrace$
\end{enumerate}
\begin{itemize}
 \item These moves are only usefull for not too low $T$
\item $\Delta$ is adjusted for acceptance probability of $25\%-75\%$
\end{itemize}
\item single-particle single-slice moves
\begin{enumerate}
 \item[]\begin{enumerate}\item select particle $i$ and slice $m$\end{enumerate}
\end{enumerate}
\begin{enumerate}
\item classic rule
\begin{enumerate}
 \item[ii.] random displacement $\vec{\delta}$ uniformly from cube
$\lbrack-\frac{\Delta}{2},\frac{\Delta}{2}\rbrack^3$
\item[iii.] propose move $\vec{r}_i(\tau_m)\longrightarrow\vec{r}_i(\tau_m)
+\vec{\delta}$ for single $m$
\item[iv.] accept with probability $p=\min\lbrace1,e^{-\beta\Delta\tilde{V}}\rbrace$
\end{enumerate}
\begin{itemize}
 \item now also kinetic contributions to $\Delta\tilde{V}$
\item adjust $\Delta\approx\sqrt{\lambda\Delta\tau}$ to $\sim50\%$ acceptance
prob.
\end{itemize}
\item free-particle sampling\\
In general, the probability density for position vector $\vec{R}_m$ at fixed
$\vec{R}_i$ for $i\neq m$ can be written as
\begin{equation}
 P\lbrack\vec{R}_m\rbrack\propto\rho(\vec{R}_{m-1},\vec{R}_m;\Delta\tau)\rho(
\vec{R}_m,\vec{R}_{m+1};\Delta\tau)
\end{equation}
specifically for $\hat{V}=0$ we have
\begin{equation}
\begin{split}
 P\lbrack\vec{R}_m\rbrack&\propto\exp\bigg\lbrack-\frac{\vert\vec{R}_{m-1}-
\vec{R}_m\vert^2}{4\lambda\Delta\tau}\bigg\rbrack\exp\bigg\lbrack-\frac{
\vert\vec{R}_m-\vec{R}_{m+1}\vert^2}{4\lambda\Delta\tau}\bigg\rbrack\\
\color{darkgreen}(x-a)^2+(x-b)^2&\color{darkgreen}=2x^2-2x(a+b)+a^2+b^2\\
&\color{darkgreen}=2\Big\lbrack x^2-2x\frac{a+b}{2}+\Big(\frac{a+b}{2}\Big)^2
\Big\rbrack+\frac{(a-b)^2}{2}\\
\Rightarrow P\lbrack\vec{R}_m\rbrack&\propto\exp\bigg\lbrack-\frac{\vert
\vec{R}_m-\vec{R}_{av}\vert^2}{2\lambda\Delta\tau}\bigg\rbrack;\quad
\vec{R}_{av}=\frac{\vec{R}_{m-1}+\vec{R}_{m+1}}{2}
\end{split}
\end{equation}
in case of single particle move:
\begin{equation}
 P\lbrack r_i(\tau_m)\rbrack\propto\exp\Big\lbrack-\frac{\vert\vec{r}_i(\tau_m)
-\vec{r}_i^{av}\vert^2}{2\lambda\Delta\tau}\Big\rbrack;\quad\vec{r}_i^{av}=
\frac{\vec{r}_i(\tau_{m+1})+\vec{r}_i(\tau_{m-1})}{2}
\end{equation}
This can be sampled directly (diffusion step)
\begin{enumerate}
 \item[ii.] draw $\vec{\delta}$ from Gaussian distribution with width $\sqrt{\lambda
\Delta\tau}$ (in each dimension)
\item[iii.] propose move $\vec{r}_i(\tau_m)\rightarrow\frac{\vec{r}_i(\tau_{m+1})+
\vec{r}_i(\tau_{m-1})}{2}+\vec{\delta}$
\item[iv.] accept with probability $p=\min\lbrace1,e^{-\beta\Delta V}\rbrace$\\
(this expression is valid for primitive approximation; generalize for better
actions)
\end{enumerate}
procedure analogous to diffusion Monte Carlo
\begin{equation}
 \ket{O}=\lim_{\beta\rightarrow\infty}e^{-\beta\hat{\hami}}\ket{\Psi_0}\approx
\prod_{m=1}^M\:\;\big(
\sideset{}{^{-\Delta\tau\hat{T}}}
{\operatorname{e}}_{\substack{
\uparrow\\ \textrm{diffusion}\\ \textrm{of "`walkers"'}}}
\quad\;
\sideset{}{^{-\Delta\tau\hat{V}}}
{\operatorname{e}}_{\substack{
\uparrow\\ \textrm{birth/death}\\ \textrm{of "`walkers"'}}}
\big)\;\ket{\Psi_0}
\end{equation}
\end{enumerate}
\end{enumerate}
1) and 2) are in principle sufficient for distinguishable particles. However:
single-slice moces inefficient for large $M$
\begin{enumerate}
 \item[2'.] Multilevel sampling\\
% \begin{figure}[hpt]
%  \centering
\includegraphics[width=120mm]{bilder/3_3.png}
% \end{figure}

On fine grid: $\vec{r}_4$ can not move very far, idea:
\begin{itemize}
 \item consider coarse grid with $\tilde{\Delta\tau}=4\Delta\tau$ first
\item construct move iteratively
\item keep track of intermediate proposal prob.
\item correct everything in final acceptance rule
\end{itemize}
\item[3.] Exchange moves: for distinguishable particles, we must allow
exchanges (discrete move - not possible in MD)\\
\includegraphics[width=120mm]{bilder/3_4.png}
\\
In practice: combine with multilevel moves.
\end{enumerate}

\subsection{Observables}

\begin{itemize}
 \item density, pair correlation function, structure factor: directly from
$\vec{R}_m$
\item already the energy is nontrivial
\begin{enumerate}
 \item "`direct"' or "`Hamiltonian"' estimator:
\begin{equation}
 E_{\hami}=\frac{\operatorname{Tr}(\hat{\hami}e^{-\beta\hat{\hami}})}{
\operatorname{Tr}(e^{-\beta\hat{\hami}}}=\ldots
\end{equation}
\item thermodynamic estimator:
\begin{equation}
\begin{split}
 E_T&=-\frac{1}{Z}\frac{dZ}{d\beta}\\
&=\frac{1}{M}\sum_{m=1}^M\big\lbrace\frac{3N}{2\Delta\tau}-\frac{\vert\vec{R}_m
-\vec{R}_{m-1}\vert^2}{4\lambda\Delta\tau^2}+V(\vec{R}_m)\big\rangle
\end{split}
\end{equation}
\item virial estimator:
\begin{equation}
 E_V=\langle\ldots-\frac{1}{2}\hspace{-15pt}
\sideset{}{^m}{
\operatorname{F}}_{\substack{\color{darkgreen}\uparrow\\ \color{darkgreen} 
\textrm{generalized force}}}\hspace{-35pt}
\sideset{}{_m}{\operatorname{\Delta}}^{\substack{\color{darkgreen}
\textrm{derivaion from}\\ \color{darkgreen}\textrm{center of mass}\\
\color{darkgreen}\downarrow}}\hspace{-15pt}
\ldots\rangle
\end{equation}
\end{enumerate}
\item condensate fraction\\
insert two free ends $\vec{r}_0,\vec{r}_+$\\
measure prob$(\vec{r}_0-\vec{r}_+)$\\
no exchange: \\
{
 \centering
\includegraphics[width=100mm]{bilder/3_5.png}
}\\
with exchange: \\
{\centering
\includegraphics[width=60mm]{bilder/3_6.png}
\includegraphics[width=50mm]{bilder/3_7.png}}
\item superfluid density\\
\includegraphics[width=110mm]{bilder/3_8.png}
\end{itemize}

\subsection{Applicaion to vacancies in solid ${}^4$He}

see slides: LINK

\section{Weltlinien-Quanten-Monte-Carlo-Methode (world line QMC)}
\label{sec:3.3}

Weltlinien-QMC ist ein konzeptionell relativ einfacher Zugang für Modelle
wechselwirkender Quanten-Spins oder Fermionen, der erst Simulationen erlaubt
und Basis für effizientere Algorithmen (Schleifen-, Wurm-Algorithmen) ist.\\
Als Anwendungsbeispiel wählen wir das XXZ-Modell:
\begin{equation}
 \hami=\sum_{i=1}^L\hami^{(i)};\;\hami^{(i)}=\frac{J_x}{a}(S_i^+S_{i+1}^-+
S_i^-S_{i+1}^+)+J_z(S_i^z+S_{i+1}^z)
\end{equation}
mit $S_i^x=\frac{1}{2}\left(\begin{smallmatrix}0&1\\1&0\end{smallmatrix}\right),\;
S_i^y=\frac{1}{2}\left(\begin{smallmatrix}0&-i\\i&0\end{smallmatrix}\right),\;
S_i^z=\frac{1}{2}\left(\begin{smallmatrix}1&0\\0&-1\end{smallmatrix}\right)$\\
und periodischen Randbedingungen: $\vec{S}_{L+1}\equiv\vec{S}_1$\\
Basis des Hilbert-Raums sind $\sigma_i\in\big\lbrace\ket{\uparrow}=
\left(\begin{smallmatrix}1\\0\end{smallmatrix}\right);\;\ket{\downarrow}=
\left(\begin{smallmatrix}0\\1
\end{smallmatrix}\right)\big\rbrace$\\
Mit $S^+=S^x+iS^y,\;S^-=S^x-iS^y$ gilt:
\begin{equation}
 \hami=\frac{J_x}{2}\sum_{i=1}^L(S_i^+S_{i+1}^-+S_i^-S_{i+1}^+)+J_z\sum_{i=1}^L
(S_i^zS_{i+1}^z)
\end{equation}
Für 2 Teilchen mit offenen Randbed.:
\begin{equation}
\begin{split}
 &\;\;\;\uparrow\uparrow\;\;\;\;\uparrow\downarrow\;\;\;\;\;\downarrow\uparrow\;
\;\;\;\downarrow\downarrow\\
\hami_{12}\hat{=}&\begin{pmatrix}
\frac{J_z}{4}&0&0&0\\0&-\frac{J_z}{4}&\frac{J_x}{2}&0\\0&\frac{J_x}{2}&-\frac{
J_z}{4}&0\\0&0&0&\frac{J_z}{4}\end{pmatrix}\begin{array}{c}
 \uparrow\uparrow\\ \uparrow\downarrow\\ \downarrow\uparrow\\
\downarrow\downarrow
\end{array}
\end{split}
\end{equation}
Damit ergeben sich für Eigenvektoren und Eigenwerte des 2-Spin-Systems:
\begin{equation}
 \begin{matrix}
\vec{S}^2&S^z&\textrm{EV}&\textrm{EW}\\
1&1&\ket{\uparrow\uparrow}&\frac{J_z}{4}\\
1&-1&\ket{\downarrow\downarrow}&\frac{J_z}{4}\\
1&0&\frac{1}{\sqrt{2}}(\ket{\uparrow\downarrow}+\ket{\downarrow\uparrow})&
\frac{J_x}{2}-\frac{J_z}{4}\\
0&0&\frac{1}{\sqrt{2}}(\ket{\uparrow\downarrow}-\ket{\downarrow\uparrow})&
-\frac{J_x}{2}-\frac{J_z}{4}
 \end{matrix}
\tag*{\color{red}(\textbf{*})}
\end{equation}
Grundidee von Weltlinien-Verfahren ist es nun, den Hamilton-Operator so
aufzuteilen, dass jeder der Teile in kommutierende 2-Platz-Terme zerfällt.
Bei Modellen, in denen $\hami$ maximal nächste Nachbarn auf einem
hyperkubischen Gitter verbindet, sind das 2d Teile:
$\hami=\hami_1^x+\hami_2^x+\hami_1^y+\hami^y_2+\ldots$
\\ \\
\begin{swort}{Wichtig}
 Notwendige Voraussetzung ist, dass die Systemlänge in jeder Richtung gerade
ist: $L^(\alpha)=2n_{\alpha}$
\end{swort}
Speziell im vorliegenden 1-d-Fall:
\begin{equation}
 \hami={\color{darkgreen}\underbrace{\color{black}\sum_{n=1}^{\frac{L}{2}}
\hami^{(2n)}}_{\color{red}\hami_1}}+{\color{darkgreen}\underbrace{
\color{black}\sum_{n=1}^{\frac{L}{2}}
\hami^{(2n)}}_{\hami_2}}
\end{equation}
{
\centering
\includegraphics[width=50mm]{bilder/3_9.png}}
\newpage
Diese Aufteilung wird nun in einer Trotter-Zerlegung der Zustandssumme
benutzt ($\Delta\tau=\frac{\beta}{M}$):
\begin{equation}
\begin{split}
 Z=\operatorname{Sp}\lbrack e^{-\beta\hami}\rbrack&=\operatorname{Sp}\big
\lbrack(e^{-\Delta\tau\hami})^M\big\rbrack\\
&\approx\operatorname{Sp}\big\lbrack(e^{-\Delta\tau\hami_1}e^{-\Delta\tau\hami
_2})^M\big\rbrack\\
&=\sum_{\vec{\sigma}_1\ldots\vec{\sigma}_2}\langle\vec{\sigma}_1\vert
e^{-\beta\hami_1}\vert\vec{\sigma}_{2M}\rangle\langle\vec{\sigma}_{2M}\vert
e^{-\beta\hami_2}\vert\vec{\sigma}_{2M-1}\rangle\\
&\qquad\langle\vec{\sigma}_3\vert e^{-\beta\hami_1}\vert\vec{\sigma}_2\rangle
\langle\vec{\sigma}_2\vert r^{-\beta\hami_2}\vert\vec{\sigma}_1\rangle
\end{split}
\end{equation}
mit dem Gesamtspinvektor $\vec{\sigma}_{\tau}$ zur imaginären Zeit $\tau$:
\begin{equation}
 \vec{\sigma}_{\tau}=(\sigma_{1,\tau},\sigma_{2,\tau},\ldots,\sigma_{L,\sigma})
\end{equation}
Graphische Repräsentation: Weltlinien-Evolution von $\uparrow$-Zuständen
(bei Fermion- Simulationen: von besetzten Zuständen) in der imaginären Zeit:\\
\includegraphics[width=\linewidth]{bilder/3_10.png}
$L=6$ Spins, davon 2 up-Spins ($S^z=-1$)\\
$M=2$ Diskretisierungsschritte der imaginären Zeit.\\
\\
Wir führen noch einige Bezeichnungen ein:\\
Konfiguration: $w=(\vec{\sigma}_1,\vec{\sigma}_2,\ldots,
\vec{\sigma}_{2M})$\\
Gewicht einer Konfiguration:
\begin{equation}
\begin{split}
\Omega(w)=\langle\vec{\sigma}_1\vert e^{-\beta\hami_1}&\vert
\vec{\sigma}_{2M}\rangle\langle\vec{\sigma}_{2M}\vert e^{-\beta\hami_2}\vert
\vec{\sigma}_{2M-1}\rangle\ldots\\
&\langle\vec{\sigma}_3\vert e^{-\beta\hami_1}\vert\vec{\sigma}_2\rangle\langle
\vec{\sigma}_2\vert e^{-\beta\hami_2}\vert\vec{\sigma}_1\rangle
\end{split}
\end{equation}
Damit gilt: $Z=\sum_{w}\Omega(w)$\\
Die Matrixwlemente in $\Omega(w)$ faktorisieren:
\begin{equation}
 \langle\vec{\sigma}_{\tau+1}\vert e^{-\Delta\tau\hami_2}\vert\vec{\sigma}
_{\tau}\rangle=\prod_{i=1}^{\frac{L}{2}}\langle\sigma_{2i,\tau+1},\sigma_{2i+1,
\tau+1}\vert e^{-\Delta\tau\hami^{(2i)}}\vert\sigma_{2i,\sigma},\sigma
_{2i+1,\tau}\rangle
\end{equation}
Damit ist die Berechnung der Gewichte auf das 2-Platz-Problem zurückgeführt.
Z.B. erhalten wir mit $\color{red}*$:
\begin{equation}
\begin{split}
 \langle\downarrow\uparrow\vert e^{-\Delta\tau\hami_2}\vert\uparrow\downarrow
\rangle&=\frac{1}{\sqrt{2}}\langle\downarrow\uparrow\vert e^{-\Delta\tau
\hami_2}\big\lbrack\frac{1}{\sqrt{2}}(\ket{\uparrow\downarrow}-\ket{\downarrow
\uparrow})+\frac{1}{\sqrt{2}}(\ket{\uparrow\downarrow}+\ket{\downarrow\uparrow}
)\big\rbrack\\
&=\frac{1}{\sqrt{2}}\bra{\downarrow\uparrow}e^{-\Delta\tau(-\frac{J_z}{4}
-\frac{J_x}{2})}\frac{1}{\sqrt{2}}(\ket{\uparrow\downarrow}-\ket{\downarrow
\uparrow})\\
&\qquad+e^{-\Delta\tau(-\frac{J_z}{4}+\frac{J_x}{2})}\frac{1}{\sqrt{2}}
(\ket{\uparrow\downarrow}+\ket{\downarrow\uparrow})\big\rbrack\\
&=-e^{\Delta\tau\frac{J_z}{4}}\sinh\Big(\frac{\Delta\tau J_x}{2}\Big)
\end{split}
\end{equation}
Analog erhält man für die Plaketten: \\
{\centering\includegraphics[width=100mm]{bilder/3_11.png}}
% $e^{-\Delta\tau\frac{J_z}{4}}$\\
% $e^{\Delta\tau\frac{J_z}{4}}\cosh\big(\Delta\tau\frac{J_x}{2}\big)$\\
% $-e^{\Delta\tau\frac{J_z}{4}}\sinh\big(\Delta\tau\frac{J_x}{2}\big)$\\
\\
Beachte: negative Gewichte würden i.A. zu einem Vorzeichenproblem führen.
Bei nächst-Nachbar-Hüpfen auf bipartiten Gittern ist die Anzahl der
Diagonalplaketten jedoch in jeder Konfiguration gerade, das Vorzeichen also
irrelevant.\\
\\
Weltlinien-Algorithmus
\begin{enumerate}
 \item initialisiere Weltlinien (z.B. senkrecht) für gewünschtes $L,S_z$
\item Wähle schattierte Plakette, schlage Update vor (möglich, falls schatt.
Plakette AF, stationär) \\
{\centering\includegraphics[width=90mm]{bilder/3_12.png}}
\item Berechne Verhältnis der Gewichte aus 4 umgebenden weissen Plaketten
für $w_{alt},w_{neu}$:
\begin{equation}
 r=\frac{\Omega(w_{neu})}{\Omega(w_{alt})}=\frac{\Big(\fbox{N}\cdot\fbox{O}
\cdot\fbox{S}\cdot\fbox{W}\Big)_{neu}}{\Big(\fbox{N}\cdot\fbox{O}\cdot
\fbox{S}\cdot\fbox{W}\Big)_{alt}}
\end{equation}
Akzeptiere Update mit W' $p=\frac{r}{1+r}$ (heat bath)\\
(Alternative: Metropolis - $p=\min\lbrace r,1\rbrace$)
\end{enumerate}
Berechnung von Observablen $\langle O\rangle=\frac{\operatorname{Sp}\lbrack
e^{-\beta\hami}O\rbrack}{\operatorname{Sp}\lbrack e^{-\beta\hami}\rbrack}$\\
Für Observablen, die lokal in der imaginären Zeit sind, kann man schreiben:
$\langle O\rangle=\frac{\sum_w\Omega(w)O(w)}{\sum_w\Omega(w)}$\\
Die Kunst ist es nun, einen Schätzer $O(w)$ zu bestimmen, der (i) korrekt ist
und (ii) möglichst kleine Varianz hat.

\subsection{Berechnung der Energie}

\begin{equation}
\begin{split}
 \langle\hami\rangle&=\frac{1}{Z}\operatorname{Sp}\big\lbrack(e^{-\Delta\tau
\hami_1}e^{-\Delta\tau\hami_2})^M(\hami_1+\hami_2)\big\rbrack\\
&=\frac{1}{Z}\operatorname{Sp}\big\lbrack(e^{-\Delta\tau
\hami_1}e^{-\Delta\tau\hami_2})^{M-1}(e^{-\Delta\tau
\hami_1}\hami_1e^{-\Delta\tau\hami_2}+e^{-\Delta\tau
\hami_1}e^{-\Delta\tau\hami_2}\hami_2)\big\rbrack
\end{split}
\end{equation}
Durch Einschieben von Einheitsoperatoren $1=\sum_{\sigma}\ket{\vec{\sigma}}
\bra{\vec{\sigma}}$ erhalten wir:
\begin{equation}
\begin{split}
 \langle\hami\rangle&=\frac{1}{Z}\sum_{\vec{\sigma}_1,\ldots,\vec{\sigma}_{2M}}
\Big\rbrack\langle\vec{\sigma}_{1}\vert e^{-\Delta\tau\hami_{1}}\vert
\vec{\sigma}_{2M}\rangle\ldots\langle\vec{\sigma}_{3}\vert e^{-\Delta\tau
\hami_{1}}\vert\vec{\sigma}_{2}\rangle\langle\vec{\sigma}_{2}\vert e^{-\Delta
\tau\hami_{2}}\hami_2\vert\vec{\sigma}_{2}\rangle\\
&\qquad+\langle\vec{\sigma}_{1}\vert e^{-\Delta\tau\hami_{1}}\vert\vec{\sigma}
_{2M}\rangle\ldots\langle\vec{\sigma}_{3}\vert e^{-\Delta\tau\hami_{1}}\hami_1
\vert\vec{\sigma}_{2}\rangle\langle\vec{\sigma}_{2}\vert e^{-\Delta\tau
\hami_{2}}\vert\vec{\sigma}_{1}\rangle\Big\rbrack\\
&=\frac{1}{Z}\sum_{\vec{\sigma}_1,\ldots,\vec{\sigma}_{2M}}
\langle\vec{\sigma}_{1}\vert e^{-\Delta\tau\hami_{1}}\vert\vec{\sigma}_{2M}
\rangle\ldots\langle\vec{\sigma}_{3}\vert e^{-\Delta\tau\hami_{1}}\vert
\vec{\sigma}_{2}\rangle\langle\vec{\sigma}_{2}\vert e^{-\Delta\tau\hami_{2}}
\vert\vec{\sigma}_{1}\rangle\\
&\qquad\underbrace{\Big\lbrack\frac{\langle\vec{\sigma}_{3}\vert e^{-\Delta
\tau\hami_{1}}\hami_1\vert\vec{\sigma}_{2}\rangle}{\langle\vec{\sigma}_{3}
\vert e^{-\Delta\tau\hami_{1}}\vert\vec{\sigma}_{2}\rangle}+\frac{\langle
\vec{\sigma}_{2}\vert e^{-\Delta\tau\hami_{2}}\hami_2\vert\vec{\sigma}_{1}
\rangle}{\langle\vec{\sigma}_{2}\vert e^{-\Delta\tau\hami_{2}}\vert
\vec{\sigma}_{1}\rangle}\Big\rbrack}_{E(w)}\\
&=\frac{\sum_w\Omega(w)E(w)}{\sum_w\Omega(w)}\equiv\lbrace E(w)\rbrace_w
\end{split}
\end{equation}
Die Matrixelemente lassen sich als Summen über (weisse) Plaketten schreiben,
z.B.:
\begin{equation}
\begin{split}
\frac{\langle
\vec{\sigma}_{2}\vert e^{-\Delta\tau\hami_{2}}\hami_2\vert\vec{\sigma}_{1}
\rangle}{\langle\vec{\sigma}_{2}\vert e^{-\Delta\tau\hami_{2}}\vert
\vec{\sigma}_{1}\rangle}&=\sum_{i=1}^{\frac{L}{2}}\frac{
\langle\sigma_{2i,\tau+1},\sigma_{2i+1,\tau+1}\vert e^{-\Delta\tau
\hami^{(2i)}}\hami^{(2i)}\vert\sigma_{2i,\sigma},\sigma{2i+1,\tau}\rangle}{
\langle\sigma_{2i,\tau+1},\sigma_{2i+1,\tau+1}\vert e^{-\Delta\tau
\hami^{(2i)}}\vert\sigma_{2i,\sigma},\sigma{2i+1,\tau}\rangle}\\
&=\sum_{i=1}^{\frac{L}{2}}-\frac{\partial}{\partial\Delta T}
\big\lbrack
\underbrace{\langle\sigma_{2i,\tau+1},\sigma_{2i+1,\tau+1}\vert e^{-\Delta\tau
\hami^{(2i)}}\vert\sigma_{2i,\sigma},\sigma{2i+1,\tau}\rangle}_w\big\rbrack
\end{split}
\end{equation}
\newpage
Zuerst berechnen wir das neue Matrixelement direkt (d.h. gemäß erster Zeile)
für \includegraphics[width=5mm]{bilder/3_13.png}:
\begin{equation}
\begin{split}
 \langle\downarrow\uparrow\vert e^{-\Delta\tau\hami_2}{\color{magenta}\hami_2}
\vert\uparrow\downarrow\rangle&=
\frac{1}{\sqrt{2}}\langle\downarrow\uparrow\vert e^{-\Delta\tau\hami_2}
{\color{magenta}\hami_2}\big\lbrack
\frac{1}{\sqrt{2}}(\ket{\uparrow\downarrow}-\ket{\downarrow\uparrow})+
\frac{1}{\sqrt{2}}(\ket{\uparrow\downarrow}+\ket{\downarrow\uparrow})\big
\rbrack\\
&=\frac{1}{\sqrt{2}}
\langle\downarrow\uparrow\vert\big\lbrack e^{-\Delta\tau(-\frac{J_z}{4}-
\frac{J_x}{2})}
\frac{1}{\sqrt{2}}(\ket{\uparrow\downarrow}-\ket{\downarrow\uparrow})\\
&\qquad+e^{-\Delta\tau(-\frac{J_z}{4}+
\frac{J_x}{2})}
\frac{1}{\sqrt{2}}(\ket{\uparrow\downarrow}+\ket{\downarrow\uparrow})\big
\rbrack\\
&=-e^{\Delta\tau\frac{J_z}{4}}{\color{magenta}\bigg\lbrack}
\sinh\Big(\frac{\Delta\tau J_x}{2}\Big)
{\color{magenta}\Big(-\frac{J_z}{4}\Big)+\cosh\Big(\frac{\Delta\tau J_x}{2}\Big)
\Big( -\frac{J_x}{2}\Big)\bigg\rbrack }
\end{split}
\end{equation}
Allerdings ist es praktischer, die logarithmische Ableitung für die Berechnung
des Quotienten zu benutzen:\\
\includegraphics[width=\linewidth]{bilder/3_14.png}
% \begin{tabular}{l|l}
% \color{darkgreen}$w$&\color{darkgreen}$E_{loc}=-\frac{\partial}{\partial
% \Delta\tau}\ln\lbrack w\rbrack$\\
% $e^{-\Delta\tau\frac{J_z}{4}}$ &$\frac{J_z}{4}$\\
% $e^{\Delta\tau\frac{J_z}{4}}\cosh\big(\Delta\tau\frac{J_x}{2}\big)$&
% $-\frac{J_z}{4}-\frac{J_x}{2}\tanh\big(\Delta\tau\frac{J_x}{2}\big)$\\
% $(-)e^{\Delta\tau\frac{J_z}{4}}\sinh\big(\Delta\tau\frac{J_x}{2}\big)$&
% $-\frac{J_z}{4}-\frac{J_x}{2}\coth\big(\Delta\tau\frac{J_x}{2}\big)$
% \end{tabular}\\
Insgesamt erhalten wir den Erwartungswert der Energie für ein Sample (eine
Konfiguration $w$) als Summe über die Plaketten in den unteren beiden Reihen:
\begin{equation}
 \langle\hami\rangle_{w,einzel}=\sum_{\substack{\textrm{Plakette $p$ mit}\\
\textrm{Zeiten $\tau_1,\tau_2$}\\ \textrm{bzw. }\tau_2,\tau_3}}
E_{loc}(p)
\end{equation}
Um die Varianz zu reduzieren, dürfen wir wegen der Äquivalenz der imaginären
Zeiten über alle Plaketten summieren. Für $N$ Messungen (z.B. nach jeweils
einem vollen Sweep, bei dem alle schwarzen Plaketten betrachtet wurden)
erhalten wir:
\redbox{
\langle\hami\rangle=\frac{1}{N}\sum_{\substack{\textrm{sweeps}\\s=1}}^N
\frac{1}{M}\sum_{\substack{\textrm{alle}\\ \textrm{Platten }p}}
E_{loc}(p,s)
}
\hl{NN-Spinkorrelation:} $\langle\sigma_i^2\sigma_{i+1}^2\rangle=\frac{1}{N}
\sum_s\frac{1}{M}\sum_p C_{loc}(p,s)$ mit\\
$C_{loc}=1$ für \includegraphics[width=20mm]{bilder/3_15.png}\\
$C_{loc}=-1$ sonst

\section{WL-QMC mit Loop-Updates}

Probleme des in \ref{sec:3.3} vorgestellten Weltlinien-QMC-Verfahrens mit
lokalen Updates:
\begin{enumerate}
 \item konstante Anzahl der Weltlinien ($\rightarrow\;S^z$ konstant im
XXZ-Modell, $n$ konstant in Hubbard-Modellen)
\item Autokorellationszeiten divergieren für $\Delta\tau\rightarrow0$
($\Rightarrow\;\Delta\tau$-Extrapolation schwierig+teuer)
\end{enumerate}
Beide Nachteile werden in Schleifen (=Loop)-Algorithmen überwunden.

\subsection{Abbildung auf 6-Vertex-Modell}

Im ersten Schritt werden die Weltlinien-Konfigurationen auf Konfigurationen
des 6-Vertex-Modells abgebildet. Dieses besteht aus Vertizes auf einem um 45
Grad rotierten Gitter (Mittelpunkt der weissen Plaketten), die divergenzfrei
durch Pfeile verbunden sind. Auf der Ebene der Plaketten gilt folgende
Zuordnung:\\
{\centering\includegraphics[width=100mm]{bilder/3_16.png}}\\
Grundsätzlich gilt also, dass aufwärtsgerichtete Pfeile zu Weltlinien
korrespondieren.\\
\\
Anwendungsbeispiel (der Klarheit halber sind beide Systeme übereinander
gezeichnet):
\begin{enumerate}
 \item nur aufwärtsgerichtete Pfeile\\
\includegraphics[width=120mm]{bilder/3_17.png}
\item vollständig vor Update \\
\includegraphics[width=120mm]{bilder/3_18.png}
\item Auswahl einer Schleife \\
\includegraphics[width=120mm]{bilder/3_19.png}\newpage
\item Spin-Flip der Schleife und Weltlinien-Update: \\
\includegraphics[width=120mm]{bilder/3_20.png}\\
Mit dem Update hat sich also (i) die Zahl der Weltlinien erhöht und (ii) ist
jetzt die Windungszahl (in $L$-Richtung) endlich.
\end{enumerate}
Im Prinzip könnte man ein Loop-QMC-Verfahren wie im Beispiel beschrieben
implementieren:
\begin{enumerate}
 \item Initialisierung etc.
\item Für Konfiguration $w$: wähle Loop zufällig $\Rightarrow\;w_{neu}$
\item Berechne $\frac{\Omega(w_{neu})}{\Omega(w)}$
\item Akzeptiere z.B. nach heat bath-Regel
\item Messe Observablen
\item zurück zu 2.
\end{enumerate}
$\;$
\begin{problem}
 \begin{itemize}
 \item Schritt 3. ist i.A. sehr teuer
\item Akzeptanzrate in 4. ist oft niedrig
\end{itemize}
\end{problem}
\begin{swort}{Lösung}
 Verlagere Akzeptanzschritt auf Konstruktion der Schleife (analog zur
Clusterkonstruktion bei Fortuin-Kasteleyn). Führe dazu neue Variablen ein,
sogenannte Graphen, die gegebenen WL- bzw Vertex-Konfigurationen mit einer zu
bestimmenden Wahrscheinlicheit zugewiesen werden. Dabei definiert ein Satz von
zulässigen Graphen eindeutig Schleifen, analog den F-K-Clustern im Ising-Fall.
Umgekehrt sollen zu einem gegebenen Graphen alle kompatiblen
Vertex-Konfigurationen gleichwahrscheinlich sein.
\end{swort}\\ \\
Konkret muss man in diesem Ansatz maximal 4 Plaketten-Graphen einführen:
\\
\includegraphics[width=\linewidth]{bilder/3_21.png}
% "`vertikal"', "`horizontal"', "`eingefroren"', "`diagonal"'\\
Dabei gilt die Kompatibilitätstabelle:\\
\includegraphics[width=\linewidth]{bilder/3_22.png}\\
\\
Zur Logik: die Linien in den Graphen sind ungerichtet und müssen möglichen
Vertex-Pfaden entsprechen. Dies gilt allerdings nicht für den eingefrorenen
Graphen, der dem Umklappen aller Spins entspricht und immer kompatibel ist.
Als Beispiel betrachten wir die 1. Zeile:\\
\includegraphics[width=20mm]{bilder/3_23.png}
Es gibt Pfade $\color{red}A\rightarrow B, A\rightarrow C, D\rightarrow B,
D\rightarrow C$, aber z.B. keine Verbindung $\color{red}A\rightarrow D$ oder
$\color{red}D\rightarrow A$ $\longrightarrow$ kein Graph 
\includegraphics[width=10mm]{bilder/3_24.png}\\
Formalismus: für die betrachtete Plakette sei $S$ eine Spinkonfiguration, $G$
ein Graph und $W(S)$ das in \ref{sec:3.3} bestimmte Gewicht zu $S$.\\
Es sollen nun Gewichte $W(S,G)$ für jede kompatible Konfiguration
$\lbrace S,G\rbrace$ bestimmt werden, so dass gilt:
\begin{enumerate}
 \item $\sum_GW(S,G))W(S)$
\item $W(S,G)=W(S',G)\;\forall S,S'$ kompatibel mit $G$
\end{enumerate}
Wenn dieses lineare Gleichungssystem gelöst ist, lautet der Algorithmus wie
folgt:
\begin{enumerate}
 \item Initialisiere Weltlinien
 \item Weise jeder (weissen) Plakette mit Spinkonfiguration $S$ einen
Graphen $G$ zu mit W':
\begin{equation}
 P(S\rightarrow\lbrace S,G\rbrace)=\frac{W(S,G)}{W(S)}
\end{equation}
\item Identifiziere alle Schleifen (jeweils verbundene Linien). Lege für jede
Schleife die Spinrichtung neu fest (jeweils W' $\frac{1}{2}$).
\item Messe Observablen
\item zurück zu 2.
\end{enumerate}
Dies entspricht dem Swendson-Wang-Algorithmus.\\
Alternative: Wolld-artige Modifikation
\begin{enumerate}
 \item[3'.] Wähle einen Spin und flippe die assoziierte Schleife 
\end{enumerate}
Offensichtlich erfüllen beide Varianten das detaillierte Gleichgewicht.\\
Die Gewichte $W(S,G)$ sind i.A. nicht eindeutig festgelegt. Es empfiehlt sich,
die gefrorenen Graphen zu minimieren, da sonst die Cluster zu groß werden
(und im Fall eines einzign Clusters effektiv nichts passiert).\\
Speziell für das Heisenberg-Modell (XXZ-Modell mit $J_x=J_z$) können
wir auf den gefrorenen Graphen ganz verzichten und auch den Diagonalgraphen
eliminieren:\\
\includegraphics[width=\linewidth]{bilder/3_25.png}
% \begin{equation}
% \begin{split}
% e^{\Delta\tau\frac{J_z}{4}}\cosh\big(\Delta\tau\frac{J_x}{2}\big)
% &\equiv W_1=W_{1,1}+W_{1,2}\\
% e^{\Delta\tau\frac{J_z}{4}}\sinh\big(\Delta\tau\frac{J_x}{2}\big)
% &\equiv W_2=W_{2,2}\\
% e^{-\Delta\tau\frac{J_z}{4}} &\equiv W_3=W_{3,1}
% \end{split}
% \end{equation}
Dabei folgt aus der Gleichverteilungsbedingung:
\begin{equation}
\begin{split}
W_{1,1}&\stackrel{!}{=}W_{3,1}=W_3\\
W_{1,2}&\stackrel{!}{=}W_{2,2}=W_2
\end{split}
\end{equation}
Dies ist konsistent, da
\begin{equation}
\begin{split}
W_{1,1}+W_{1,2}&=W_3+W_2\\
&=e^{\Delta\tau\frac{J}{4}}\big(e^{-\Delta\tau\frac{J}{2}}+
\sinh(\Delta\tau\frac{J}{2})\big)\\
&=e^{\delta\tau\frac{J}{4}}\frac{2e^{-\Delta\tau\frac{J}{2}}
+e^{\Delta\tau\frac{J}{2}}-e^{-\Delta\tau\frac{J}{2}}}{2}\\
&=e^{\delta\tau\frac{J}{4}}\cosh\big(\frac{\Delta\tau J}{2}\big)=W_1
\end{split}
\end{equation}
Abschließend noch ein Beispiel:
\begin{enumerate}
 \item erzwungene Graphen\\
\includegraphics[width=140mm]{bilder/3_26.png}
\item alle Graphen mit \includegraphics[width=30mm]{bilder/3_29.png}\\
\includegraphics[width=140mm]{bilder/3_27.png}\newpage
\item Update hier: größter Cluser umgeklappt\\
\includegraphics[width=140mm]{bilder/3_28.png}
\end{enumerate}