Notes.tex

\documentclass[12pt]{article}

% packages
\usepackage{setspace}
\usepackage{tikz}
\usepackage{array}
\usepackage[margin=0.75in]{geometry}
\usepackage{amsmath,bm}
\usepackage{amssymb}
\usepackage{bbold}
\usepackage{physics}
\usepackage{xcolor}
\usepackage{indentfirst}
\usepackage{enumerate}
\usepackage{mathtools}
\usepackage{cancel}

\doublespacing
\allowdisplaybreaks

\title{Notes on SSE for TFIM + $Z$-Field (LTFIM) at Finite Temperature}
\author{Isaac De Vlugt}
\date{\today}
\begin{document}

\maketitle

\section{Introduction}

This is motivated by our need to make an SSE QMC for the Rydberg Hamiltonian.
\begin{align}
    H = & \sum_{i=1}^N \frac{\Omega}{2} \left(\ketbra{g}{r}_i + \ketbra{r}{g}_i\right) - \delta \sum_{i=1}^N \ketbra{r}{r}_i + \frac{1}{2} \sum_{\langle i,j \rangle} V_{ij} \ketbra{r}{r}_i \bigotimes \ketbra{r}{r}_j
\end{align}
Upon identifying $\ket{g},\ket{r}$ with eigenstates of $\sigma^z$ ($\ket{0},\ket{1}$), we can translate this Hamiltonian into spin language.
%\begin{align}
    %H = & \sum_{i=1}^N \frac{\Omega}{2} \left(\ketbra{0}{1}_i + \ketbra{1}{0}_i\right) - \delta \sum_{i=1}^N \ketbra{1}{1}_i + \frac{1}{2} \sum_{\langle i,j \rangle} V_{ij} \ketbra{1}{1}_i \bigotimes \ketbra{1}{1}_j \\
    %=& \sum_{i=1}^N \frac{\Omega}{2} \sigma^x_i + \delta \sum_{i=1}^N \left(\sigma^z_i - \ketbra{0}{0}_i\right) + \frac{1}{2} \sum_{\langle i,j \rangle} V_{ij} \left(\sigma^z_i - \ketbra{0}{0}_i\right) \bigotimes \left(\sigma^z_j - \ketbra{0}{0}_j\right) \\
    %=& 
%\end{align}
%Here, $n_i$ is simply an occupation number (i.e. either 0 or 1). We can therefore trivially map this to spin-$\frac{1}{2}$ with $n_i = \frac{1}{2}\left(\sigma_i^z + \mathbb{1} \right)$. The hamiltonian then becomes the following.
\begin{align}
    H = \frac{\Omega}{2} \sum_{i=1}^N \left( \ketbra{0}{1}_i + \ketbra{1}{0}_i \right) - \delta \sum_{i=1}^N \ketbra{1}{1}_i + \frac{1}{2} \sum_{(i,j)} V_{ij} \ketbra{1_i 1_j}{1_i 1_j}
\end{align}
This is very similar to a transverse-field Ising model with an additional longitudinal field.

\section{Longitudinal Field TFIM (LTFIM)}

The Rydberg Hamiltonian is extremely similar to a LTFIM (general $J_{ij}$ for all pairs, assuming $J_{ij} > 0 \quad \forall \quad i,j$ as in the Rydberg case since $V_{ij} = \frac{C}{R_{ij}^6}$).
\begin{align} \label{eq:LTFIM_hamiltonian}
    H = & - h_x \sum_{i = 1}^N \sigma^x_i - h_z \sum_{i=1}^N \sigma_i^z + \sum_{\expval{i,j}} J_{ij} \sigma_i^z \sigma_j^z \qquad \mathrm{LTFIM}
\end{align}
Thus, if an efficient SSE formalism exists for such a Hamiltonian, then it is foreseeable that one may exist for the Rydberg Hamiltonian. So, let's work with the LTFIM as a toy model for the Rydberg Hamiltonian.

\subsection{Decomposition of Hamiltonian} \label{sec:H_decomp}

We first must decompose the Hamiltonian in Eq.~\eqref{eq:LTFIM_hamiltonian} in the form
\begin{align}
    H = - \sum_{t,a} H_{t,a},
\end{align}
where $t$ denotes the the type of $H_{t,a}$ (usually off-diagonal versus diagonal operators) and $a$ denotes the lattice unit over which the operator acts (i.e. sites, bonds, plaquettes, etc).
Inspiration can be taken from Sandvik's directed loops paper, wherein the XXZ model is explored in great detail. 
\begin{subequations} \label{eq:decomposition}
    \begin{align}
        H_{0,0} =& \mathbb{1} \\
        H_{-1,a} =& h_x \sigma_i^x \\
        H_{1,a} =& h_x \\
        H_{1,b} =& J \sigma_i^z \sigma_j^z + h_{z,b} (\sigma_i^z + \sigma_j^z) + C
    \end{align}
\end{subequations}
Here, 
\begin{equation}
    h_{z,b} = \frac{h_z}{2D}
\end{equation}
is the modified logitudinal field strength to account for over counting when the sum over sites is put into the sum over bonds ($D$ is the lattice dimension).
This general expression accounts for PBCs and OBCs.
$C = \max(J, 2 h_{z,b} - J) + \varepsilon$ is a constant to aleviate sign problems with $H_{1,b}$, and $\varepsilon$ is a small non-negative constant to aid with numerics. 
In all, we have added a constant shift to Eq.~\eqref{eq:LTFIM_hamiltonian} of $N_b C + Nh$ that must be accounted for when calculating the energy. 

The non-zero matrix elements of all parts of the Hamiltonian (Eq.~\eqref{eq:decomposition}) are
\begin{subequations} \label{eq:mat_elems}
    \begin{align}
        \matrixel{\uparrow}{H_{-1,a}}{\downarrow} = \matrixel{\downarrow}{H_{-1,a}}{\uparrow} =& h_x \\
        \matrixel{\uparrow}{H_{1,a}}{\uparrow} = \matrixel{\downarrow}{H_{1,a}}{\downarrow} =& h_x \\
        W_1 \equiv \matrixel{\uparrow \uparrow}{H_{1,b}}{\uparrow \uparrow} =& J + 2h_{z,b} + C \\
        W_2 \equiv \matrixel{\downarrow \downarrow}{H_{1,b}}{\downarrow \downarrow} =& J - 2h_{z,b} + C \\
        W_3 \equiv \matrixel{\downarrow \uparrow}{H_{1,b}}{\downarrow \uparrow} = \matrixel{\uparrow \downarrow}{H_{1,b}}{\uparrow \downarrow} =& -J + C
    \end{align}
\end{subequations}

\subsection{Diagonal Updates}

We may have to toy around with the following procedures, but the general flow of things should be the same no matter what.
For proposal modifications to the simulation cell, I will use primes.

\subsubsection{Finite Temperature}

Traverse the operator list and do the following.

\begin{enumerate}
    \item If $H_{1,a}$ or $H_{1,b}$ is encountered remove it ($n \rightarrow n - 1$) with probability
    \begin{align*}
        P_{[1,a]_p / [1,b]_p \rightarrow [0,0]_p} =& \mathrm{min}\left( \frac{W^\prime(\alpha(p), S_M)}{ W(H_{[1,a]_p} \text{ or } H_{[1,b]_{p,\uparrow\uparrow}} \text{ or } H_{[1,b]_{p,\downarrow\downarrow}} \text{ or } H_{[1,b]_{p,\uparrow\downarrow / \downarrow\uparrow}}) }, 1 \right) \\
        =& \mathrm{min}\left( \frac{W^\prime(\alpha(p), S_M)}{NW(H_{[1,a]_p}) + N_b(W(H_{[1,b]_{p,\uparrow\uparrow}}) + W(H_{[1,b]_{p,\downarrow\downarrow}}) + W(H_{[1,b]_{p,\uparrow\downarrow / \downarrow\uparrow}})) }, 1 \right) \\
        =& \mathrm{min}\left(\frac{M-n+1}{\beta \left(Nh + N_b\left(W_1 + W_2 + 2W_3\right)\right) },1\right).
    \end{align*}
    The factor of $N$ or $N_b$ in the denominator is from the fact that where the operator $H_{1,a}$ or $H_{1,b}$ is when we remove it does not matter. 
    \item If the unity operator $H_{0,0}$ is encountered, decide whether or not to accept inserting a diagonal operator ($n \rightarrow n + 1$) with the probability 
    \begin{align}
        P_{H_{[0,0]_p} \rightarrow H_{[1,a]_p} / H_{[1,b]_p}} =& \mathrm{min}\left(\frac{\beta \left(Nh + N_b\left(W_1 + W_2 + 2W_3\right)\right) }{M-n},1\right)
    \end{align}
    \item If it was decided to insert a diagonal operator, choose the operator to insert with probabilities
    \begin{subequations} \label{eq:heat_bath_1}
        \begin{align} 
            P_{H_{[1,a]_p}} = & \frac{Nh}{Nh + N_b\left(W_1 + W_2 + 2W_3\right)}, \\
            P_{H_{[1,b]_p}} = & \frac{N_b\left(W_1 + W_2 + 2W_3\right)}{Nh + N_b\left(W_1 + W_2 + 2W_3\right)}.
        \end{align}
    \end{subequations}
    \item If it was decided to insert $H_{1,a}$, then choose a random site (probability $1 / N$) to put the operator.
    \item If it was decided to insert $H_{1,b}$, first sample a bond, $b^\prime$, based on the matrix elements $W_1$, $W_2$, and $W_3$. \label{step:bond_insertion}
    Specifically, 
    \begin{align*}
        P_{\uparrow \uparrow} =& \frac{W_1}{W_1 + W_2 + 2W_3} \\
        P_{\downarrow \downarrow} =& \frac{W_2}{W_1 + W_2 + 2W_3} \\
        P_{\uparrow \downarrow} =& \frac{W_3}{W_1 + W_2 + 2W_3} \\
        P_{\downarrow \uparrow} =& \frac{W_3}{W_1 + W_2 + 2W_3}
    \end{align*}
    Then, randomly select a bond $b^{\prime\prime}$ in the actual propagated spin state.
    Accept the insertion of $H_{1,b}$ at the randomly chosen bond $b^{\prime\prime}$ with probability
    \begin{align*}
        \frac{W_{b^{\prime\prime}}}{W_{b^{\prime}}}.
    \end{align*}
                
    {\color{red} {\bf This doesn't work:}
    If it was decided to insert $H_{1,b}$, choose a random bond to insert it.
    In the case when one of $W_1$, $W_2$, or $W_3$ are zero and that corresponding bond is chosen as the insertion location, the insertion must be rejected in its entirety.}
\end{enumerate}

\subsubsection{Zero Temperature}

Since we want to maximize the expansion order here to project out the ground state of the Hamiltonian, we want to avoid using identity operators and force our simulation cell to be filled with local terms of our Hamiltonian.
So, traverse the operator list and do the following.
\begin{enumerate}
    \item Remove the current operator $H_{1,a}$ or $H_{1,b}$.
    \item A new operator type is selected corresponding to probabilities in Eq.~\eqref{eq:heat_bath_1}.
    \item If $H_{1,a}$ is chosen, select a random site and place it there.
    \item If $H_{1,b}$ is chosen, follow step~\ref{step:bond_insertion} as in the finite temperature case.
\end{enumerate}

\subsection{Loop Updates}

In the spirit of the TFIM SSE formalism, we may still terminate on site operators $H_{1,a}$ and $H_{-1,a}$. 
However, loop formation does not have to be deterministic since the matrix elements of $H_{1,b}$ are not all the same and a Swendsen-Wang (SW) update procedure would be invalid.
Therefore, we must either devise a non-deterministic loop formation algorithm wherein detailed balance (DB) is satisfied, or form deterministic loops as in the case of the TFIM SSE formalism but flip the cluster by accounting for the weight change (unlike a traditional SW procedure).

\subsubsection{Detailed balance for non-deterministic loops}

Recall the DB requirement given a configuration $s$ and $s^\prime$.
\begin{align} \label{eq:detailed_balance}
    P(s \rightarrow s^\prime)W(s) = P(s^\prime \rightarrow s)W(s^\prime)
\end{align}
It turns out that for detailed balance to be satisfied when loop segments reach a vertex associated with $H_{1,b}$, the following two equations must be satisfied.
\begin{subequations} \label{eq:db_requirements}
    \begin{align}
        W(s,\ell_1,\ell_2) =& W(s^\prime, \ell_2, \ell_1) \label{eq:db_1} \\
        \sum_x W(s,e,x) =& W_s \label{eq:db_2}
    \end{align}
\end{subequations}
\begin{itemize}
    \item $s$, $s^\prime$: vertex with it's four legs specified (legs being the spin states surrounding the operator $H_{1,b}$) 
    \item $e$: the leg that the loop segment enters
    \item $x$: the leg that the loop segment exits
    \item $\sum_x$: a sum over all exit paths (given a vertex $s$ and an entrace leg $e$)
    \item $W_s$: the matrix element of a single vertex with its four legs encoded as $s$
    \item $W(s,e,x) = W_s P(s,e \rightarrow s^\prime, x)$
\end{itemize}
Eq.~\eqref{eq:db_requirements} define the direted loop equations that must be satisfied for DB. 
It should also be noted that although multi-branching loop segments should be allowed for non-deterministic loop formation, they are excluded in this discussion for simplicity. 

\subsubsection{Directed loop equations}

To obtain directed loop equations and devise ways to satisfy them, we need to first look at possible loop segments through vertices corresponding to matrix elements of $H_{1,b}$.
Following Sandvik's directed loops paper, any ``switch'' loop segment (i.e. a segment whose exit leg is not the same lattice site as the entrance leg) will not be allowed since flipping the entrance and exit legs will yield off-diagonal matrix element of $H_{1,b}$, which is non-existant ($H_{1,b}$ is a diagonal operator).
Therefore, the only two loop segments that are allowed are ``continue straight'' or ``bounce''. See Fig.~\ref{fig:continue_bounce}. 

\begin{figure}
    \centering
    \begin{tikzpicture}

        % new

        \filldraw[black] (0.5, 1) node[anchor=center] {continue straight};

        \draw[thick] (0.1,0.4) circle (0.1); % UL
        \draw[thick] (0.9,0.4) circle (0.1); % UR
        \draw[thick] (0.1,-0.2) circle (0.1); % LL
        \draw[thick] (0.9,-0.2) circle (0.1); % LR

        \draw[black, fill=black] (0,0) rectangle (1,0.2);

        \draw[red, thick, ->] (0.1,-0.5) -- (0.1,0.7);

        % new

        \filldraw[black] (4.5, 1) node[anchor=center] {bounce};

        \draw[thick] (4.1,0.4) circle (0.1); % UL
        \draw[thick] (4.9,0.4) circle (0.1); % UR
        \draw[thick] (4.1,-0.2) circle (0.1); % LL
        \draw[thick] (4.9,-0.2) circle (0.1); % LR

        \draw[black, fill=black] (4,0) rectangle (5,0.2);

        \draw[red, thick] (4,-0.5) -- (4,-0.1);
        \draw[red, thick] (4,-0.1) arc (180:0:0.1);
        \draw[red, thick, ->] (4.2,-0.1) -- (4.2,-0.5);

    \end{tikzpicture}    
\caption{...} \label{fig:continue_bounce}
\end{figure}

For loop segments to satisfy DB, we must relate vertices in which the two spins that are connected by the segment are flipped and the direction of the segment if reversed. 
An example is in Fig.~\ref{fig:related_vertices}.

\begin{figure} [!h]
    \centering
    \begin{tikzpicture}

        % new

        \draw[thick] (0.1,0.4) circle (0.1); % UL
        \draw[thick] (0.9,0.4) circle (0.1); % UR
        \draw[thick] (0.1,-0.2) circle (0.1); % LL
        \draw[thick] (0.9,-0.2) circle (0.1); % LR

        \draw[black, fill=black] (0,0) rectangle (1,0.2);

        \draw[red, thick, ->] (0.1,-0.5) -- (0.1,0.7);

        % new

        \draw[thick, fill=black] (4.1,0.4) circle (0.1); % UL
        \draw[thick] (4.9,0.4) circle (0.1); % UR
        \draw[thick, fill=black] (4.1,-0.2) circle (0.1); % LL
        \draw[thick] (4.9,-0.2) circle (0.1); % LR

        \draw[black, fill=black] (4,0) rectangle (5,0.2);

        \draw[red, thick, ->] (4.1,0.7) -- (4.1,-0.5);

    \end{tikzpicture}    
\caption{...} \label{fig:related_vertices}
\end{figure}
Mathematically, these graphs represent the equivalence relation Eq.~\ref{eq:db_1}.
Note that the bounce loop segments map to themselves under this transformation (flipping spins, reversing loop segment direction).
Furthermore, Eq.~\ref{eq:db_2} relates vertices with different exit legs having the same spin configuration and entrance leg.
Graphically, the sum of the graphs in Fig.~\ref{fig:continue_bounce} must be equal to $W_3$ (here, we use filled/unfilled dots to represent spin-up/down).

Now, let's look at the directed loop equations. 
It's at this point that where we refer to the directed loops paper, specifically Fig. 8.
This figure outlines that all possible vertex configurations can be divided into 8 subsets that don't transform into each other as dictated by Eq.~\ref{eq:db_1}.
It's a little different/simpler in our case (i.e. XXZ model versus LTFIM), since we cannot have any "switch" loop segments.
The other 4 subsets of related vertices are simply those in Fig.~\ref{fig:vertex_subsets} but with the spins flipped.
Each configuration in a given quadrant maps to another configuration in the same quadrant (i.e. similar to Fig.~\ref{fig:related_vertices}). 
Again, bounce operations map to themselves.

\begin{figure}
    \centering
    \begin{tikzpicture}

        \draw[gray, thick] (0,-3) -- (8,-3);
        \draw[gray, thick] (4, 1) -- (4, -7);

        % new 1,1

        \draw[thick] (0.1,0.4) circle (0.1); % UL
        \draw[fill=black,thick] (0.9,0.4) circle (0.1); % UR
        \draw[thick] (0.1,-0.2) circle (0.1); % LL
        \draw[fill=black,thick] (0.9,-0.2) circle (0.1); % LR

        \draw[black, fill=black] (0,0) rectangle (1,0.2);

        \draw[red, thick] (0.8,-0.5) -- (0.8,-0.1);
        \draw[red, thick] (0.8,-0.1) arc (180:0:0.1);
        \draw[red, thick, ->] (1,-0.1) -- (1,-0.5);

        % new 1,2

        \draw[thick] (2.1,0.4) circle (0.1); % UL
        \draw[fill=black,thick] (2.9,0.4) circle (0.1); % UR
        \draw[thick] (2.1,-0.2) circle (0.1); % LL
        \draw[fill=black,thick] (2.9,-0.2) circle (0.1); % LR

        \draw[black, fill=black] (2,0) rectangle (3,0.2);

        \draw[red, thick, ->] (2.9,-0.5) -- (2.9,0.7);

        % new 1,3

        \draw[thick] (0.1,-1.6) circle (0.1); % UL
        \draw[thick] (0.9,-1.6) circle (0.1); % UR
        \draw[thick] (0.1,-2.2) circle (0.1); % LL
        \draw[thick] (0.9,-2.2) circle (0.1); % LR

        \draw[black, fill=black] (0,-2) rectangle (1,-1.8);

        \draw[red, thick, ->] (0.9,-1.3) -- (0.9,-2.5);

        % new 1,4

        \draw[thick] (2.1,-1.6) circle (0.1); % UL
        \draw[thick] (2.9,-1.6) circle (0.1); % UR
        \draw[thick] (2.1,-2.2) circle (0.1); % LL
        \draw[thick] (2.9,-2.2) circle (0.1); % LR

        \draw[black, fill=black] (2,-2) rectangle (3,-1.8);

        \draw[red, thick] (2.8,-1.7) -- (2.8,-1.3);
        \draw[red, thick] (2.8,-1.7) arc (180:360:0.1);
        \draw[red, thick, ->] (3,-1.7) -- (3,-1.3);

        % new 2,1

        \draw[thick] (5.1,0.4) circle (0.1); % UL
        \draw[fill=black,thick] (5.9,0.4) circle (0.1); % UR
        \draw[thick] (5.1,-0.2) circle (0.1); % LL
        \draw[fill=black,thick] (5.9,-0.2) circle (0.1); % LR

        \draw[black, fill=black] (5,0) rectangle (6,0.2);

        \draw[red, thick] (5,-0.5) -- (5,-0.1);
        \draw[red, thick] (5,-0.1) arc (180:0:0.1);
        \draw[red, thick, ->] (5.2,-0.1) -- (5.2,-0.5);

        % new 2,2

        \draw[thick] (7.1,0.4) circle (0.1); % UL
        \draw[fill=black,thick] (7.9,0.4) circle (0.1); % UR
        \draw[thick] (7.1,-0.2) circle (0.1); % LL
        \draw[fill=black,thick] (7.9,-0.2) circle (0.1); % LR

        \draw[black, fill=black] (7,0) rectangle (8,0.2);

        \draw[red, thick, ->] (7.1,-0.5) -- (7.1,0.7);

        % new 2,3

        \draw[fill=black,thick] (5.1,-1.6) circle (0.1); % UL
        \draw[fill=black,thick] (5.9,-1.6) circle (0.1); % UR
        \draw[fill=black,thick] (5.1,-2.2) circle (0.1); % LL
        \draw[fill=black,thick] (5.9,-2.2) circle (0.1); % LR

        \draw[black, fill=black] (5,-2) rectangle (6,-1.8);

        \draw[red, thick, ->] (5.1,-1.3) -- (5.1,-2.5);

        % new 2,4

        \draw[fill=black,thick] (7.1,-1.6) circle (0.1); % UL
        \draw[fill=black,thick] (7.9,-1.6) circle (0.1); % UR
        \draw[fill=black,thick] (7.1,-2.2) circle (0.1); % LL
        \draw[fill=black,thick] (7.9,-2.2) circle (0.1); % LR

        \draw[black, fill=black] (7,-2) rectangle (8,-1.8);

        \draw[red, thick] (7,-1.7) -- (7,-1.3);
        \draw[red, thick] (7,-1.7) arc (180:360:0.1);
        \draw[red, thick, ->] (7.2,-1.7) -- (7.2,-1.3);

        % new 3,1

        \draw[fill=black,thick] (0.1,-3.6) circle (0.1); % UL
        \draw[thick] (0.9,-3.6) circle (0.1); % UR
        \draw[fill=black,thick] (0.1,-4.2) circle (0.1); % LL
        \draw[thick] (0.9,-4.2) circle (0.1); % LR

        \draw[black, fill=black] (0,-4) rectangle (1,-3.8);

        \draw[red, thick] (0.8,-3.7) -- (0.8,-3.3);
        \draw[red, thick] (0.8,-3.7) arc (180:360:0.1);
        \draw[red, thick, ->] (1,-3.7) -- (1,-3.3);

        % new 3,2

        \draw[fill=black,thick] (2.1,-3.6) circle (0.1); % UL
        \draw[thick] (2.9,-3.6) circle (0.1); % UR
        \draw[fill=black,thick] (2.1,-4.2) circle (0.1); % LL
        \draw[thick] (2.9,-4.2) circle (0.1); % LR

        \draw[black, fill=black] (2,-4) rectangle (3,-3.8);

        \draw[red, thick, ->] (2.9,-3.3) -- (2.9,-4.5);

        % new 3,3

        \draw[fill=black,thick] (0.1,-5.6) circle (0.1); % UL
        \draw[fill=black,thick] (0.9,-5.6) circle (0.1); % UR
        \draw[fill=black,thick] (0.1,-6.2) circle (0.1); % LL
        \draw[fill=black,thick] (0.9,-6.2) circle (0.1); % LR

        \draw[black, fill=black] (0,-6) rectangle (1,-5.8);

        \draw[red, thick, ->] (0.9,-6.5) -- (0.9,-5.3);

        % new 3,4

        \draw[fill=black,thick] (2.1,-5.6) circle (0.1); % UL
        \draw[fill=black,thick] (2.9,-5.6) circle (0.1); % UR
        \draw[fill=black,thick] (2.1,-6.2) circle (0.1); % LL
        \draw[fill=black,thick] (2.9,-6.2) circle (0.1); % LR

        \draw[black, fill=black] (2,-6) rectangle (3,-5.8);

        \draw[red, thick] (2.8,-6.5) -- (2.8,-6.1);
        \draw[red, thick] (2.8,-6.1) arc (180:0:0.1);
        \draw[red, thick, ->] (3,-6.1) -- (3,-6.5);

        % new 4,1

        \draw[fill=black,thick] (5.1,-3.6) circle (0.1); % UL
        \draw[thick] (5.9,-3.6) circle (0.1); % UR
        \draw[fill=black,thick] (5.1,-4.2) circle (0.1); % LL
        \draw[thick] (5.9,-4.2) circle (0.1); % LR

        \draw[black, fill=black] (5,-4) rectangle (6,-3.8);

        \draw[red, thick] (5,-4.5) -- (5,-4.1);
        \draw[red, thick] (5,-4.1) arc (180:0:0.1);
        \draw[red, thick, ->] (5.2,-4.1) -- (5.2,-4.5);

        % new 4,2

        \draw[fill=black,thick] (7.1,-3.6) circle (0.1); % UL
        \draw[thick] (7.9,-3.6) circle (0.1); % UR
        \draw[fill=black,thick] (7.1,-4.2) circle (0.1); % LL
        \draw[thick] (7.9,-4.2) circle (0.1); % LR

        \draw[black, fill=black] (7,-4) rectangle (8,-3.8);

        \draw[red, thick, ->] (7.1,-4.5) -- (7.1,-3.3);

        % new 4,3

        \draw[thick] (5.1,-5.6) circle (0.1); % UL
        \draw[thick] (5.9,-5.6) circle (0.1); % UR
        \draw[thick] (5.1,-6.2) circle (0.1); % LL
        \draw[thick] (5.9,-6.2) circle (0.1); % LR

        \draw[black, fill=black] (5,-6) rectangle (6,-5.8);

        \draw[red, thick, ->] (5.1,-5.3) -- (5.1,-6.5);

        % new 4,4

        \draw[thick] (7.1,-5.6) circle (0.1); % UL
        \draw[thick] (7.9,-5.6) circle (0.1); % UR
        \draw[thick] (7.1,-6.2) circle (0.1); % LL
        \draw[thick] (7.9,-6.2) circle (0.1); % LR

        \draw[black, fill=black] (7,-6) rectangle (8,-5.8);

        \draw[red, thick] (7,-5.7) -- (7,-5.3);
        \draw[red, thick] (7,-5.7) arc (180:360:0.1);
        \draw[red, thick, ->] (7.2,-5.7) -- (7.2,-5.3);

    \end{tikzpicture}    
\caption{...} \label{fig:vertex_subsets}
\end{figure}

Now this is where I am a little shady on some of the details. 
In the directed loops paper, Sandvik notes that there are two symmetries we can take advantage of in order to reduce the number of directed loop equations that we need to solve.
\begin{enumerate}
    \item Permutation of two spins acted on by $H_{1,b}$ (i.e. interchanging spins while keeping the orientation of the loop segment)
    \item Imaginary-time symmetry (swap spins above and below the black-filled bar that represents $H_{1,b}$)
\end{enumerate}
In all, this simplifies our directed loop equations to just those associated to the upper-left and lower-left quadrants (i.e. all other quadrants can be mapped back to the upper-left quadrant or lower-left quadrant). 

The directed loop equations for the upper-left quadrant are
\begin{subequations} \label{eq:directed_loop_eq1}
    \begin{align}
        W_3 = -J + C =& W(3,4,4) + W(3,4,2) = b_1 + a, \\
        W_2 = J - 2h_{z,b} + C =& W(2,2,4) + W(2,2,2) = a + b_2.
    \end{align}
\end{subequations}
Note that we used Eq.~\ref{eq:db_1} to equate $W(3,4,2) = W(2,2,4) = a$.

For the lower-left quadrant, 
\begin{subequations} \label{eq:directed_loop_eq2}
    \begin{align}
        W_3 = -J + C =& W(3,2,2) + W(3,2,4) = b_3 + d, \\
        W_1 = J + 2h_{z,b} + C =& W(1,4,2) + W(1,4,4) = d + b_4.
    \end{align}
\end{subequations}
Note that we used Eq.~\ref{eq:db_1} to equate $W(3,2,4) = W(1,4,2) = d$.
The variables $b_i$ denote bounce weights.

\subsubsection{Heat-bath solutions to directed loop equations}

If we reach the operator $H_{1,b}$ at a given leg, the probability of choosing an exit is given by the weight of the vertex with the spins connected by the loop segment being flipped. 
This is a standard heat-bath update scheme that satisfies detailed balance (proof in directed loops paper).
So, let's verify that these are indeed solutions to our directed loop equations.
I.e., it should be that
\begin{align*}
    \frac{b_1}{W_3} =& \frac{W_3^2}{W_3 + W_2} \\
    \frac{b_2}{W_2} =& \frac{W_2^2}{W_3 + W_2} \\
    \frac{a}{W_3} =& \frac{W_3 W_2}{W_3 + W_2},
\end{align*}
and
\begin{align*}
    \frac{b_3}{W_3} =& \frac{W_3^2}{W_1 + W_3} \\
    \frac{b_4}{W_1} =& \frac{W_1^2}{W_1 + W_3} \\
    \frac{d}{W_3} =& \frac{W_3 W_1}{W_1 + W_3},
\end{align*}
You can check that indeed, these are solutions.

\subsubsection{General solutions to directed loop equations}

While the heat-bath probabilties should work, it is not preferred that the bounce probabilities are non-zero since bouncing in effect does nothing to update the simulation cell.
Let's analyze how $b_i$ can be minimized whilst Eqs.~\ref{eq:directed_loop_eq1} and \ref{eq:directed_loop_eq2} are satisfied and $b_i, a, d \geq 0$ (no sign problems).
Let's solve for the $b_i$'s,

\begin{subequations}
    \begin{align}
        b_1 =& -a - J + C, \\
        b_2 =& -a + J - 2h_{z,b} + C,
    \end{align}
\end{subequations}
and
\begin{subequations}
    \begin{align}
        b_3 =& -d - J + C, \\
        b_4 =& -d + J + 2h_{z,b} + C.
    \end{align}
\end{subequations}

Clearly, the lines defining $b_{1,2}$ vs $a$ are parallel and will never intersect, similarily for the lines defining $b_{3,4}$ vs $d$.
So there doesn't exist a solution to Eqs.~\ref{eq:directed_loop_eq1} and \ref{eq:directed_loop_eq2} wherein $b_i = 0$ $\forall$ $b_i$.
Clearly, this is not favorable. 
Bounce probabilities will always persist and formed loops, flipped or not, won't update the simulation cell and advance the simulation. 

\subsubsection{Deterministic loops}

Instead of the forseeably inefficient non-deterministic directed loops, we may instead form deterministic clusters in the same exact way as is outlined in the well-known SSE formalism for the transverse-field Ising model.
\begin{itemize}
    \item Clusters terminate on site operators.
    \item When a bond operator is encountered, the cluster branches from every possible leg of the bond vertex.
\end{itemize}
However, the difference here is that instead of flipping with probability 1/2, we now flip each cluster, labelled by an index $c$, in the simulation cell with the following Metropolis probability.
\begin{align}
    P_{\text{flip},c} =& \min\left(\frac{W_c^\prime}{W_c}, 1\right)
\end{align}
Here, $W_c$ primed (unprimed) denotes the weight of the local cluster $c$ being flipped (unflipped).
Specifically, $W_c$ (primed or unprimed) is the product of all matrix element weights ($W_1$, $W_2$, $W_3$, or $h$) present in the local cluster $c$. 

\section{Rydberg Hamiltonian}

\begin{align*}
    H =& \frac{\Omega}{2} \sum_{i=1}^N \left( \ketbra{0}{1}_i + \ketbra{1}{0}_i \right) - \delta \sum_{i=1}^N \ketbra{1}{1}_i + \sum_{(i,j)} V_{ij} \ketbra{1_i 1_j}{1_i 1_j} \\
    =& \Omega \sum_{i=1}^N \sigma_i^x + \sum_{(i,j)} V_{ij} \ketbra{1_i 1_j}{1_i 1_j} - \delta_b \left(\ketbra{1}{1}_i + \ketbra{1}{1}_j\right)
\end{align*}
Here, $\delta_b$ is a modified version of $\delta$ that accounts for over-counting when moving the original sum over sites into the sum over all-to-all bonds.
Let's decompose this Hamiltonian into local operators in a similar spirit to that of the LTFIM. 
Also, for simplicity, we redefined $\Omega$ as $2\Omega$ and $V$ as $2V$ just to get rid of pesky factors of 1/2.
This gives us the following local operator decomposition.
\begin{subequations} \label{eq:decomposition_ryd}
    \begin{align}
        H_{0,0} =& \mathbb{1} \\
        H_{-1,a} =& \Omega \sigma_i^x \\
        H_{1,a} =& \Omega \\
        H_{1,b} =& V_{ij} \ketbra{1_i 1_j}{1_i 1_j} - \delta_b \left(\ketbra{1}{1}_i + \ketbra{1}{1}_j\right) + C_{ij}
    \end{align}
\end{subequations}
$C_{ij} = \max(2\delta_b - V_{ij}, \delta_b)$ is a constant shift that is needed in order to alleviate sign problems.
The relevant matrix elements of this local decomposition are (taking $1 (0) = \uparrow (\downarrow)$)
\begin{subequations} \label{eq:mat_elems_ryd}
    \begin{align}
        \matrixel{\uparrow}{H_{-1,a}}{\downarrow} = \matrixel{\downarrow}{H_{-1,a}}{\uparrow} =& \Omega \\
        \matrixel{\uparrow}{H_{1,a}}{\uparrow} = \matrixel{\downarrow}{H_{1,a}}{\downarrow} =& \Omega \\
        W_{ij}^{(1)} \equiv \matrixel{\uparrow \uparrow}{H_{1,b}}{\uparrow \uparrow} =& V_{ij} - 2\delta_b + C_{ij} \\
        W_{ij}^{(2)} \equiv \matrixel{\downarrow \downarrow}{H_{1,b}}{\downarrow \downarrow} =& C_{ij}\\
        W_{ij}^{(3)} \equiv \matrixel{\downarrow \uparrow}{H_{1,b}}{\downarrow \uparrow} = \matrixel{\uparrow \downarrow}{H_{1,b}}{\uparrow \downarrow} =& -\delta_b + C_{ij}
    \end{align}
\end{subequations}

\subsection{Diagonal Updates}

We may have to toy around with the following procedures, but the general flow of things should be the same no matter what.
For proposal modifications to the simulation cell, I will use primes as before.

\subsubsection{Finite Temperature}

Traverse the operator list and do the following.

\begin{enumerate}
    \item If $H_{1,a}$ or $H_{1,b}$ is encountered remove it ($n \rightarrow n - 1$) with probability
    \begin{align*}
        P_{[1,a]_p / [1,b]_p \rightarrow [0,0]_p} =& \mathrm{min}\left( \frac{W^\prime(\alpha(p), S_M)}{ W(H_{[1,a]_p} \text{ or } H_{[1,b]_{p,\uparrow\uparrow}} \text{ or } H_{[1,b]_{p,\downarrow\downarrow}} \text{ or } H_{[1,b]_{p,\uparrow\downarrow / \downarrow\uparrow}}) }, 1 \right) \\
        =& \mathrm{min}\left( \frac{W^\prime(\alpha(p), S_M)}{NW(H_{[1,a]_p}) + \sum_{ij} W_{ij}^{(1)} + W_{ij}^{(2)} + 2W_{ij}^{(3)} }, 1 \right) \\
        =& \mathrm{min}\left(\frac{M-n+1}{\beta \left(N\Omega - 4N_b\delta_b + \sum_{ij} V_{ij} + 4C_{ij}\right) },1\right).
    \end{align*}
    \item If the unity operator $H_{0,0}$ is encountered, decide whether or not to accept inserting a diagonal operator ($n \rightarrow n + 1$) with the probability 
    \begin{align}
        P_{H_{[0,0]_p} \rightarrow H_{[1,a]_p} / H_{[1,b]_p}} =& \mathrm{min}\left(\frac{\beta \left(N\Omega - 4N_b\delta_b + \sum_{ij} V_{ij} + 4C_{ij}\right) }{M-n},1\right)
    \end{align}
    \item If it was decided to insert a diagonal operator, choose the operator to insert with probabilities
    \begin{subequations} \label{eq:heat_bath_1_ryd}
        \begin{align} 
            P_{H_{[1,a]_p}} = & \frac{N\Omega}{N\Omega - 4N_b\delta_b + \sum_{ij} V_{ij} + 4C_{ij}}, \\
            P_{H_{[1,b]_p}} = & \frac{ - 4N_b\delta_b + \sum_{ij} V_{ij} + 4C_{ij}}{N\Omega - 4N_b\delta_b + \sum_{ij} V_{ij} + 4C_{ij}}.
        \end{align}
    \end{subequations}
    \item To choose which type of diagonal operator (site or bond), we follow the algorithm outlined in Sandvik's paper on an SSE framework for the TFIM with arbitrary interactions. 
    We wish to insert an operator at bond $(i,j)$ (note that a site operator is a bond operator when $i=j$). 
    Let the matrix $M$ have elements $M_{ij}$ such that
    \begin{align} \label{eq:choices}
        M_{ij} = \begin{cases}
            \Omega & \text{ for } i = j \\
             W^{(1,2,3)}_{ij} & \text{ for } i \neq j
        \end{cases}
    \end{align}
    We then store the cumulative probabilities
    \begin{align}
        P_c(k=1,...,N) =& \frac{\sum_{i=1}^k P(i)}{\sum_{i=1}^N P(i)}
    \end{align}
    where $P(i) = \sum_j M_{ij}$.
    Select the index $i$ in $(i,j)$ by choosing a random number $R$ between 0 and 1 and search $P_c$ for the smallest $k$ such that $P(k) = \sum_j M_{kj} \geq R$. 
    Then, the first index of $(i,j)$ is $i = k$.
    The index $j$ is then chosen probabilistically by noting that the relative probability for $j$ given $i=k$ is $M_{kj}$.
    So, given this choice of $(i,j)$, if $i = j$ we insert $H_{1,a}$. 
    If $i \neq j$, we insert $H_{1,b}$.
\end{enumerate}

\subsubsection{Zero Temperature}

Since we want to maximize the expansion order here to project out the ground state of the Hamiltonian, we want to avoid using identity operators and force our simulation cell to be filled with local terms of our Hamiltonian.
The steps for the diagonal update here are therefore the exact same as with the diagonal updates for finite-temperature with the only difference being that we {\it always} remove the current diagonal operator and keep attempting to replace it with a new one.

\subsubsection{Deterministic loops}

Instead of the forseeably inefficient non-deterministic directed loops, we may instead form deterministic clusters in the same exact way as is outlined in the well-known SSE formalism for TFIM.
\begin{itemize}
    \item Clusters terminate on site operators.
    \item When a bond operator is encountered, the cluster branches from every possible leg of the bond vertex.
\end{itemize}
However, the difference here is that instead of flipping with probability 1/2, we now flip each cluster, labelled by an index $c$, in the simulation cell with the following Metropolis probability.
\begin{align}
    P_{\text{flip},c} =& \min\left(\frac{W_c^\prime}{W_c}, 1\right)
\end{align}
Here, $W_c$ primed (unprimed) denotes the weight of the local cluster $c$ being flipped (unflipped).
Specifically, $W_c$ (primed or unprimed) is the product of all matrix element weights ($W_1$, $W_2$, $W_3$) present in the local cluster $c$. 
Note that the site operators' matrix elements do not factor into $P_{\text{flip},c}$ since they are symmetric upon flipping (there is no weight change). 

\section{Lattices}

Here, I will just give a quick description of the lattice types that can be found in \texttt{src/lattice\_types.jl}.
Currently we have 1D chains, rectangular, triangular, and ruby lattices available, all of which contain periodic boundary condition options except for the ruby lattice. 
I've taken the site enumeration convention for the lattice types as follows (see Fig.~\ref{fig:2Dlattices}).

\begin{enumerate}
    \item 1D chains: site labels increase from left to right.
    \item 2D lattices with one site per unit cell: site labels increase like periodic boundaries in the x-direction.
    \item 2D lattices with a multisite basis: site labels within a unit cell increase left to right, then top to bottom. 
\end{enumerate}

\begin{figure}
    \centering
    \includegraphics[width=0.3\linewidth]{2Dsimple.png}
    \includegraphics[width=0.5\linewidth]{2Dmulti.png}
    \caption{The site enumeration convention for 2D simple lattices (left) and 2D lattices with a many-site basis (right).}
    \label{fig:2Dlattices}
\end{figure}
\end{document}