#  Spontaneous charge current in time-reversal-symmetry breaking phase in
kagome metals

Koki Shimura1, Rina Tazai2, Youichi Yamakawa1, Seiichiro Onari1, and Hiroshi
Kontani1 1Department of Physics, Nagoya University, Furo-cho, Nagoya 464-8602,
Japan
2 Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502,
Japan.

###### Abstract

The charge loop current (cLC) state has attracted increasing attention in
kagome metals. Here, we calculate the spontaneous currents along the nearest
sites $i$ and $j$, $J_{i,j}$, induced by the cLC order that is the imaginary
and odd-parity hopping integral modulation $\delta t_{i,j}$. We reveal that
the magnitude of $J_{i,j}$ strongly depends on the nearest sites $i$ and $j$
in the $2\times 2$ cLC state, where $\eta\equiv|\delta t_{i,j}|$ is equivalent
for all nearest sites. The obtained $J_{i,j}$ becomes large near the van-Hove
singularity (vHS) filling ($n\sim n_{\rm vHS}$) even when $\eta$ is fixed.
Interestingly, the obtained $J_{i,j}$ exhibits the logarithmic divergence
behavior at low temperatures for $n\sim n_{\rm vHS}$ with a fixed $\eta$, by
reflecting the vHS points that are the characteristic of kagome metals. The
present study provides useful information for local electronic state
measurements, such as the site-selective NMR and STM experiments.

Introduction \- In recent years, kagome lattice metals have become one of
central topics in the field of strongly correlated electron systems. Multiple
exotic metallic phase transitions, such as the star-of-David bond-order (BO),
superconducting (SC) state, and time-reversal-symmetry-breaking (TRSB) phase
without spin polarization, are realized thanks to the strong geometrical
frustration that prevents simple spin and charge orders. In $A$V3Sb5
($A$=Cs,Rb,K), the $2\times 2$ BO appears at $T_{\rm BO}\approx 100$K [1, 2].
Inside the BO phase, the SC state appears at $T_{\rm SC}\approx 1-10$K [3, 4,
5, 6, 7]. In addition, the electronic nematic state has been observed inside
[8, 9, 10] and outside [11] of the BO phase. To explain the observed
nematicity, the single-${\bm{Q}}$ odd-parity BO [12] and $E_{g}$ symmetry
nematic BO [13] have been proposed theoretically.

To understand these exotic electronic states in kagome metals, various
theoretical studies have been performed, based on the extended mean-field
theory [14, 15], the functional renormalization group (RG) theory [16, 17],
the parquet RG theory [19, 18], and the density-wave (DW) equation analysis
with the higher-order vertex corrections (VCs) [20]. Both the RG theory and
the DW equation theory have been successfully applied to Fe-based [21, 22, 23,
38, 39, 43] and cuprate superconductors [44]. The nematic order and the BO in
these systems originate from the higher-order VCs [38]. Recently, the present
authors proposed a unified explanation for the BO and the heavily anisotropic
$s$ wave SC state in kagome metals based on the DW equation analysis [20].
However, theoretical understanding for the TRSB state has been limited until
recently.

Recently, the TRSB state has been observed by various experimental methods.
Inside the BO phase, the TRSB has been observed by the STM [1], the $\mu$SR
[27, 28, 26], the anomalous Hall effect [24, 25], and the magnetochiral
anisotropy (eMChA) [45, 29] measurements. The spontaneous charge loop-current
(cLC) order is a natural candidate of the TRSB order. In kagome metals, the
$2\times 2$ cLC order has finite orbital magnetization $M_{\rm orb}$; see Ref.
[12]. References [45, 29] report that the TRSB domains with random chiralities
are detwinned by the magnetic field $h_{z}\sim 1{\rm T}$, and $T_{\rm TRSB}$
is drastically enhanced by small $h_{z}$ as well as the uniaxial pressure
${\epsilon}$. Such drastic $h_{z}$\- and ${\epsilon}$-dependences of the TRSB
states are naturally explained based on the Ginzburg-Landau (GL) theory under
$(h_{z},{\epsilon})$ in Ref. [12], owing to the finite $M_{\rm orb}$. Very
interestingly, the TRSB state outside the BO phase ($T_{\rm TRSB}\approx
130$K) has been discovered by the recent magnetic torque measurement [11]. The
single ${\bm{Q}}$ cLC order ($M_{\rm orb}=0$) is a natural candidate.
Actually, the cLC order can emerge above $T_{\rm BO}$ theoretically [35].

The cLC order was originally studied in cuprate superconductors [31, 32, 33].
The cLC order is defined as the imaginary and odd-parity hopping integral
modulation between sites $i$ and $j$, $\delta t_{i,j}$ [20, 35, 36]. The Berry
curvature due to $\delta t_{i,j}$ gives rise to the spontaneous current [30].
To understand the microscopic origin of the cLC order in kagome metals,
various theoretical studies have been performed, such as the mean-field theory
[14, 15], the parquet RG theory [19, 18], and the DW equation analysis [35].
Importantly, the BO fluctuations in kagome metals mediate not only the SC
pairing, but also the cLC order that is the TRSB particle-hole pair
condensation; see Ref. [35]. Therefore, the cLC state is naturally expected to
occur in kagome metals with the BO instability. Recently, the competition
between the BO, cLC and the SC states is intensively studied based on the GL
theories [19, 12, 40, 47].

In the cLC state in kagome metal model, the uniform orbital magnetization
$M_{\rm orb}$ was recently studied in Ref. [12]. However, the nanoscale charge
current flowing between the nearest bonds, $J_{i,j}$, has not been performed
to our knowledge. The knowledge of $J_{i,j}$ is significant to understand the
local magnetization measurements, such as the $\mu$SR [27, 28, 26] and the NMR
[46] measurements.

In this Letter, we analyze the spontaneous current $J_{i,j}$ in the $2\times
2$ cLC order phase. We reveal that the magnitude of $J_{i,j}$ strongly depends
on the nearest sites $(i,j)$. The obtained $J_{i,j}$ becomes large near the
van-Hove singularity (vHS) filling ($n\sim n_{\rm vHS}$). Interestingly,
$J_{i,j}$ exhibits the logarithmic divergence behavior at low temperatures for
$n\sim n_{\rm vHS}$, by reflecting the vHS points that are the characteristic
of kagome metals. The present study provides useful information for local
electronic state measurements, such as the site-selective NMR and STM
experiments.

Model Hamiltonian with $3\bm{Q}$ cLC order \- First, we show a 3-site kagome
lattice model in Fig. 1 (a). Here, the nearest-neighbor hopping integral is
$t=-0.5$[eV]. Additionally, to prevent perfect nesting, a next-nearest-
neighbor hopping integral is introduced with $t^{\prime}=-0.02$[eV]. There are
three sites, A, B, and C in the unit cell. The vectors between nearest-
neighboring sites are defined as $\bm{a}_{AB},\bm{a}_{BC}$, and $\bm{a}_{CA}$.
Fig. 1 (b) shows the Fermi surface (FS) of the 3-site kagome lattice model and
nesting vectors ${{\bm{q}}}_{n}$ $(n=1,2,3)$. There exist three
${{\bm{q}}}_{n}$ that connect vHS at the M point.

Figure 1: (a) Tight-binding model of the kagome lattice model. (b) FS of the
3-site kagome model at $n=5.0$. (c) The $3\bm{Q}$ cLC
$\bm{\eta}=(\eta,\eta,\eta)/\sqrt{3}$ depicted in the kagome lattice.

When the nonzero wave vectors ${{\bm{q}}}_{n}$ of BO or cLC emerge in kagome
lattice, the translational symmetry of the system is violated. In particular,
by $3\bm{Q}$ order that three ${{\bm{q}}}_{n}$ coexist, the unit cell expands
to 12 sites. Fig. 1 (c) shows the $3\bm{Q}$ current order
$\bm{\eta}=(\eta_{1},\eta_{2},\eta_{3})$. We assign $\delta
t_{i,j}^{c}=i\eta_{1}$ $[\delta t_{i,j}^{c}=-i\eta_{1}]$ for
$(l,m)=(1,2),(2,4),(4,5),(5,1)$ $[(l,m)=(7,8),(8,10),(10,11),(11,7)]$, $\delta
t_{i,j}^{c}=i\eta_{2}$ $[\delta t_{i,j}^{c}=-i\eta_{2}]$ for
$(l,m)=(2,3),(3,11),(11,12),(12,2)$$[(l,m)=(5,6),(6,8),(8,9),(9,5)]$, and
$\delta t_{i,j}^{c}=i\eta_{3}\ [\delta t_{i,j}^{c}=-i\eta_{3}]$ for
$(l,m)=(4,6),(6,10),(10,12),(12,4)$$[(l,m)=(1,3),(3,7),(7,9),(9,1)]$. Here
$\delta t^{c}_{i,j}$ is the purely imaginary ($(\delta
t_{i,j}^{c})^{*}=-\delta t_{i,j}^{c}$) symmetry breaking term of the hopping
integral, which is odd parity ($\delta t^{c}_{i,j}=-\delta t^{c}_{j,i}$) [35].

The Hamiltonian with the current order is given by

$\hat{H}=\sum_{{{\bm{k}}},l,m}h_{lm}({{\bm{k}}})c^{\dagger}_{{{\bm{k}}},l}c_{{{\bm{k}}},m},$
(1)

where $h_{lm}({{\bm{k}}})$ is the Fourier transformation of
$t_{i,j}=t^{(0)}_{i,j}+\delta t^{c}_{i,j}$. The hopping integral is
constructed with the original hopping integral $t^{(0)}_{i,j}$ and the current
order $\delta t^{c}_{i,j}$. $l$ and $m$ denote the site number within the unit
cell, $l,m=1\sim 12$. We use a index $i=\\{\alpha,l\\}$, where $\alpha$
denotes the unit cell. We introduce the notation
${{\bm{R}}}_{i}={{\bm{R}}}_{\alpha}+{{\bm{r}}}_{l}$ to designate the position
of atoms, where ${{\bm{R}}}_{\alpha}$ is the position of the unit cell and
${{\bm{r}}}_{l}$ is the relative position of the atom within the unit cell. In
the 12-site kagome lattice model, the FS folds, and the M and $\Gamma$ points
in the lattice become equivalent. Therefore, the wave vector of the cLC order
is ${{\bm{q}}}_{n}=\bm{0}$ for $n=1,2,3$ within folded BZ.

Here, we introduce the form factor of the cLC, which is given by the Fourier
transformation of $\delta t^{c}_{i,j}$ shown in Fig. 1 (c). In the 3-site
kagome lattice model, the cLC form factor between the nearest site $S$ and
site $S^{\prime}$ ($S,S^{\prime}=A,B,C$) is given as
blue$f_{SS^{\prime}}^{{{\bm{q}}}_{n}}({{\bm{k}}})=2\sin({{\bm{k}}}\cdot\bm{a}_{SS^{\prime}})$,
where inter-site vector $\bm{a}_{SS^{\prime}}$ and the wavevector
${{\bm{q}}}_{n}$ are shown in Fig. 1 (a) and (b), respectively. Here,
$n=1,2,3$ for
$\left\\{S,S^{\prime}\right\\}=\left\\{A,B\right\\},\left\\{B,C\right\\},\left\\{C,A\right\\}$,
respectively. The set of the current order functions is given by
$(\eta_{1}f_{AB},\eta_{2}f_{BC},\eta_{3}f_{CA})$, where $\eta_{n}$ is the nth
cLC order parameter. A more detailed explanation is presented in Ref. [20].
Hereafter, we set $\eta_{1}=\eta_{2}=\eta_{3}=\eta$ in the numerical study. We
can also present the cLC form factor in the 12-site model,
${\tilde{f}}_{lm}(k)$, where $l,m=1\sim 12$ [12]. For example,
${\tilde{f}}_{12}({{\bm{k}}})=ie^{-i\bm{k}\cdot\bm{a}_{12}}$. In this
representation, the cLC order wavevector becomes uniform
$({{\bm{q}}}=\bm{0})$, so the incoming and outgoing momenta of
${\tilde{f}}_{lm}({{\bm{k}}})$ are the same. Here, we introduce the $12\times
12$ matrix expression of the current order function
$[{\hat{C}}_{{{\bm{k}}}}]_{l,m}\equiv{\tilde{f}}_{l,m}({{\bm{k}}})$, which is
convenient in the present numerical study. Here,
$[{\hat{C}}_{{{\bm{k}}}}]_{l,m}\neq 0$ only when $(l,m)$ is the pair of the
nearest sublattice.

The $3\bm{Q}$-cLC form factor of kagome metals is derived microscopically from
the analysis of the DW equation with the BO fluctuation exchange term as the
kernel, as shown in ref. [35]. The $\delta t_{i,j}$ derived in ref. [35] has
the same absolute value for all nearest-neighbor bonds $(i,j)$, meaning
$|\delta t_{i,j}^{c}|=\eta$.  It is worth noting that the $\delta t_{i,j}$
derived in ref. [35] also has long-range components, but for simplicity in
this paper, we consider only the nearest-neighbors.

Figure 2: Band structure, DOS, and FS of the 12-site kagome lattice model with
$\eta=0.1$.

We show the band structure, density of states (DOS), and FS of the 12-site
kagome lattice model that cLC order is introduced at $\eta=0.1$ in Fig 2. In
the 12-site kagome lattice, vHS exists at the M-point, leading to a large DOS.
The introduction of cLC order resolves the strong degeneracy near the M point,
while the large magnitude of the DOS remains.

Next, to calculate the spontaneous charge current in the cLC phase, we
introduce the current order operator from Heisenberg eq.
$\dot{n}_{i}=-i[n_{i},H]=-i\left(\sum_{j}t_{i,j}c_{i}^{\dagger}c_{j}-\sum_{j}t_{j,i}c_{j}^{\dagger}c_{i}\right)$.
Then, the current operator is derived from the continuity equation as

$\displaystyle j_{i,j}=-i(t_{i,j}c^{\dagger}_{i}c_{j}-(i\leftrightarrow j)).$
(2)

By taking the expectation value with respect to grand canonical ensemble, the
current from site $j$ to site $i$ is given as

$\displaystyle J_{i,j}=\langle
j_{i,j}\rangle=it_{i,j}g_{j,i}-(i\leftrightarrow j),$ (3)

where the equal time Green function is given as
$g_{i,j}=\lim_{u\rightarrow-0}\langle Tc^{\dagger}_{i}(u)c_{j}(0)\rangle$.
$\langle\cdots\rangle$ represents the mean value for the grand canonical
ensemble. The Green function in momentum representation is

$\displaystyle
G_{l,m}({{\bm{k}}},i\epsilon_{n})=((i\epsilon_{n}+\mu)\hat{\bm{1}}-\hat{h}({{\bm{k}}}))^{-1}_{l,m}.$
(4)

We derive the equal time Green function by the Fourier transform of Eq. (4) as
$g_{i,j}=T\sum_{n,{{\bm{k}}}}G_{l,m}({{\bm{k}}},i\epsilon_{n})e^{i{{\bm{k}}}\cdot({{\bm{R}}}_{\alpha}-{{\bm{R}}}_{\beta})}e^{i\epsilon_{n}\delta}$.
Note that the factor $e^{i\epsilon_{n}\delta}$ is not necessary for $l\neq m$.
$\epsilon_{n}=(2n-1)\pi T$ is the fermion Matsubara frequency. Here n is
integer.

Figure 3:  (a) Obtained distribution of the charge currents. Here,
$J(J^{\prime\prime})$ is the current along the triangular (hexagonal) loop
where the direction of each cLC order is aligned. The current on other bonds
is given as $J^{\prime\prime}$ (b) Obtained currents as a functions of $\eta$.
(c) 3-site cluster models with different cLC orders. The circulating current
in $\rm(\hskip 1.79993pti\hskip 1.79993pt)$ differs from that in $\rm(\hskip
0.80002ptii\hskip 0.80002pt)$; see the text.

Analysis of spontaneous current \- Now, we present the results of the current
distribution induced by the cLC orders. Fig. 3 (a) illustrates the real-space
current distribution within the kagome lattice. The obtained currents as a
function of $\eta$ are shown in Fig. 3 (b) ($T=2.5$meV, $n=5.0$) . As $\eta$
increases, the difference in the magnitude of the currents is enlarged.

To understand the obtained nontrivial relation between $J_{i,j}$ and $\delta
t_{i,j}$, here we analyze a simple 3-site cluster model with cLC order
oriented in different directions. Fig. 3 (c) represents two 3-site cluster
models, $\rm(\hskip 1.79993pti\hskip 1.79993pt)$ and $\rm(\hskip
0.80002ptii\hskip 0.80002pt)$. In the case of $\rm(\hskip 1.79993pti\hskip
1.79993pt)$, $\delta t^{c}_{ij}$ is given as $\delta t^{c}_{21}=\delta
t^{c}_{31}=\delta t^{c}_{32}=i\eta$. The current $J^{1}_{1,2}$ that flows in
model $\rm(\hskip 1.79993pti\hskip 1.79993pt)$ is derived as
$J^{1}_{1,2}=-2\eta(t^{2}+\eta^{2})A$, where
$A=T\sum_{n}\left[(i\epsilon_{n}+\mu)^{3}+3(t^{2}+\eta^{2})(i\epsilon_{n}+\mu)+2t(t^{2}+\eta^{2})\right]^{-1}$.
In the case of $\rm(\hskip 1.79993ptii\hskip 1.79993pt)$, $\delta t^{c}_{ij}$
is given as $\delta t^{c}_{12}=\delta t^{c}_{23}=\delta t^{c}_{31}=i\eta$.
Then, the current $J^{2}_{1,2}$ is caluculated as
$J^{2}_{1,2}=6\eta(t^{2}-\frac{\eta^{2}}{3})B$, where
$B=T\sum_{n}\left[(i\epsilon_{n}+\mu)^{3}+3(t^{2}+\eta^{2})(i\epsilon_{n}+\mu)\right.$$\left.+2t(t^{2}-3\eta^{2})\right]^{-1}$.
Thus, different magnitudes of current emerge even when the same order
parameter $\eta$ is provided in three directions
${{\bm{q}}}_{1},{{\bm{q}}}_{2},{{\bm{q}}}_{3}$. It holds that
$|J_{12}^{m}|=|J_{23}^{m}|=|J_{31}^{m}|$ for each m = 1,2. In the case that
$|a|\ll|t|$, $|J_{12}^{2}|=3|J_{12}^{1}|$ holds. In this finite site model,
the current converges to zero at $T=0$ and $T\rightarrow\infty$.

The mathematical reason for the difference in current between $\rm(\hskip
1.79993pti\hskip 1.79993pt)$ and $\rm(\hskip 1.79993ptii\hskip 1.79993pt)$ is
that the $J_{1,2}$ given in Eq. (3) contributes not only to $\delta
t_{1,2}^{c}$ but also to $\delta t_{1,3}^{c}$, $\delta t_{2,3}^{c}$, which are
included in $g_{1,2}$. When the phases of $\delta t_{i,j}^{c}$ on the closed
loop are aligned, i.e., $\delta t_{1,2}^{c}=\delta t_{2,3}^{c}=\delta
t_{3,1}^{c}$, the induced current $|J_{i,j}|$ becomes the largest. This
consideration leads us to understand why $J$ ,$J^{\prime}$, and
$J^{\prime\prime}$ are different in kagome metals as shown in Fig. 3(a)-(b).

Figure 4:  (a) $n$-dependence of the currents at $T=2.5-10.0$ meV. (b)
T-dependence of the currents at $n=5.0-5.2$.

Next, we analyze the filling and temperature dependence of the current. The
filling dependence of the current is depicted in Fig. 4 (a). It shows that the
current reaches a maximum near the vHS filling. Fig. 4 (b) illustrates the
temperature dependence of the current. The magnitude of the current increases
at the low-tempratures, with a pronounced increase near the vHS filling.

Susceptibility of the current \- Next stage, we consider susceptibility that
describes the behavior of the current. First, we expand the Green function in
terms of the order parameter $\eta$ as
$\hat{G}\left({{\bm{k}}},i\epsilon_{n}\right)\sim\hat{G}^{(0)}\left({{\bm{k}}},i\epsilon_{n}\right)+\hat{G}^{(0)}\left({{\bm{k}}},i\epsilon_{n}\right)\eta{\hat{C}}_{{{\bm{k}}}}\hat{G}^{(0)}\left({{\bm{k}}},i\epsilon_{n}\right)$,
where $\hat{G}^{(0)}({{\bm{k}}},i\epsilon_{n})$ is the Green function with
$\eta=0$. The contribution from the term
$\hat{G}^{(0)}({{\bm{k}}},i\epsilon_{n})\eta{\hat{C}}_{{{\bm{k}}}}\hat{G}^{(0)}({{\bm{k}}},i\epsilon_{n})$
to the current is dominant. The diagramatic representations of this term is
shown in Fig. 5 (a). Meanwhile, the static irreducible susceptibility is

$\displaystyle\hat{\chi}^{J-\eta}_{lml^{\prime}m^{\prime}}({{\bm{q}}})$
$\displaystyle=-\frac{T}{N}\sum_{{{\bm{k}}},n}G^{(0)}_{ll^{\prime}}({{\bm{k}}}+{{\bm{q}}},i\epsilon_{n})G^{(0)}_{m^{\prime}m}({{\bm{k}}},i\epsilon_{n}).$
(5)

As an example of form factor,
$[\hat{C}_{{{{\bm{k}}}}}]_{12}=ie^{i{{\bm{k}}}\cdot\bm{a}_{AB}}$, and
$[\hat{C}_{{{\bm{k}}}}]_{21}=-ie^{i{{\bm{k}}}\cdot\bm{a}_{BA}}$. Thus, at
${{\bm{k}}}=0$ which corresponds to $\Gamma$ point,
$[\hat{C}_{\bm{0}}]_{12}=i$ holds. Here, we recall that sites $1$,$4$,$7$, and
$10$ belong to the A site, while site $2$,$5$,$8$, and $11$ belong to B site.
In the case of the 12-site kagome lattice model, the form factor exhibits the
largest value at the $\Gamma$ point. Therefore, after approximating it as
$[\hat{C}_{{{\bm{k}}}}]_{lm}\sim[\hat{C}_{\bm{0}}]_{lm}$, it is possible to
perform the following transformations as

$\displaystyle J_{1,2}$ $\displaystyle\propto
T\sum_{n,l,m,{{\bm{k}}}}[G^{(0)}]_{1,l}\eta[\hat{C}_{\bm{0}}]_{l,m}[G^{(0)}]_{m,2}e^{i{{\bm{k}}}\cdot({{\bm{R}}}_{\alpha}-{{\bm{R}}}_{\beta})}$
(6)
$\displaystyle\sim\sum_{l,m=1-12}\chi_{12,lm}^{J-\eta}(\bm{0})[\hat{C}_{\bm{0}}]_{lm}\eta.$

Importantly, $\chi_{12,lm}^{J-\eta}(\bm{0})$ becomes large when
$\left\\{l,m\right\\}=\left\\{1,2\right\\},\left\\{4,5\right\\},\left\\{7,8\right\\},\left\\{10,11\right\\}$,
by reflecting the sublattice interference of the pure-type FS. If one puts a
restriction $\left\\{l,m\right\\}=\left\\{1,2\right\\}$, $J_{1,2}$ in Eq. (6)
is proportional to $\chi_{12,12}^{J-\eta}-\chi_{12,21}^{J-\eta}$.

Figure 5: (a) Diagrammatic representation of the first-order term
$\left(\hat{G}^{(0)}\eta\hat{C}\hat{G}^{(0)}\right)_{m,l}$. (b) Diagrammatic
representation of the sum of susceptibility
$\chi_{lml^{\prime}m^{\prime}}^{J-\eta}-\chi_{lmm^{\prime}l^{\prime}}^{J-\eta}$.
(c) Filling dependence of $\phi_{12}^{J-\eta}$ at $T=5.0,10,30$ meV. (d)
Temperature dependence of $\phi_{12}^{J-\eta}$ at $n=5.0-5.2$.

Figure 5 (b) represents diagrammatic representations of
$\chi_{1l,m2}^{J-\eta}(\bm{0})$. Taking account of the form factors depicted
in Fig. 1 (c), the summation of contributions from
$\left\\{l,m\right\\}=\left\\{1,2\right\\}$, $\left\\{4,5\right\\}$,
$\left\\{7,8\right\\}$, and $\left\\{10,11\right\\}$ leads to
$\phi_{12}^{J-\eta}=[\chi_{12,12}^{J-\eta}-\chi_{12,21}^{J-\eta}+\chi_{12,42}^{J-\eta}-\chi_{12,24}^{J-\eta}+\chi_{12,45}^{J-\eta}-\chi_{12,54}^{J-\eta}+\chi_{12,14}^{J-\eta}-\chi_{12,41}^{J-\eta}]-[\chi_{12,78}^{J-\eta}-\chi_{12,87}^{J-\eta}+\chi_{12,10\
8}^{J-\eta}-\chi_{12,8\ 10}^{J-\eta}+\chi_{12,10\ 11}^{J-\eta}-\chi_{12,11\
10}^{J-\eta}+\chi_{12,7\ 10}^{J-\eta}-\chi_{12,10\ 7}^{J-\eta}]$.

The filling and temperature dependence of $\phi_{12}^{J-\eta}$ are shown
respectively in Fig. 5 (c) and (d). $\phi_{12}^{J-\eta}$ reaches its maximum
value at the vHS filling, with a notable increase in the low temperature
regime. In kagome metals, $\phi_{12}^{J-\eta}$ significantly increases at low
temperatures for $n\sim n_{VHS}$. Analytically, at vHS filling,
$\phi_{12}^{J-\eta}\propto 1/E_{0}\left(\ln E_{0}/T\right)^{2}+\chi_{0c}$ when
$E_{0}\sim 0.1$[eV] is the energy range where the vHS bandstructure is well
described by the quadratic expression, as discussed in Ref. [19]. Here,
$\chi_{0c}$ is the constraint susceptibility, which excludes the contribution
for $|E|<E_{0}$ from Eq. (5), which is defined in SM in Ref. [37]. This
temperature dependence is attributed to the singularity at vHS in the kagome
lattice model, in stark contrast to usual metals, where $\phi_{12}^{J-\eta}$
remains constant.

Effect of self-energy \- In addition, we analyze how the current is modified
by the self-energy. Fig. 6 (a) shows the filling dependence of the current in
the presence of self-energy. Fig. 6 (b) shows the obtained damping rate
$\gamma$ dependence of current. As $\gamma$ increases, the current is
gradually suppressed.

Figure 6:  (a) Filling dependence of the current $J^{\rm
min}=\min\\{J,J^{\prime},J^{\prime\prime}\\}$ for various $\gamma$. (b)
$\gamma$ dependence of $J^{\rm min}$ for $n=5.0-5.2$.

Summary and discussions \- In this letter, we employed the Green function
method to calculate the real-space current distribution in the $3{\bm{Q}}$-cLC
phase of kagome metal. The obtained magnitude of $J_{i,j}$ strongly depends on
the nearest sites ($i,j$) even when $|\delta t_{i,j}|$ is constant for any
bonds. The obtained $J_{i,j}$ becomes large near the vHS filling, while it is
suppressed when considering self-energy. Interestingly, $J_{i,j}$ exhibits the
logarithmic divergence behavior at low temperatures for $n\sim n_{\rm vHS}$.
The present study provides useful information for local electronic state
measurements, such as the site-selective NMR and STM experiments.

We have employed the Green function method to calculate $J_{i,j}$. Based on
this method, one can calculate the current by including the beyond-mean-field
electron correlations, such as the self-energy and the VCs. This is a great
advantage of the present method.

###### Acknowledgements.

This study has been supported by Grants-in-Aid for Scientific Research from
MEXT of Japan (JP20K03858, JP20K22328, JP22K14003, JP23K03299), and by the
Quantum Liquid Crystal No. JP19H05825 KAKENHI on Innovative Areas from JSPS of
Japan.

## References

  * [1] Yu-Xiao Jiang, Jia-Xin Yin, M. Michael Denner, Nana Shumiya, Brenden R. Ortiz, Gang Xu, Zurab Guguchia, Junyi He, Md Shafayat Hossain, Xiaoxiong Liu, Jacob Ruff, Linus Kautzsch, Songtian S. Zhang, Guoqing Chang, Ilya Belopolski, Qi Zhang, Tyler A. Cochran, Daniel Multer, Maksim Litskevich, Zi-Jia Cheng, Xian P. Yang, Ziqiang Wang, Ronny Thomale, Titus Neupert, Stephen D. Wilson, Zahid Hasan Unconventional chiral charge order in kagome superconductor KV3Sb5, Nat. Mat. 20, 1353-1357 (2021).
  * [2] Hong Li, He Zhao, Brenden R. Ortiz, Takamori Park, Mengxing Ye, Leon Balents, Ziqiang Wang, Stephen D. Wilson, Ilija Zeljkovic Rotational symmetry breaking in the normal state of a kagome superconductor KV3Sb5, Nat. Phys. 18, 265 (2022)
  * [3] B. R. Ortiz, L. C. Gomes, J. R. Morey, M. Winiarski, M. Bordelon, J. S. Mangum, I. W. H. Oswald, J. A. Rodriguez-Rivera, J. R. Neilson, S. D. Wilson, E. Ertekin, T. M. McQueen, and E. S. Toberer, New kagome prototype materials: discovery of ${\mathrm{KV}}_{3}{\mathrm{Sb}}_{5},{\mathrm{RbV}}_{3}{\mathrm{Sb}}_{5}$, and ${\mathrm{CsV}}_{3}{\mathrm{Sb}}_{5}$, Phys. Rev. Materials 3, 094407 (2019).
  * [4] B. R. Ortiz, S. M. L. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, R. Seshadri, L. Balents, J. He, and S. D. Wilson, $\mathrm{Cs}{\mathrm{V}}_{3}{\mathrm{Sb}}_{5}$: A ${\mathbb{Z}}_{2}$ Topological Kagome Metal with a Superconducting Ground State, Phys. Rev. Lett. 125, 247002 (2020).
  * [5] F. H. Yu, D. H. Ma, W. Z. Zhuo, S. Q. Liu, X. K. Wen, B. Lei, J. J. Ying, and X. H. Chen, Unusual competition of superconductivity and charge-density-wave state in a compressed topological kagome metal, Nat. Commun. 12, 3645 (2021).
  * [6] K. Y. Chen, N. N. Wang, Q. W. Yin, Y. H. Gu, K. Jiang, Z. J. Tu, C. S. Gong, Y. Uwatoko, J. P. Sun, H. C. Lei, J. P. Hu, and J.-G. Cheng, Double Superconducting Dome and Triple Enhancement of ${T}_{c}$ in the Kagome Superconductor ${{CsV}}_{3}{{Sb}}_{5}$ under High Pressure, Phys. Rev. Lett. 126, 247001 (2021).
  * [7] B. R. Ortiz, P. M. Sarte, E. M. Kenney, M. J. Graf, S. M. L. Teicher, R. Seshadri, and S. D. Wilson, Superconductivity in the ${\mathbb{Z}}_{2}$ kagome metal ${{KV}}_{3}{{Sb}}_{5}$, Phys. Rev. Materials 5, 034801 (2021).
  * [8] L. Nie, K. Sun, W. Ma, D. Song, L. Zheng, Z. Liang, P. Wu, F. Yu, J. Li, M. Shan, D. Zhao, S. Li, B. Kang, Z. Wu, Y. Zhou, K. Liu, Z. Xiang, J. Ying, Z. Wang, T. Wu, and X. Chen, Charge-density-wave-driven electronic nematicity in a kagome superconductor, Nature 604, 59 (2022).
  * [9] Y. Xu, Z. Ni, Y. Liu, B. R. Ortiz, S. D. Wilson, B. Yan, L. Balents, and L. Wu, Universal three-state nematicity and magneto-optical Kerr effect in the charge density waves in AV3Sb5 (A=Cs, Rb, K), arXiv:2204.10116.
  * [10] H. Li, H. Zhao, B. R. Ortiz, T. Park, M. Ye, L. Balents, Z. Wang, S. D. Wilson, and I. Zeljkovic, Rotation symmetry breaking in the normal state of a kagome superconductor KV3Sb5, Nat. Phys. 18, 265 (2022).
  * [11] T. Asaba, A. Onishi, Y. Kageyama, T. Kiyosue, K. Ohtsuka, S. Suetsugu, Y. Kohsaka, T. Gaggl, Y. Kasahara, H. Murayama, K. Hashimoto, R. Tazai, H. Kontani, B. R. Ortiz, S. D. Wilson, Q. Li, H. -H. Wen, T. Shibauchi, Y. Matsuda, Evidence for an odd-parity nematic phase above the charge density wave transition in kagome metal CsV3Sb5. accepted for publication in Nat. Phys. (arXiv:2309.16985).
  * [12] R. Tazai, Y. Yamakawa, and H. Kontani, Drastic magnetic-field-induced chiral current order and emergent current-bond-field interplay in kagome metals, accepted for publication in Proceedings of the National Academy of Sciences (PNAS) (available at https://arxiv.org/abs/2303.00623).
  * [13] arXiv:2305.18087
  * [14] M. M. Denner, R. Thomale, and T. Neupert, Analysis of Charge Order in the Kagome Metal $A{{V}}_{3}{{Sb}}_{5}$ ($A={K},{Rb},{Cs}$), Phys. Rev. Lett. 127, 217601 (2021).
  * [15] Y.-P. Lin and R. M. Nandkishore, Complex charge density waves at Van Hove singularity on hexagonal lattices: Haldane-model phase diagram and potential realization in the kagome metals AV3Sb5 (A = K, Rb, Cs), Phys. Rev. B 104, 045122 (2021).
  * [16] M. L. Kiesel, C. Platt, and R. Thomale, Unconventional Fermi Surface Instabilities in the Kagome Hubbard Model, Phys. Rev. Lett. 110, 126405 (2013).
  * [17] W.-S. Wang, Z.-Z. Li, Y.-Y. Xiang, and Q.-H. Wang, Competing electronic orders on kagome lattices at van Hove filling, Phys. Rev. B 87, 115135 (2013).
  * [18] H. D. Scammell, J. Ingham, T. Li, and O. P. Sushkov, Chiral excitonic order from twofold van Hove singularities in kagome metals, Nat. Commun. 14, 605 (2023)
  * [19] T. Park, M. Ye, and L. Balents, Electronic instabilities of kagome metals: Saddle points and Landau theory, Phys. Rev. B 104, 035142 (2021).
  * [20] R. Tazai, Y. Yamakawa, S. Onari, and H. Kontani, Mechanism of exotic density-wave and beyond-Migdal unconventional superconductivity in kagome metal AV3Sb5 (A = K, Rb, Cs), Sci. Adv. 8, eabl4108 (2022).
  * [21] S. Onari, Y. Yamakawa, and H. Kontani, Phys. Rev. Lett. 112, 187001 (2014).
  * [22] S. Onari, Y. Yamakawa, and H. Kontani, Phys. Rev. Lett. 116, 227001 (2016).
  * [23] Y. Yamakawa, S. Onari and H. Kontani, Phys. Rev. X 6, 021032 (2016).
  * [24] S.-Y. Yang, Y. Wang, B. R. Ortiz, D. Liu, J. Gayles, E. Derunova, R. Gonzalez-Hernandez, L. $\check{S}$mejkal, Y. Chen, S. S. P. Parkin, S. D. Wilson, E. S. Toberer, T. McQueen, and M. N. Ali, Giant, unconventional anomalous Hall effect in the metallic frustrated magnet candidate, KV3Sb5, Sci. Adv. 6, eabb6003 (2020).
  * [25] F. H. Yu, T. Wu, Z. Y. Wang, B. Lei, W. Z. Zhuo, J. J. Ying, and X. H. Chen, Concurrence of anomalous Hall effect and charge density wave in a superconducting topological kagome metal, Phys. Rev. B 104, L041103 (2021).
  * [26] C. Mielke, D. Das, J.-X. Yin, H. Liu, R. Gupta, Y.-X. Jiang, M. Medarde, X. Wu, H. C. Lei, J. Chang, P. Dai, Q. Si, H. Miao, R. Thomale, T. Neupert, Y. Shi, R. Khasanov, M. Z. Hasan, H. Luetkens, and Z. Guguchia, Time-reversal symmetry-breaking charge order in a kagome superconductor, Nature 602, 245 (2022).
  * [27] R. Khasanov, D. Das, R. Gupta, C. Mielke, M. Elender, Q. Yin, Z. Tu, C. Gong, H. Lei, E. T. Ritz, R. M. Fernandes, T. Birol, Z. Guguchia, and H. Luetkens, Time-reversal symmetry broken by charge order in ${\mathrm{CsV}}_{3}{\mathrm{Sb}}_{5}$, Phys. Rev. Research 4, 023244 (2022).
  * [28] Z. Guguchia, C. Mielke, D. Das, R. Gupta, J.-X. Yin, H. Liu, Q. Yin, M. H. Christensen, Z. Tu, C. Gong, N. Shumiya, M. S. Hossain, T. Gamsakhurdashvili, M. Elender, P. Dai, A. Amato, Y. Shi, H. C. Lei, R. M. Fernandes, M. Z. Hasan, H. Luetkens, and R. Khasanov, Tunable unconventional kagome superconductivity in charge ordered RbV3Sb5 and KV3Sb5, Nat. Commun. 14, 153 (2023).
  * [29] C. Guo, G. Wagner, C. Putzke, D. Chen, K. Wang, L. Zhang, M. G. Amigo, I. Errea, M. G. Vergniory, C. Felser, M. H. Fischer, T. Neupert, and P. J. W. Moll, Correlated order at the tipping point in the kagome metal CsV3Sb5, arXiv:2304.00972.
  * [30] F. D. M. Haldane, Model for a Quantum Hall Effect without Landau Levels: Condensed-Matter Realization of the ”Parity Anomaly”, Phys. Rev. Lett. 61, 2015 (1988).
  * [31] I. Affleck and J. B. Marston, Large-n limit of the Heisenberg-Hubbard model: Implications for high-$T_{c}$ superconductors, Phys. Rev. B 37, 3774(R) (1988).
  * [32] C. M. Varma, Theory of the pseudogap state of the cuprates, Phys. Rev. B 73, 155113 (2006).
  * [33] W. H. P. Nielsen, W. A. Atkinson, and B. M. Andersen, Signatures of orbital loop currents in the spatially resolved local density of states, Phys. Rev. B 86, 054510 (2012).
  * [34] R. Tazai, Y. Yamakawa, and H. Kontani, Emergence of charge loop current in the geometrically frustrated Hubbard model: A functional renormalization group study, Phys. Rev. B 103, L161112 (2021).
  * [35] R. Tazai, Y. Yamakawa, and H. Kontani, Charge-loop current order and Z3 nematicity mediated by bond-order fluctuations in kagome metals, Nat. Commun. 14, 7845 (2023).
  * [36] H. Kontani, Y. Yamakawa, R. Tazai, and S. Onari, Odd-parity spin-loop-current order mediated by transverse spin fluctuations in cuprates and related electron systems, Phys. Rev.Research 3, 013127 (2021)
  * [37] M. Tsuchiizu, Yusuke Ohno, Seiichiro Onari, and H. Kontani Orbital Nematic Instability in the Two-Orbital Hubbard Model: Renormalization-Group + Constrained RPA Analysis, Phys. Rev. B 93, 155148.
  * [38] H. Kontani, R. Tazai, Y. Yamakawa, and S. Onari, Unconventional density waves and superconductivities in Fe-based superconductors and other strongly correlated electron systems, Adv. Phys. 70, 355 (2021).
  * [39] R. Tazai, S. Matsubara, Y. Yamakawa, S. Onari, and H. Kontani, Rigorous formalism for unconventional symmetry breaking in Fermi liquid theory and its application to nematicity in FeSe, Phys. Rev. B 107, 035137 (2023).
  * [40] F. Grandi, A. Consiglio, M. A. Sentef, R. Thomale, and D. M. Kennes, Theory of nematic charge orders in kagome metals, Phys. Rev. B 107, 155131 (2023)
  * [41] X. Wu, T. Schwemmer, T. Müller, A. Consiglio, G. Sangiovanni, D. Di Sante, Y. Iqbal, W. Hanke, A. P. Schnyder, M. M. Denner, M. H. Fischer, T. Neupert, and R. Thomale, Nature of Unconventional Pairing in the Kagome Superconductors $A{{V}}_{3}{{Sb}}_{5}$ ($A={K},{Rb},{Cs}$), Phys. Rev. Lett. 127, 177001 (2021).
  * [42] J. Huang, Y. Yamakawa, R. Tazai, and H. Kontani, Odd-parity intra-unit-cell bond-order and induced nematicity in kagome metal CsTi3Bi5 driven by quantum interference mechanism, arXiv:2305.18093 (available at https://arxiv.org/abs/2305.18093).
  * [43] A. V. Chubukov, M. Khodas, and R. M. Fernandes, Phys. Rev. X 6, 041045 (2016).
  * [44] M. Tsuchiizu, K. Kawaguchi, Y. Yamakawa, and H. Kontani, Phys. Rev. B 97, 165131 (2018).
  * [45] C. Guo, C. Putzke, S. Konyzheva, X. Huang, M. Gutierrez-Amigo, I. Errea, D. Chen, M. G. Vergniory, C. Felser, M. H. Fischer, T. Neupert, and P. J. W. Moll, Switchable chiral transport in charge-ordered Kagome metal CsV3Sb5, Nature 611, 461 (2022).
  * [46] J. Luo, Z. Zhao, Y. Z. Zhou, J. Yang, A. F. Fang, H. T. Yang, H. J. Gao, R. Zhou, and G.-q. Zheng, Possible star-of-David pattern charge density wave with additional modulation in the kagome superconductor CsV3Sb5, npj Quantum Materials 7, 30 (2022).
  * [47] M. H. Christensen, T. Biro, B. M. Andersen, and R. M. Fernandes, Loop currents in AV3Sb5 kagome metals: Multipolar and toroidal magnetic orders, Phys. Rev. B 106, 144504 (2022)


∎

11institutetext: Saha Institute of Nuclear Physics, HBNI, 1/AF Bidhannagar,
Kolkata 700064, India.
11email<EMAIL_ADDRESS>
# Recent progress on cluster and meron algorithms for strongly correlated
systems

Debasish Banerjee

(Received: date / Accepted: date)

###### Abstract

Ab-initio studies of strongly interacting bosonic and fermionic systems is
greatly facilitated by efficient Monte Carlo algorithms. This article
emphasizes this requirement, and outlines the ideas behind the construction of
the cluster algorithms and illustrates them via specific examples. Numerical
studies of fermionic systems at finite densities and at real-times are
sometimes hindered by the infamous sign problem, which is also discussed. The
construction of meron cluster algorithms, which can solve certain sign
problems are discussed. Cluster algorithms which can simulate certain pure
Abelian gauge theories (realized as quantum link models) are also discussed.

###### Keywords:

Cluster Algorithms Fermionic and Spin models Lattice Gauge Theories

##  1 Introduction: why cluster algorithms?

The ab-initio studies of strongly correlated systems occurring in Nature,
whether in particle physics or in condensed matter physics, is an extremely
challenging topic. Analytical solutions are difficult to find in most cases,
and perturbation theory fails for strong couplings. Certain weak coupling
methods (such as the epsilon-expansion Card2008 ; Peli2002 ) have been
successful in addressing the existence of fixed points in renormalization
group flow, as well as compute critical exponents at phase transitions by
systematically improving over the mean field estimates. Large-N methods have
provided another analytic handle on some interesting quantum field theories
(QFTs) Mosh2002 . In the case of conformal field theories (CFTs) there has
been exciting developments through the use of AdS/CFT Mald1998 , and the more
recent conformal bootstrap Rata2008  and the large charge expansion Hell2015 .
Tensor Networks Scho2011  have significantly contributed to computing
paradigms in lower dimensional systems, by opening up the possibility to
simulate a large range of strongly interacting systems, even in real-time.

However, in the overwhelmingly large majority of cases, the Markov Chain Monte
Carlo (MCMC) methods, starting from the initial proposal of Metropolis et al.
Metr1953  have provided an unbiased ab-initio method to numerically sample
multidimensional integrals and evaluate the expectation values of physical
operators. Let us denote the degrees of freedom in a (classical or quantum)
system by $\\{q_{i}\\},i=1,\cdots,V$, and the partition function for the
system as

${\cal Z}=\int{\cal
D}q~{}\mathrm{e}^{-S(\\{q_{i}\\})}=\mathrm{Tr~{}~{}exp}(-\beta H),$ (1)

where $S(\\{q_{i}\\})$ is an action functional corresponding to the
Hamiltonian $H$ at an inverse temperature $\beta$. The quantity
$\mathrm{e}^{-S(\\{q_{i}\\})}=W(\\{q_{i}\\})$ is the Boltzmann weight and is
typically positive for the chosen computational basis of $\\{q_{i}\\}$. If
this is not the case, we encounter the sign problem and importance sampling
fails. We will assume the positivity of the Boltzmann weight for the moment,
and will come back to the exceptions later. The expectation value of a
physical operator $\mathcal{O}$ (for example an order parameter or a
correlation function) is

$\langle\mathcal{O}\rangle=\frac{1}{{\cal Z}}\int{\cal
D}q~{}\mathcal{O}~{}\mathrm{e}^{-S(\\{q_{i}\\})}$ (2)

Starting from an initial probability distribution, a Markov Chain is a series
of probability distributions $\Pi_{k}(\\{q_{i}\\})$, which steadily converge
to the fixed point distribution $\Pi^{\star}(\\{q_{i}\\})$, which is the
equilibrium distribution $W(\\{q_{i}\\})$. We will assume that the reader is
familiar with the basics of Markov Chains, and local Monte Carlo algorithms
and point to excellent textbooks which cover this topic in depth MC1 ; MC2 .
While it is not always difficult to construct a Markov Chain with the
necessary properties, the difficulty lies ensuring that the rate of
convergence is insensitive to V. Further, even if a Markov chain were to
converge to the equilibrium distribution, one needs to sample this
distribution to obtain uncorrelated samples on which expectation values and
correlation functions of relevant operators can be measured. This is exactly
where improved algorithms, such as the cluster algorithm, assert their
importance. To make this notion quantitative, we first introduce the concept
of autocorrelation time.

Operationally, after equilibrium is reached, (local) Monte Carlo algorithm
generates a set of configurations of $(\\{q_{i}\\})$ according to the
Boltzmann weight, on which physical operators are measured. Let us denote the
distribution at the $l$-th Monte-Carlo time as $(\\{q_{i}\\})_{l}$. Typically
the distributions $(\\{q_{i}\\})_{l}$ and $(\\{q_{i}\\})_{l+1}$ are highly
correlated since only a small change is involved in one step to the next.
Therefore, measurement of the physical operator on these two subsequent
configurations, $\mathcal{O}_{l}(\\{q_{i}\\})$ and
$\mathcal{O}_{l+1}(\\{q_{i}\\})$ are also not independent.

Numerical calculations necessarily deal with finite data, let us denote this
by $N_{\rm conf}$. The sample mean of the dataset is
$\overline{\mathcal{O}}=\frac{1}{N_{\rm conf}}\sum_{l}\mathcal{O}_{l}$.
According to the central limit theorem, the sample mean
$\overline{\mathcal{O}}$ has a Gaussian distribution about the exact
expectation value $\langle\mathcal{O}\rangle$, such that
$\langle\mathcal{O}\rangle=\overline{\mathcal{O}}\pm\sigma_{\overline{\mathcal{O}}}$.
With uncorrelated measurements, the best unbiased estimate of the error from
the finite data-set is

$\sigma_{\overline{\mathcal{O}}}=\left[\frac{1}{N_{\rm conf}(N_{\rm
conf}-1)}\sum_{k}(\mathcal{O}_{k}-\overline{\mathcal{O}})^{2}\right]^{\frac{1}{2}}.$
(3)

Note that we have used $\overline{\mathcal{O}}$ as the best estimate of the
true mean, $\langle\mathcal{O}\rangle$. For large $N_{\rm conf}$,

$\sigma^{2}_{\overline{\mathcal{O}}}=\left[\frac{1}{N_{\rm
conf}}\sum_{m}(\mathcal{O}_{m}-\overline{\mathcal{O}})\frac{1}{N_{\rm
conf}}\sum_{n}(\mathcal{O}_{n}-\overline{\mathcal{O}})\right]=\frac{1}{N^{2}_{\rm
conf}}\sum_{m,n}^{N_{\rm conf}}\Gamma_{\rm\mathcal{O}}(m-n).$ (4)

The autocorrelation function, $\Gamma_{\rm\mathcal{O}}(m)$, is a property of
the Monte Carlo algorithm and is defined through the unequal time-correlator

$\Gamma_{\rm\mathcal{O}}(m)=\langle\mathcal{O}_{m^{\prime}}\mathcal{O}_{m+m^{\prime}}\rangle-\langle\mathcal{O}\rangle^{2}=\Gamma_{\rm\mathcal{O}}(-m)\overset{m\gg
1}{\sim}C\exp(-m/\tau_{\rm exp})$ (5)

and decays exponentially at asymptotic times, with the decay time $\tau_{\rm
exp}$ corresponding to the slowest mode in the Monte Carlo dynamics. The
integrated autocorrelation time, $\tau_{\rm int}$, defined as

$2\tau_{\rm
int}\approx\sum_{m=-M}^{M}\Gamma_{\rm\mathcal{O}}(m)=\Gamma_{\rm\mathcal{O}}(0)+2\sum_{m=1}^{M}\Gamma_{\rm\mathcal{O}}(m)$
(6)

and the cut $M$ is due to finiteness of the dataset, and we have ignored
corrections of order $\tau_{\rm exp}/N_{\rm conf}$, for large $N_{\rm conf}$.
Note that $\Gamma_{\rm\mathcal{O}}(0)$ is the sample variance. The expression
for the error on the mean is then

$\sigma^{2}_{\rm\mathcal{O}}=\frac{2\tau_{\rm
int}\Gamma_{\rm\mathcal{O}}(0)}{N_{\rm
conf}}=\frac{\Gamma_{\rm\mathcal{O}}(0)}{N_{\rm eff}}$ (7)

implying that there are only $N_{\rm eff}$ independent configurations in the
Markov chain. Unlike the equilibrium values, $\tau_{\rm int}$ depends on the
algorithm. Near a critical point, for example, autocorrelation times diverge
with the physical correlation length $\xi$ as $\tau\sim c\xi^{z}$, and $z\sim
2$ for local algorithms ($z$ is called the dynamical critical exponent). Such
studies are therefore numerically challenging especially near the critical
point, and the phenomena is called critical slowing down. Cluster algorithms
can guarantee extremely small $\tau_{\rm int}$, and sometimes even $z\sim 0$,
thereby practically eliminating this problem.

##  2 Quantum spin models

Using the formulation of Fortuin and Kasteyln Kast1969 ; Fort1972 , the first
cluster algorithm for simulating classical spin (the Ising and the Potts)
models were constructed by Swendsen and Wang Swen1987 , and then extended by
Niedermayer Nied1988  and Wolff Wolf1989  for continuous spins systems. We
refer to the review Nied1997  for more details. Cluster algorithms for
continuous systems, such as hard spheres, have also been developed. We do not
discuss them here, but refer to Krau2003  for a more detailed exposition. We
instead focus on cluster algorithms for quantum systems, starting with quantum
spins.

Cluster algorithms for quantum spin systems, in the world line formulation of
quantum spins, were first developed in Ever1993 ; Ever21993 , and extended to
the continuous time version in Bear1996 . We note that the Stochastic Series
Expansion (SSE) of the quantum spin Hamiltonian as developed in Sand1991
leads to a loop Monte Carlo updating method Sand2002 , similar to the cluster
algorithm. There has been significant developments in generalizing this
algorithm for a variety of models, larger quantum spins, and we refer the
reader to Ever2002  for details and a more complete set of references. We note
that the worm algorithm is another powerful idea which has led to the
development of efficient algorithm for a large class of models, and in many
cases is as powerful as the cluster algorithms. We refer to the interested
reader to Prok2001 ; Boni2006  for details.

Let us illustrate the construction of a cluster algorithm for the Heisenberg
model, which not only has a very rich physical pedagogy behind it, but is also
useful to understand the magnetic properties in certain electronic systems
Heis1928 . This method is independent of spatial dimensions, but we will
consider $(2+1)-$d, anti-ferromagnetic version $J>0$ for definiteness.
Starting from the Hamiltonian on a square lattice of linear extent $L$, the
partition function ${\cal Z}$ at an inverse temperature $\beta$ is:

${\cal Z}=\mathrm{Tr~{}~{}exp}(-\beta
H);~{}~{}~{}H=J\sum_{x,\hat{i}=1,2}\vec{S}_{x}\cdot\vec{S}_{x+\hat{i}}.$ (8)

Using the Suzuki-Trotter formula Suzu1976 , we separate the Hamiltonian into 4
parts ($2d$ parts in $d-$spatial dimensions), $H=H_{1}+H_{2}+H_{3}+H_{4}$,
such that the spins contained in each part mutually commute.

$\displaystyle H_{1}$
$\displaystyle=J\sum_{x=(2m,n)}\vec{S}_{x}\cdot\vec{S}_{x+\hat{1}};~{}~{}H_{2}=J\sum_{x=(2m+1,n)}\vec{S}_{x}\cdot\vec{S}_{x+\hat{1}};$
(9) $\displaystyle H_{3}$
$\displaystyle=J\sum_{x=(m,2n)}\vec{S}_{x}\cdot\vec{S}_{x+\hat{2}};~{}~{}H_{4}=J\sum_{x=(m,2n+1)}\vec{S}_{x}\cdot\vec{S}_{x+\hat{2}}.$

We construct ${\cal Z}$ by inserting intermediate time-slices, such that
$\beta=4N\epsilon$, and $\epsilon$ is the temporal lattice spacing. Using
${\cal I}=|n_{k}\rangle\langle n_{k}|$, we expand ${\cal Z}$ as a product over
the matrix elements of the transfer matrix, expressed in the
$S^{z}=\pm\frac{1}{2}$ basis, which is chosen as the computational basis:

$\displaystyle{\cal Z}$ $\displaystyle=\mathrm{Tr~{}~{}exp}(-\beta
H)=\lim_{\begin{subarray}{c}N\to\infty\\\ \epsilon\to
0\end{subarray}}\sum_{n}\langle n|\left[{\rm
e}^{-\epsilon(H_{1}+H_{2}+H_{3}+H_{4})}\right]^{N}|n\rangle$ (10)
$\displaystyle=\lim_{\begin{subarray}{c}N\to\infty\\\ \epsilon\to
0\end{subarray}}\sum_{\begin{subarray}{c}\\{n\\}\\\
n_{4N}=n_{0}\end{subarray}}\prod_{k=0}^{N-1}\langle n_{k}|{\rm e}^{-\epsilon
H_{1}}|n_{k+1}\rangle\langle n_{k+1}|{\rm e}^{-\epsilon
H_{2}}|n_{k+2}\rangle\langle n_{k+2}|{\rm e}^{-\epsilon
H_{3}}|n_{k+3}\rangle\langle n_{k+3}|{\rm e}^{-\epsilon H_{4}}|n_{k+4}\rangle$

We have used the Trotter formula in the second line, and the matrix elements
connect only the specific bonds. Finally, we can express the matrix elements
via the local action

${\cal Z}=\prod_{x,t}\sum_{s(x,t)=\pm 1}{\rm e}^{-S[s_{1},s_{2},s_{3},s_{4}]}$
(11)

where the term $S[s_{1},s_{2},s_{3},s_{4}]$ connects two neighboring spins,
$s_{1},s_{2}$ at a time-slice, $t$, with their forward-in-time partners
$s_{3},s_{4}$. Note that only a set of bonds are active at a given time-slice,
and all others are passive, and this 4-spin interaction traces an active
plaquette. A schematic figure for the Trotter decomposition in $(1+1)-$d is
shown in Fig 1 (left), where the shaded plaquettes are the active ones, and
indicate the spin pairs which are interacting at that given time. A
representative $s_{1},s_{2},s_{3},s_{4}$ are also shown in the same figure.
For the anti-ferromagnetic Heisenberg model, the explicit values of the
transfer matrix (in the basis where the states of the two spins at sites $x$
and $y$ are
$|\uparrow\uparrow\rangle,|\uparrow\downarrow\rangle,|\downarrow\uparrow\rangle,|\downarrow\downarrow\rangle$)
are:

${\rm e}^{-S[s_{1},s_{2},s_{3},s_{4}]}=\langle s_{1}s_{2}|{\rm e}^{-\epsilon
J\vec{S}_{x}\vec{S}_{y}}|s_{3}s_{4}\rangle={\rm e}^{-\frac{\epsilon
J}{4}}\begin{bmatrix}1&0&0&0\\\ 0&\frac{1}{2}(1+{\rm e}^{\epsilon
J})&\frac{1}{2}(1-{\rm e}^{\epsilon J})&0\\\ 0&\frac{1}{2}(1-{\rm e}^{\epsilon
J})&\frac{1}{2}(1+{\rm e}^{\epsilon J})&0\\\ 0&0&0&1\\\ \end{bmatrix}$ (12)

Note that there are only two off-diagonal elements, both of which are
negative. For a bipartite lattice, the negative sign can be eliminated by
doing a unitary transformation on every alternative spin, such that the off-
diagonal elements become $\frac{1}{2}({\rm e}^{\epsilon J}-1)$. For a
triangular lattice, for example (or for any non-bipartite lattice, in
general), this does not work, and we encounter the first example of a sign
problem. In the chosen basis the probability weights are not positive
definite, and Monte Carlo simulation cannot be performed.

Figure 1: (left) An example of a world-line configuration of the anti-
ferromagnetic Heisenberg model. There is one transition and one anti-
transition which spatially displaces a spin-up (down) and brings it back.
(right) Possible bond-types for this model: (a) the A-type breakup connects
identical spins forward in time, (b) the B-type breakup connects opposite
(spin-up with spin-down) spins sideways. We can figure out the probability of
the weights by computing which breakup can be applied for which interaction
plaquette. For those where the spin-states do not, and cannot change (top),
only A can be applied, whereas where the spin state flips (middle) only B can
be applied. The last case is where the spin state can change but did not.
Here, the A-breakup can be applied with the probability $W_{A}=A/(A+B)=2/({\rm
exp}(\epsilon J)+1)$.

Next, to construct a cluster algorithm we have to expand the configuration
space by including bond variables, together with the quantum spins. This
implies choosing the bond variables $[b]$ such that the following equation is
satisfied for the different possibilities of $s_{1},s_{2},s_{3},s_{4}$,

${\rm e}^{-S[s_{1},s_{2},s_{3},s_{4}]}=\sum_{[b=A,B]}{\rm
e}^{-S[s_{1},s_{2},s_{3},s_{4};b]}.$ (13)

In this particular example, as shown in Fig 1, the linear equations can be
satisfied by the use of two bond breakups, $A$ and $B$. The clusters are
constructed as follows: (a) start from an initial site, (b) follow the chosen
breakups as explained in the caption of Fig 1 (right) until the initial point
is reached. In this case, the resulting cluster is the same as a loop, and
hence the relation to the loop algorithm. Flipping the cluster implies the
operation $|\uparrow_{k}\rangle\leftrightarrow|\downarrow_{k}\rangle$, which
preserves the Boltzmann weight of the configuration, and yet changes the
configuration globally, thereby decorrelating them very fast. This is the
reason for the efficiency of the cluster algorithm — the clusters correspond
to correlated degrees of freedom, which get efficiently updated by flips. In
fact, representing the configurations in terms of clusters allows one to
construct improved estimators for quantities such as the magnetization,
susceptibility and the correlation function.

Let us illustrate the construction of improved estimators through examples.
The total (uniform) magnetization of the spins is given as ${\cal
M}[\\{S^{z}\\}]=\sum_{i}(S^{z})_{i}=\sum_{C}{\cal M}_{C}$, where $i$ runs over
all the space-time lattice sites, and $C$ runs over all the clusters into
which a single configuration of spins has been decomposed. The equality
follows since each spin must uniquely belong to a cluster. The cluster
magnetization is thus given as ${\cal M}_{C}=\sum_{i\in C}(S^{z})_{i}$. Since
each cluster flip changes the sign of ${\cal M}_{C}$, we conclude that
$\langle{\cal M}\rangle=0$. Note that this is consistent with the physical
result that the magnetization always vanishes in a finite volume. Such results
are extremely difficult to observe, especially in the broken phase, with local
update algorithms. Another example is the (uniform) susceptibility:

$\chi=\frac{\beta}{V}\langle{\cal
M}^{2}\rangle=\frac{\beta}{V}\langle(\sum_{C}{\cal
M}_{C})^{2}\rangle=\frac{\beta}{V}\langle\sum_{C_{i},C_{j}}{\cal
M}_{C_{i}}{\cal M}_{C_{j}}\rangle=\frac{\beta}{V}\langle\sum_{C}{\cal
M}^{2}_{C}\rangle$ (14)

Note that different clusters $C_{i}$ and $C_{j}$ are uncorrelated and hence
the cross-terms vanish under the operation of all cluster flips, leaving
behind only diagonal terms. Further in a cluster $C$, since all the spins are
pointing in the same direction, the magnetization is given by the total
cluster size, ${\cal M}_{C}=\pm|C|$. This allows us to rewrite $\chi$ as
follows, $\chi=\frac{\beta}{V}\langle\sum_{C}|C|^{2}\rangle$. Thus it is
evident that clusters are relevant physical objects, since their size is
related to a physical quantity.

Two variants of this cluster algorithm are well-known: multi-cluster
algorithm, which proceeds to construct all the possible clusters on a given
configuration of spins and bonds, and then flips each cluster with a $50\%$
probability, and the single-cluster, or the Wolff algorithm which constructs a
single cluster at always flips it. It is also possible to construct this
cluster algorithm in continuous time.

It is to be noted that both types of cluster algorithms suffer from Trotter
errors, caused due to the finite $\epsilon$ essential in the numerical
simulations. In particular, the errors on the physical quantities go as
$\epsilon^{2}$, and thus one typically simulates for several different
$\epsilon$ values to take the time-continuum limit. It is however possible to
construct the cluster algorithms directly in the time-continuum limit Bear1996
, which eliminate the Trotter errors completely.

It must be emphasized that for the cluster algorithm to be successful, certain
special configurations called reference configurations need to exist. In the
case of the Heisenberg anti-ferromagnet the reference configuration(s) are the
ones where the spins are staggered on the bi-partite lattice. Any
configuration can then be decomposed into a certain set of clusters, and the
clusters are accordingly flipped to bring the configuration into (one of) the
reference configurations. This trivially proves ergodicity. In many cases
(examples are frustrated magnets, spin glasses, gauge theories), while it is
conceivable to build clusters, absence of reference configurations causes the
clusters to grow too big, almost filling up the whole volume. In such a case,
the cluster algorithm does not perform any better than a local update
algorithm.

##  3 Fermions

Attempting to simulate fermions brings us head on with the fermion sign
problem. This is most directly seen in the world-line representation of the
fermions in the occupation number basis. In dimensions $d\geq 2$, it is easy
to construct world-lines of identical fermions which twist around each other
as they evolve in Euclidean time. Configurations which have fermions
exchanging positions an odd number of times have an overall negative sign
compared with those where fermions exchange positions an even number of times
(everything else being identical), due to the Pauli exclusion principle. A
particularly clear statement of this fermion sign problem and what it entails
to solve the problem is given in Troy2005 .

As in the case of quantum spins, to simulate fermionic systems we construct
the partition function as usual. In the Lagrangian formulation, one uses
Grassmann variables, which is typically integrated out to leave behind a
determinant involving auxiliary fields (if the theory has four-Fermi
interactions, or Yukawa couplings), or even gauge fields (as in quantum
electrodynamics, or quantum chromodynamics). Updating procedures either
involve updating the determinant directly Blan1981 ; Assa2008 , or recast the
determinant in terms of bosonic fields, which are then updated Duan1987 . To
construct cluster algorithms for fermions, we first obtain a bosonic
representation using the Jordan-Wigner transformation. If the resulting
bosonic system can be efficiently updated using a cluster algorithm, then one
can identify types of interactions for which the fermion sign problem can be
solved Chan2002 . When such an approach succeeds, it is known as the meron
algorithm first introduced in Chan1999 , and will be discussed below. Other
novel approaches, such as the fermion bag Chan2010 ; Huff2017  or the
diagrammatic Monte Carlo Boni2006 ; Houc2010  are getting increasingly popular
in simulating fermionic systems. The former method expands the fermionic
action and groups regions where the fermions can propagate freely against
regions where they are bound into monomers and dimers. The Monte Carlo
procedure then updates these regions, which are denoted as fermion bags. The
latter approach directly samples the Feynman diagrams associated with the
interactions, and in certain regimes can be related to the bag approach.

Before discussing the construction of the meron algorithm, we elaborate the
nature of the fermionic sign problem, which will also serve to understand the
solution. First, we set up the fermionic path integral in the occupation
number basis $\\{n_{k}\\}$:

${\cal Z}_{\rm f}=\mathrm{Tr~{}~{}exp}(-\beta H_{\rm
f})=\sum_{\\{n_{k}\\}}p(\\{n_{k}\\});\langle A\rangle=\frac{1}{{\cal Z}_{\rm
f}}\sum_{\\{n_{k}\\}}A(\\{n_{k}\\})p(\\{n_{k}\\}),$ (15)

where $p(\\{n_{k}\\})$ is the probability for a given configuration, and $A$
is an operator in the occupation number basis. To qualify the severity of this
problem, one can consider the sign of the configuration as a part of the
observable $A$ through the decomposition $p(\\{n_{k}\\})={\rm
sign}(\\{n_{k}\\})|p(\\{n_{k}\\})|$, the sign being $\pm 1$ depending on
whether the fermions in the configuration swap their positions even or odd
number of times. It follows that

$\displaystyle\langle A\rangle_{\rm f}$
$\displaystyle=\frac{\sum_{\\{n_{k}\\}}A(\\{n_{k}\\})p(\\{n_{k}\\})}{\sum_{\\{n_{k}\\}}p(\\{n_{k}\\})}$
(16) $\displaystyle=\frac{\sum_{\\{n_{k}\\}}A(\\{n_{k}\\}){\rm
sign}(\\{n_{k}\\})|p(\\{n_{k}\\})|}{\sum_{\\{n_{k}\\}}|p(\\{n_{k}\\})|}\times\frac{\sum_{\\{n_{k}\\}}|p(\\{n_{k}\\})|}{\sum_{\\{n_{k}\\}}{\rm
sign}(\\{n_{k}\\})|p(\\{n_{k}\\})|}$ $\displaystyle=\frac{\langle A\cdot{\rm
sign}\rangle_{\rm b}}{\langle{\rm sign}\rangle_{\rm b}}$

The subscript ${\bf b}$ indicates that the averaging is done over a bosonic
system whose weights, $|p(\\{n_{k}\\}|$, are positive definite by
construction. The quantity $\langle{\rm sign}\rangle_{\rm b}$ measures the
severity of the sign problem. It can shown that $\langle{\rm sign}\rangle_{\rm
b}={\rm exp}(-\beta V\Delta f)$, where $\Delta f=f_{\rm f}-f_{\rm b}$ is the
difference in free energy density between the fermion and bosonic ensembles.
At low temperatures and large volumes, this quantity is exponentially small:

$\frac{(\sigma_{\rm sign})_{\rm b}}{\langle\rm
sign\rangle_{b}}=\frac{\left[\langle{\rm sign}^{2}\rangle_{\rm b}-\langle{\rm
sign}\rangle^{2}_{\rm b}\right]^{\frac{1}{2}}}{\sqrt{N}\langle\rm
sign\rangle_{b}}\approx\frac{{\rm exp}(\beta V\Delta f)}{\sqrt{N}},$ (17)

where $N$ is the number of uncorrelated data, and we have used $\langle{\rm
sign}^{2}\rangle_{\rm b}=1,\langle{\rm sign}\rangle_{\rm b}\approx 0$.
Clearly, $N\sim{\rm exp}(2\beta V\Delta f)$ to determine the sign. Thus even
if $\langle A\rangle_{\rm f}\sim 1$, it is measured as the ratio of two
exponentially small signals, and requires exponential effort to extract from
statistical noise.

One could aim to cancel the negative signs by pairing the configurations which
carry an $-1$ sign with those that carry a $1$ sign, such that only a few
unmatched configurations remain. In fact, this is the idea behind the meron
algorithm. It is however important to realize that this solves only one-half
of the sign problem. With the matching, we effectively have ${\rm
sign}^{2}={\rm sign}$, since ${\rm sign}=0,1$ and hence

$\frac{(\sigma_{\rm sign})_{\rm b}}{\langle\rm
sign\rangle_{b}}=\frac{\left[\langle{\rm sign}\rangle_{\rm b}-\langle{\rm
sign}\rangle^{2}_{\rm b}\right]^{\frac{1}{2}}}{\sqrt{N}\langle\rm
sign\rangle_{\rm b}}\approx\frac{{\rm exp}(\beta V\Delta f/2)}{\sqrt{N}},$
(18)

Note that we have achieved an exponential gain in statistics, but one still
needs $N\sim{\rm exp}(\beta V\Delta f)$ for any meaningful result. The meron
algorithm incorporates the two steps together: it identifies the
configurations which can be analytically cancelled, and never generates them.
Thus, one always generates configurations which contribute non-trivially to
${\cal Z}_{\rm f}$. The algorithm, however, only works for a restricted class
of interactions.

Let us illustrate this for the ${\rm t}-{\rm V}$ model in $(2+1)-$d,
extensively used in the study of certain electronic systems (though the
algorithm can be constructed in any dimensions). The Hamiltonian is

${\rm H}_{\rm f}=-\frac{{\rm
t}}{2}\sum_{x,\hat{i}=1,2}(c^{\dagger}_{x}c_{x+\hat{i}}+c^{\dagger}_{x+\hat{i}}c_{x})+{\rm
V}\sum_{x,\hat{i}=1,2}(n_{x}-\frac{1}{2})(n_{x+\hat{i}}-\frac{1}{2}),$ (19)

with $V\geq t>0$. The fermion creation (annhiliation) operators at site x,
$c_{x}(c^{\dagger}_{x})$ satisfy the following anti-commutation relations:
$\\{c_{x},c_{y}\\}=\\{c^{\dagger}_{x},c^{\dagger}_{y}\\}=0;~{}~{}\\{c^{\dagger}_{x},c_{y}\\}=\delta_{xy}$.
The corresponding bosonic system is obtained by using the Jordan-Wigner
transformation. An ordering of the lattice points is defined: $l=x+y\cdot L$
(in $d=2$), and the fermionic operators are expressed as follows:

$c^{\dagger}_{x}=\sigma^{3}_{1}\cdot\sigma^{3}_{2}\cdots\sigma^{3}_{l-1}\sigma^{+}_{l};~{}~{}c_{x}=\sigma^{3}_{1}\cdot\sigma^{3}_{2}\cdots\sigma^{3}_{l-1}\sigma^{-}_{l};~{}~{}n_{x}=c^{\dagger}_{x}c_{x}=\frac{1}{2}(\sigma^{3}_{l}+1).$
(20)

This preserves the anti-commutation relations. We use the (fermion) occupation
number basis ($n_{x}|0\rangle=0,n_{x}|1\rangle=|1\rangle$, with $0$ denoting
an unoccupied and $1$ an occupied site) to construct ${\cal Z}_{\rm f}$ by
splitting the Hamiltonian into different parts and invoking the Suzuki-Trotter
formula. With some algebra to keep track of the negative signs, the transfer
matrix between the adjacent time-slices can be expressed through the following
$4\times 4$ plaquette interactions (the state labels at sites $x$ and $y$
denote $|00\rangle,|01\rangle,|10\rangle,|11\rangle$):

$\displaystyle{\rm e}^{-S[n_{1},n_{2},n_{3},n_{4}]}$ $\displaystyle=\langle
n_{1}n_{2}|{\rm e}^{\epsilon\frac{{\rm
t}}{2}(c^{\dagger}_{x}c_{y}+c^{\dagger}_{y}c_{x})-\epsilon
V(n_{x}-\frac{1}{2})(n_{y}-\frac{1}{2})}|n_{3}n_{4}\rangle$ (21)
$\displaystyle={\rm e}^{\frac{\epsilon V}{4}}\begin{bmatrix}{\rm e}^{-\epsilon
V/2}&0&0&0\\\ 0&{\rm cosh}(\epsilon{\rm t}/2)&\Sigma{\rm sinh}(\epsilon{\rm
t}/2)&0\\\ 0&\Sigma{\rm sinh}(\epsilon{\rm t}/2)&{\rm cosh}(\epsilon{\rm
t}/2)&0\\\ 0&0&0&{\rm e}^{-\epsilon V/2}\\\ \end{bmatrix}$

The structure of the transfer matrix is very similar to the one of the quantum
spins, the main difference is the string operator
$\Sigma=\sigma^{3}_{l+1}\cdot\sigma^{3}_{l+2}\cdots\sigma^{3}_{m-1}$ between
the sites $x=l$ and $x+\hat{i}=m$, where $l<m$ in the lattice ordering. This
allows the rewriting of ${\cal Z}_{\rm f}$ in the occupation number basis
$[\\{n_{x}=0,1\\}]$ as

${\cal Z}_{\rm f}=\sum_{\\{n\\}}{\rm sign}[\\{n\\}]{\rm e}^{-S[\\{n\\}]}.$
(22)

Figure 2: Example configurations which contributes to ${\cal Z}_{\rm f}$.
(left) The two fermions $f_{1}$ and $f_{2}$ are static during their entire
Euclidean time-evolution, and trace out vertical worldlines. For $t=J$, only
$A$ (vertical) and $B$ (horizontal) bonds are required. The figure shows one
bond configuration. The oriented clusters are constructed by starting from a
single point and following the bonds – upwards (downwards) for a filled
(empty) site for $A-$breakups, sideways for a filled (empty) site touching
$B-$breakups. This choice of breakups gives rise to 5 different clusters, as
marked in the figure. (middle) On flipping clusters marked 1, 3 and 4 we
obtain a completely different world-line configuration, but one where the
fermions do not exchange positions. (right) Flipping cluster 5 gives rise to a
configuration where the fermions interchange positions, and hence due to the
anticommutation relations, this configuration has an overall negative sign
compared to the one on the left. The cluster 5 is a meron.

Some example configurations which contribute to ${\cal Z}_{\rm f}$ are shown
in Fig. 2. Without the sign factor, the weights of the resulting bosonic
system are identical to that of the XXZ-Hamiltonian:

$H=\sum_{x,\hat{i}=1,2}\left[{\rm t}(S^{1}_{x}\cdot
S^{1}_{x+\hat{i}}+S^{2}_{x}\cdot S^{2}_{x+\hat{i}})+{\rm V}S^{3}_{x}\cdot
S^{3}_{x+\hat{i}}\right],$ (23)

which reduces to the anti-ferromagnet for ${\rm t}={\rm V}=J$. The ${\cal
Z}_{\rm f}$ is decomposed in terms of bonds, in addition to spins:

${\cal Z}_{\rm f}=\sum_{[n,b]}{\rm sign}[n,b]~{}{\rm
exp}(\left[-S[n,b]\right])$ (24)

Interestingly, in the limit ${\rm t}={\rm V}$ the only breakups that are
needed are the $A$ and the $B$ breakup, with same probabilities as derived in
Fig. 1 (right). Fig 2 shows a typical fermion configuration, and a particular
set of breakups, and how clusters are identified. Cluster flips are operations
$\\{n_{k}\leftrightarrow 1-n_{k}|n_{k}\in{\cal C}_{i}\\}$, where ${\cal
C}_{i}$ is the i-th cluster. Cluster flips give rise to new worldlines
significantly different from their parent ones, but carry the same weights.

The crucial consequence of these cluster rules is that the sign of the whole
configuration factorizes into a product of the signs associated with each
individual cluster. An example is already seen in Fig 2 (right). The clusters
which can change the sign of the configuration are called merons. It is
possible to reach a reference configuration by flipping clusters
appropriately. In this example, the configuration in Fig 2 (left) is the
reference configuration, which has a positive sign (${\rm sign}_{{\cal
C}_{i}}=1$) for all clusters i. When a cluster is flipped, its contribution to
sign is ${\rm sign}_{{\cal C}_{i}}=1(-1)$ if the quantity $N_{w}+N_{h}/2$ is
odd (even). $N_{w}$ is the temporal winding number of the cluster, while
$N_{h}$ is the number of spatial hops Chan1999 . The meron concept allows
exact pairing of even and odd signs:

$\displaystyle{\cal Z}_{\rm f}$ $\displaystyle=\sum_{[n,b]}\overline{{\rm
sign}[b]}~{}{\rm exp}(-S[n,b])=\sum_{[b],{\rm zero-meron}}2^{N_{\cal
C}}~{}{\rm exp}(-S[b]);$ (25) $\displaystyle\overline{{\rm sign}[b]}$
$\displaystyle=\frac{1}{2^{N_{\cal C}}}\sum_{\rm clusterflips}{\rm
sign}[n,b].$

$\overline{{\rm sign}[b]}=0$ if at least one of the clusters is a meron.
Further, the existence of the reference configuration with a positive
Boltzmann weight, guarantees that the unpaired configurations all carry
positive weight, and can be reached by flipping clusters from an initial
configuration which has zero-merons. The Monte Carlo sampling proceeds by
generating configurations which only lie in the zero-meron sector. If a
proposed flip gives rise to meron cluster, it is rejected. Thus, one is able
to completely solve the sign problem with the meron concept.

Let us emphasize that for generic interactions the meron concept does not
hold, and whether a cluster is a meron depends on the orientation of other
clusters. Therefore, much more (exponential) computational effort is needed to
identify the merons. However, there is a large class of interactions which can
be solved with the meron concept, and we refer the reader to Chan2002  for a
review.

##  4 Pure Gauge Theories

Cluster algorithms for gauge theories run into problems since the relevant
variables that need to be updated get frustrated, and very often the clusters
occupy the entire volume. The resulting algorithm is not much better than a
local update algorithm. One direction of research tried to identify the
relevant degrees of freedom, which may or may not be identical to the degrees
of freedom in which the action is constructed. For example, for a $\phi^{4}$
theory, the $\phi$ field can be separated into an Ising ($Z_{2}$) variable and
a modulus variable. The $Z_{2}$ degrees of freedom were updated using a
cluster algorithm while the modulus of the $\phi$ field was updated with a
local update algorithm Brow1989 , which resulted in an efficient algorithm
with a low ($z<1$) dynamical exponent. Similar methods have been applied to
the $SU(2)$ lattice gauge theory Ever1991 ; Kerl1994  targeting the Polyakov
loop as the embedded degree of freedom. While some success on the theories
with one and two time-slices was reported, these methods do not work well for
larger lattices.

Figure 3: (left) Layout of the plaquette, and the four spin interaction.
(right) The action of the operators in the Hamiltonian on the possible
plaquette states. The $U_{\Box}(U^{\dagger}_{\Box})$ term flips a clockwise
(anti-clockwise) oriented plaquette to an anti-clockwise (clockwise) oriented
one, and annhiliates the remaining 14 possible states. This is the $J-$ term
of the Hamiltonian, and is an off-diagonal operator. The $\lambda-$term counts
the total number of flippable (both in the clockwise and anti-clockwise
orientation) plaquettes.

A completely different formulation of gauge theories, the quantum link model
formulation Horn1981 ; Orla1990 ; Chan1997 , which realizes continuous gauge
symmetries with finite dimensional Hilbert spaces and generalizes the Wilson
construction of lattice gauge theories, has recently been more amenable to
simulation with cluster algorithms. As an illustrative example, we will
describe the $U(1)$ quantum link model (QLM). The degrees of freedom of this
Abelian lattice gauge theory are quantum links, $U_{x,\hat{i}}$ and electric
fluxes, $E_{x,\hat{i}}$, both of which are operators defined on the bonds
connecting neighboring sites, $x$ and $x+\hat{i}$ and $\hat{i}=1,2$ in $d=2$.
These operators satisfy canonical commutation relations:

$[E_{x,\hat{i}},U_{y,\hat{j}}]=U_{x,\hat{i}}\delta_{xy}\delta_{ij};~{}[E_{x,\hat{i}},U^{\dagger}_{y,\hat{j}}]=-U^{\dagger}_{x,\hat{i}}\delta_{xy}\delta_{ij};~{}[U_{x,\hat{i}},U^{\dagger}_{y,\hat{j}}]=2E_{x,\hat{i}}\delta_{xy}\delta_{ij},$
(26)

The operators $E,U$ can be chosen to be the components of a spin-$S$ object:
$U=S^{+},U^{\dagger}=S^{-},E=S^{3}$. The parameter $S$ can be thought of
regulating the local Hilbert space dimension. In the limit $S\to\infty$
Schl2000 , these models reproduce the Hamiltonian formulation of Wilson gauge
theories Kogu1974 . The Hamiltonian we will be interested in only contains
operators which commute with the Gauss’ Law:

$H=-J\sum_{\Box}\left(U_{\Box}+U^{\dagger}_{\Box}\right)+\lambda\sum_{\Box}\left(U_{\Box}+U^{\dagger}_{\Box}\right)^{2};G_{x}=\sum_{i}\left(E_{x,x+i}-E_{x-i,x}\right)$
(27)

In addition, we could have included any function of the electric fluxes, and
in particular, the kinetic energy of the gauge fields $\sum_{x,i}E^{2}_{x,i}$.
When the quantum link operators are considered in the spin-$\frac{1}{2}$
representation, this model has several interesting physical applications. In
$(2+1)-$d, the physics of this model is relevant for the understanding the
behavior of low-temperature frustrated magnets Shan2004 , and spin-liquid
phases Herm2004 ; Bale2010 . A close cousin of this model is the quantum dimer
model (QDM), especially well-known in condensed matter physics as a toy model
to describe certain aspects of superconductivity Rokh1988 ; Sach1989 ;
Moes2007 . Therefore, we continue with the gauge links in the
spin-$\frac{1}{2}$ representation, and hence ignore the kinetic energy term of
the gauge links, which is a constant. The local link Hilbert space is
2-dimensional, and in the electric flux basis they are denoted by upward
(downward) arrows for $E=\frac{1}{2}(-\frac{1}{2})$ for vertical links, and
right (left) arrow for $E=\frac{1}{2}(-\frac{1}{2})$ for horizontal links. The
action of the Hamiltonian on the electric flux states is shown in Fig 3.

Specifying the Gauss’ Law further specifies which basis states can contribute
to the partition function, and which not. For the QLM, the Gauss’ Law is
chosen as $G_{x}|\psi\rangle=0$, for every site $x$, for every eigenstate
$|\psi\rangle$ of the Hamiltonian. The QDM chooses a different set of states,
which are specified by $G_{x}|\psi\rangle=(-1)^{x_{1}+x_{2}}|\psi\rangle$.
Physically, this implies that the vacuum selected by the QLM is charge neutral
at each vertex, while the QDM chooses a vacuum with staggered positive and
negative unit charges. The link states allowed at the vertex for the QLM is
shown in Fig 4 (left) and for the QDM in Fig 4 (right). The Gauss’ Law
generates the gauge transformations, which can be generically denoted as
$V=\prod_{x}{\rm exp}(i\theta_{x}G_{x})$, where $\theta_{x}\in(0,2\pi]$ is the
parameter of the transformation. The Hamiltonian is invariant under this
transformation: $\widetilde{H}=V^{\dagger}\cdot H\cdot V=H$, since
$[G_{x},H]=0$, for all x. Physically, $G_{x}$ labels different super-selection
sectors of the theory which do not mix under unitary dynamics.

Figure 4: (left) The Gauss’ Law condition for the QLM, in $d=2$, this allows 6
possible states on each vertex, while (right) the Gauss’ Law condition for the
QDM in $d=2$ allows for 4 states for each vertex. The red circle represents
$Q=+1$, while the blue circle the $Q=-1$, which are distributed in a staggered
fashion on the square lattice. Figure 5:  Mapping of an electric flux
configuration (shown with arrows on the links) to a height configuration
(shown with $+$ and $-$ variables at the centers of the plaquettes). The $+$
and $-$ are placeholders for $0,1(\pm\frac{1}{2})$ on A (B) sub-lattices.
Every time a flux pointing right or upwards (corresponding to $E=\frac{1}{2}$)
is crossed, the height variable is changed, while it remains unchanged if a
left or downward pointing flux (corresponding to $E=-\frac{1}{2}$) is crossed.
The configuration shown has all plaquettes flippable, and the corresponding
height variables are in their reference configuration.

The partition function of this model is defined by projecting the eigenstates
of the Hamiltonian into a specified Gauss’ Law sector, ${\cal
Z}=\mathrm{Tr}[\mathrm{exp}(-\beta H){\cal P}_{G}]$. The projection operator
${\cal P}_{G}=\prod_{x}G_{x}$, therefore constrains the Hilbert space further.
To construct the cluster algorithm, we note that we can divide the Hamiltonian
into two parts $H=H_{A}+H_{B}$, such that the all terms in each part mutually
commute. Thus, the Trotterization divides the lattice into even and odd sub-
lattices, such that at a given time-step only a single sub-lattice needs to be
updated. For a fixed $\beta=\epsilon N$, we have a two-step transfer matrix

${\cal Z}={\rm Tr}~{}\left[(T_{A}T_{B})^{N}P_{G}\right].$ (28)

The single plaquette transfer matrix can be expressed as

$T_{\Box}=1+(U_{\Box}+U^{\dagger}_{\Box}){\rm e}^{-\epsilon J\lambda}{\rm
sinh}(\epsilon J)+(U_{\Box}+U^{\dagger}_{\Box})^{2}\left[{\rm e}^{-\epsilon
J\lambda}{\rm cosh}(\epsilon J)-1\right]$ (29)

The resulting transfer matrix is 16-dimensional (i.e. a $16\times 16$ matrix),
and it is easy to read off the Boltzmann weights from this equation. For
plaquettes which are not flippable, only the diagonal element is unity, and
all other off-diagonal elements are zero. For the two flippable plaquette
configurations, the diagonal contribution has the weight ${\rm
e}^{-\epsilon\lambda}{\rm cosh}(\epsilon J)$ while the off-diagonal elements
(which indicate the plaquette flips) carry the weight ${\rm
e}^{-\epsilon\lambda}{\rm sinh}(\epsilon J)$.

The idea behind the construction of the cluster algorithm is to dualize the
gauge theory in $(2+1)-$d, which gives rise to a $\mathbb{Z}(2)$ quantum
height model in $(2+1)-$d. The dualization is an exact rewriting of the
partition function in terms of new degrees of freedom which are
$\mathbb{Z}(2)$ degrees of freedom located at dual sites. As shown in Fig 5,
every flux configuration can be mapped to a height configuration. A
configuration of quantum height variables is assigned the values
$h^{A}_{\tilde{x}}=0,1$ and $h^{B}_{\tilde{x}}=\pm\frac{1}{2}$. located at the
dual sites $\tilde{x}=(x_{1}+\frac{1}{2},x_{2}+\frac{1}{2})$, and is
associated with a flux configuration

$E_{x,\hat{i}}=\left[h^{X}_{\tilde{x}}-h^{X^{\prime}}_{\tilde{x}+\hat{i}-\hat{1}-\hat{2}}\right]{\rm
mod}2=\pm\frac{1}{2};~{}~{}X,X^{\prime}\in\\{A,B\\}.$ (30)

The modulo function acts by adding or subtracting 2 until the result is in the
desired range $(-1,1]$. We note that while for each height configuration there
is exactly one flux configuration, the reverse is not true. It can be shown
that there are exactly two height configurations for each flux configuration.
In addition we note that the mapping to the height variables only works when
(a) there are no charges at the lattice sites and (b) when charges $Q=\pm 2$
are located on the lattice sites. This explicitly indicates the
$\mathbb{Z}(2)$ nature of the height variables.

Figure 6: (left top) Layout of the height variables across the timeslices t-1,
t and t+1. (left bottom) Two kinds of bonds are enough to solve the linear
equations: the A-breakup joins heights, while the B-breakup keeps the heights
separate. Different scenarios need to be considered — out-of-plane (temporal)
bonds, and in-plane (spatial) bonds. (right top) Cluster rules for temporal
bonds: with $m_{1},m_{2},m_{3},m_{4}$ in a reference configuration and
$m_{5}=m_{6}$, both A- and B-breakups can contribute, while if $m_{5}\neq
m_{6}$ then only the B-breakup contributes. To satisfy detailed balance, we
must have $P_{A}=A/(A+B)={\rm e}^{-\epsilon J\lambda}/{\rm cosh}(\epsilon J)$
and $P_{B}={\rm tanh}(\epsilon J)$. In the case when $m_{5}=m_{6}$ and
$m_{1},m_{2},m_{3},m_{4}$ are not in a reference configuration, then we have
to bind $m_{5}$ and $m_{6}$ with probability 1, otherwise a forbidden
configuration will be generated. (right bottom) Spatial breakups: if
$m_{5}\neq m_{6}$ then connect $m_{1},m_{2},m_{3},m_{4}$ with probability 1,
to prevent the generation of a forbidden configuration. If $m_{5}=m_{6}$ then
two different cases can occur: if $m_{1},m_{2},m_{3},m_{4}$ are in a reference
configuration, then the Boltzmann weight is ${\rm e}^{-\epsilon J\lambda}{\rm
cosh}(\epsilon J)$ and we can apply either the A, or the B-breakup. Otherwise,
the weight is 1 and only B breakup is applied ($m_{5}$ and $m_{6}$ are not
bound together). Solving the equations, we get $P_{A}=1-{\rm e}^{-\epsilon
J\lambda}/{\rm cosh}(\epsilon J)$ and $P_{B}=1-P_{A}$.

Before deriving the cluster rules, it is useful to understand the Boltzmann
weights of the configurations in terms of the height variables, as shown Fig
6. In terms of the height variables, the transfer matrix gives rise to a
6-height interaction between the heights $m_{1},\cdots,m_{6}$. Of these,
$m_{1},m_{2},m_{3},m_{4}$ live on a timeslice t, $m_{5}$ lives on the
timeslice t-1 and $m_{6}$ on t+1. Depending on the values of
$m_{1},m_{2},m_{3},m_{4}$, the height $m_{5}$ can undergo a transition to
$m_{6}$ (corresponding to a plaquette flip), or not. In this language if
$(m_{1},m_{2},m_{3},m_{4})=(1,1,0,0)$ or $(0,0,1,1)$ the plaquette is
flippable, and form the reference configurations. The weights of these
configurations and the derivation of the cluster rules is illustrated in Fig 6
and the figure caption.

With these cluster rules implemented in terms of the height variables, the
super-selection sectors which have charges $Q_{x}=\pm 1$ are never generated.
However, it is still possible generate charges $Q_{x}=\pm 2$, since they are
compatible with the height representation. Therefore, in addition to the
cluster rules for the Hamiltonian, the ones for the Gauss’ Law need to
additionally implemented. To implement the Gauss’ Law note that we need to
consider four height variables meeting around a site: the top-left and the
bottom-right belong to a one sub-lattice (say A), and top-right and bottom
left belong to the other sub-lattice (B). It can be easily verified that if
the site contains a charge $\pm 2$, then both the height variables in both
sub-lattices are locally out of their reference configuration. Thus, this
situation should always be avoided. This can be implemented as follows: when
updating the A sub-lattice at a given site, it is checked if the corresponding
sites of the B sub-lattice are out of the reference configuration. If this is
the case then we bind the two height variables together, so that they are
always flipped together, and a forbidden configuration is not generated. If
the relevant height variables in the B sub-lattice are in a reference
configuration then no additional bonds are put on the A sub-lattice heights.
An identical procedure is followed while updating the other sub-lattice. It is
sufficient to check this on a single timeslice, since the Gauss’ Law commutes
with the Hamiltonian.

We emphasize that this cluster algorithm for the gauge theory is somewhat
different from the usual type of cluster algorithms since clusters are
separately built on each sub-lattice depending on the values of frozen height
variables in the other sub-lattice. In the next step, the frozen sub-lattice
is updated. This algorithm was first implemented in Bane2013  to study the
physics of the U(1) QLM and it uncovered new phases of confined gauge theories
which break translation and charge conjugation symmetry Bane2013A .

It is interesting to note that, even though this algorithm is very efficient
on the QLM, extending this to the QDM is not as useful. Since the Hamiltonian
is the same, the same cluster rules can be applied. The problem is however
with the different Gauss’ Law, which now frustrates the reference
configuration. Due to the absence of the reference configuration in the chosen
basis of the quantum height model, the clusters grow quite large to occupy
about $~{}80-90\%$ of the space-time volume, which makes the cluster algorithm
not much better than a local update Metropolis algorithm. However, a
combination of both these algorithms was used to study the physics of the QDM
Bane2014 ; Bane2015  and reliably extract the phase diagram at zero
temperature. For the QDM certain other computational methods have been
reported Oake2018 ; Yan2018 , including certain novel effective theory
approaches Herz2019 .

##  5 Conclusions and Outlook

Cluster algorithms are extremely efficient tools which are very useful to
speed up the Monte Carlo sampling of strongly interacting systems wherever
they are applicable. These algorithms act by placing bonds between degrees of
freedom which are correlated, and build up patches in configuration space
which can be flipped. The flipping process gives rise to a new configuration
which is very different from the parent one, but has the same Boltzmann
weight. Thus, the correlation in between the subsequent Monte Carlo
configurations are removed, since the clusters themselves are physical degrees
of freedom. In the case of fermions, configurations often come with negative
signs associated with fermions interchanging their respective positions. A
novel cluster algorithm — the meron algorithm, is able to analytically cancel
clusters which are responsible for the negative worldlines, and only operate
in the Hilbert space which contribute non-trivially to the partition function.
Cluster algorithms for gauge theories are even more challenging than other
systems since one has to identify proper degrees of freedom which can be used
to build clusters that do not occupy the whole space-time volume. In the case
for certain quantum link models, this has been achieved by dualizing the
original Hamiltonian, and rewriting it exactly in terms of dual quantum height
variables. The resulting cluster algorithm is very efficient. Further, for
gauge theories, the Gauss’ Law needs to be implemented, which may or may not
be possible with the cluster rules that allow the sampling the full Hilbert
space. For example, in the case of QLM where the vacuum is charge neutral,
this is not a problem, but the QDM which has a staggered background charge
this does not work so well, and the clusters become large. This can be
understood as the lack of a reference configuration for the QDM, which is an
essential ingredient for the success of a cluster algorithm.

Demanding the existence of reference configurations, it is possible to
construct models which are guaranteed to have efficient cluster algorithm.
Such Hamiltonians are often called designer Hamiltonians Kaul2013 , and can be
studied in any space-time dimensions. Since the exact form of the microscopic
interactions is unimportant to study physics associated with breaking of
symmetries across thermal or quantum phase transitions, thanks to
universality, cluster algorithms can be used to study a wide range of physical
phenomena in naturally occurring strongly correlated systems.

In the past few years, there has been renewed interest with cluster and meron
algorithms. In the context of classical spin models, the meron algorithm is
being used to study the topological charge and the vorticity properties of
non-linear sigma models Biet2019 . For lattice fermions, the meron algorithm
has been used to study a Hamiltonian spin-half lattice fermions that displays
symmetry enhancement for certain coupling regimes Liu2020 . In the context of
gauge theories, the dualization technique has been extended to a non-Abelian
$SU(2)$ quantum link model on the honeycomb lattice and an efficient cluster
algorithm has been constructed in terms of the dual height variables Bane2017
. This cluster algorithm operates on a four sub-lattice structure, where the
clusters are built using the heights on one sub-lattice, while the heights on
the other three sub-lattices decide the nature of bonds. This study identified
more crystalline confined phases in non-Abelian theories, as well as showed
fractionalization of a $Z(2)$ flux.

Another very novel and interesting development with cluster algorithms is the
simulation of strongly correlated systems in real-time and far out of
equilibrium. The first such study considered interesting initial states which
are in contact with a thermal reservoir, and are completely driven via
dissipation. The authors of Bane2014A  devised a cluster algorithm to simulate
the Lindblad evolution in an anti-ferromagnetic initial state. More extensive
studies revealed what kinds of measurements give rise to sign problems, and
which measurements did not Hebe2015 . Transport phenomena in strongly
correlated spin-$\frac{1}{2}$ driven via Lindblad dynamics was studied using
cluster algorithms in Bane2015A .

###### Acknowledgements.

This article is written in the memory of the late Pushan Majumdar. Pushan
taught me the intricacies of the multi-level algorithm Lues2001 ; Maju2003 , a
highly efficient Monte Carlo method which greatly accelerates the simulation
of pure Wilson-type gauge theories. Besides, I would like to acknowledge
Pushan for the many discussions on physics, his advice on parallel and GPU
programming that we have had, and for the many cups of coffee he would prepare
after lunch. I would like to thank Shailesh Chandrasekharan and Uwe-Jens Wiese
for teaching me so much about cluster and meron algorithms. I would also like
to acknowledge Saumen Datta, Fu-Jiun Jiang, Kieran Holland, Emilie Huffman,
Ferenc Niedermayer, Stefan Schaefer, Arnab Sen, Rainer Sommer, Urs Wenger, and
Ulli Wolff for various discussion on improved algorithms.

## References

  * (1) J. Cardy, Conformal Field Theory and Statistical Mechanics, Les Houches Summer School on Exact Methods in Low-Dimensional Statistical Physics and Quantum Computing (2008).
  * (2) A. Pelissetto, E. Vicari, Critical Phenomena and renormalization group theory, Phys. Rept., 368, 549 (2002).
  * (3) M. Moshe. J. Zinn-Justin, Quantum Field Theory in the Large-N limit: A Review, Phys. Rept., 385, 69 (2003).
  * (4) J. Maldacena, The Large-N Limit of Superconformal Field Theories and Supergravity, Intl. J. of Theo. Phys., 38, 1113 (1999).
  * (5) R. Rattazzi, V. Rychkov, E. Tonni, A. Vichi, Bounding scalar operator dimensions in 4D CFT, JHEP 12, 31 (2008).
  * (6) S. Hellerman, D. Orlando, S. Reffert, M. Watanabe, On the CFT Operator Spectrum at Large Global Charge, JHEP 17, 12 (2015).
  * (7) U. Schollwoeck, The density-matrix renormalization group in the age of matrix product states, Ann. Phys. 326, 96 (2011).
  * (8) N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, E. Teller, Equation of State Calculations by Fast Computing Machines, J. Chem. Phys. 21, 1087 (1953).
  * (9) D. Landau, K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, Cambridge University Press (2014).
  * (10) M. Newman, G. Barkema, Monte Carlo Methods in Statistical Physics, Clarendon Press (1999).
  * (11) P. W. Kasteleyn, C. M. Fortuin, Phase transitions in lattice systems with random local properties, J. Phys. Soc. Jpn. Suppl., 26, 11 (1969).
  * (12) C. M. Fortuin, P. W. Kasteleyn, On the random-cluster model. I. Introduction and relation to other models, Physica 57, 536 (1972).
  * (13) R. H. Swendsen, J. -S. Wang, Non-universal critical dynamics in Monte Carlo simulations, Phys. Rev. Lett. 58(2), 86 (1987).
  * (14) F. Niedermayer, General Cluster Updating Method for Monte Carlo Simulations, Phys. Rev. Lett. 61, 2026 (1988).
  * (15) U. Wolff, Collective Monte Carlo updating for spin systems, Phys. Rev. Lett.62(4), 361 (1989).
  * (16) F. Niedermayer, Cluster Algorithms, Advances in Computer Simulations, Lecture Notes in Physics, Vol 501, 36 (1997).
  * (17) W. Krauth, Cluster Monte Carlo Algorithms, New Optimization Algorithms in Physics, (2003).
  * (18) H. G. Evertz, G. Lana, M. Marcu, Cluster algorithm for vertex models, Phys. Rev. Lett. 70, 875 (1993).
  * (19) H. G. Evertz, M. Marcu, The Loop-Cluster Algorithm for the Case of the 6 Vertex Model, Nucl. Phys. Proc. Suppl. 30, 277 (1993).
  * (20) B. B. Beard, U.-J. Wiese, Simulations of Discrete Quantum Systems in Continuous Euclidean Time, Phys. Rev. Lett. 77, 5130 (1996).
  * (21) A. W. Sandvik, J. Kurkijärvi, Quantum Monte Carlo simulation method for spin systems, Phys. Rev. B 43, 5950 (1991).
  * (22) O. F. Syljuåsen, A. W. Sandvik, Quantum Monte Carlo with Directed Loops, Phys. Rev. E 66, 046701 (2002).
  * (23) H. G. Evertz, The loop algorithm, Advances in Physics 52 (2003).
  * (24) N. V. Prokof’ev, B. V. Svistunov, Worm Algorithms for Classical Statistical Models, Phys. Rev. Lett. 87, 160601 (2001).
  * (25) W. Heisenberg, On the theory of ferromagnetism, Zeit. für Phys., 49 (9), 619 (1928).
  * (26) M. Suzuki, Relationship between d-dimensional quantum spin system and (d+1)-dimensional Ising systems, Prog. Theor. Phys., 56, 1454 (1976).
  * (27) M. Troyer, U.-J. Wiese, Computational Complexity and Fundamental Limitations to Fermionic Quantum Monte Carlo Simulations, Phys. Rev. Lett. 94, 170201 (2005).
  * (28) R. Blankenbecler, D. J. Scalapino, R. L. Sugar, Monte Carlo calculations of coupled boson-fermion systems. I, Phys. Rev. D. 24 (8), 2278 (1981).
  * (29) F. Assaad, H. G. Evertz, World-line and determinantal quantum Monte Carlo Methods for spins, phonons and electrons, Lecture Notes in Physics, Vol. 739, pp. 277–356, Springer, Berlin Heidelberg (2008).
  * (30) S. Duane, A. D. Kennedy, B. J. Pendleton, D. Roweth, Hybrid Monte Carlo, Phys. Lett. B 195 (2), 216 (1987).
  * (31) S. Chandrasekharan, J. Cox, J. C. Osborn, U.-J. Wiese, Meron cluster approach to systems of strongly correlated electrons, Nucl. Phys. B 673, 405 (2003).
  * (32) S. Chandrasekharan, U.-J. Wiese, Meron cluster solution of fermion sign problems, Phys. Rev. Lett. 83, 3116 (1999).
  * (33) S. Chandrasekharan, The Fermion bag approach to lattice field theories, Phys. Rev. D 82, 025007 (2010).
  * (34) E. Huffman, S. Chandrasekharan, Fermion bag approach to Hamiltonian lattice field theories in continuous time, Phys. Rev. D 96, 114502 (2017).
  * (35) M. Boninsegni, N. V. Prokof’ev, B. V. Svistunov, Worm algorithm and diagrammatic Monte Carlo: a new approach to continuous-space path integral Monte Carlo simulations, Phys. Rev. E 74, 036701 (2006).
  * (36) K. van Houcke, E. Kozik, N. V. Prokof’ev, B. V. Svistunov, Diagrammatic Monte Carlo, Phys. Procedia 6, 95 (2010).
  * (37) R. Brower, P. Tamayo, Embedded dynamics for $\phi^{4}$ theory, Phys. Rev. Lett. 62, 1087 (1989).
  * (38) R. Ben-Av, H. G. Evertz, M. Marcu, S. Solomon, Critical acceleration of finite temperature SU(2) gauge simulations, Phys. Rev. D 44, 2963 (1991).
  * (39) W. Kerler, T. Metz, Ising embedding for cluster algorithms in finite temperature SU(2) gauge theory, Phys. Rev. D 50, 7082 (1994).
  * (40) D. Horn, Finite Matrix Models with Continuous Local Gauge Invariance, Phys. Lett. B 100, 149 (1981).
  * (41) P. Orland, D. Rohrlich, Lattice gauge magnets: Local isospin from spin, Nucl. Phys. B 338, 647 (1990).
  * (42) S. Chandrasekharan, U.-J. Wiese, Quantum Link Models: A Discrete Approach to Gauge Theories, Nucl. Phys. B 492, 455 (1997).
  * (43) N. Shannon, G. Misguich, K. Penc, Cyclic exchange, isolated states and spinon deconfinement in an XXZ Heisenberg model on the checkerboard lattice, Phys. Rev. B 69, 220403 (2004).
  * (44) M. Hermele, M. Fisher, L. Balents, Pyrochlore Photons: The U (1) Spin-Liquid in a S=1/2 Three Dimensional Frustrated Magnet, Phys. Rev. B 69 064404 (2004).
  * (45) L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).
  * (46) B. Schlittgen, U.-J. Wiese, Low-energy effective theories of quantum spin and quantum link models, Phys. Rev. D 63, 085007 (2001).
  * (47) J. B. Kogut, L. Susskind, Hamiltonian Formulation of Wilson’s Lattice Gauge Theories, Phys. Rev. D 11, 395 (1975).
  * (48) D. S. Rokhsar, S. A. Kivelson, Superconductivity and the Quantum Hard-Core Dimer Gas, Phys. Rev. Lett. 61, 2376 (1988).
  * (49) S. Sachdev, Spin-Peierls ground states of the quantum dimer model: A finite-size study, Phys. Rev. B 40, 5204 (1989).
  * (50) R. Moessner, K. S. Raman, Quantum Dimer Models, Lecture Notes, Trieste 2007 (arXiv:0809.3051).
  * (51) D. Banerjee, F.-J. Jiang, P. Widmer, U.-J. Wiese, The (2 + 1)-d U(1) quantum link model masquerading as deconfined criticality, J. Stat. Mech. 1312, P12010 (2013).
  * (52) D. Banerjee, P. Widmer, F.-J. Jiang, U.-J. Wiese, Crystalline Confinement, PoS LATTICE2013, 333 (2014).
  * (53) D. Banerjee, M. Bögli, C. P. Hofmann, F.-J. Jiang, P. Widmer, U.-J. Wiese, Interfaces, Strings, and a Soft Mode in the Square Lattice Quantum Dimer Model, Phys. Rev. B 90, 24, 245143 (2014).
  * (54) D. Banerjee, M. Bögli, C. P. Hofmann, F.-J. Jiang, P. Widmer, U.-J. Wiese, Finite-Volume Energy Spectrum, Fractionalized Strings, and Low-Energy Effective Field Theory for the Quantum Dimer Model on the Square Lattice, Phys. Rev. B 94, 11, 115120 (2016).
  * (55) T. Oakes, S. Powell, C. Castelnovo, A. Lamacraft, J. P. Garrahan, Phases of quantum dimers from ensembles of classical stochastic trajectories, Phys. Rev. B 98, 064302 (2018).
  * (56) Z. Yan, Y. Wu, C. Liu, O. F. Syljuåsen, J. Lou, Y. Chen, Sweeping cluster algorithm for quantum spin systems with strong geometric restrictions, Phys. Rev. B 99, 165135 (2019).
  * (57) J. Herzog-Arbeitman, S. Mantilla, I. Sodemann, Solving the quantum dimer and six vertex models one electric field line at a time, Phys. Rev. B 99, 245108 (2019).
  * (58) R. K. Kaul, R. G. Melko, A. W. Sandvik, Bridging lattice-scale physics and continuum field theory with quantum Monte Carlo simulations, Annu. Rev. Con. Mat. Phys. 4, 179 (2013).
  * (59) W. Bietenholz, J. C. Pinto Barros, S. Caspar, M. Hornung, U.-J. Wiese, Meron- and Semi-Vortex-Clusters as Physical Carriers of Topological Charge and Vorticity, PoS LATTICE2019 288, (2019).
  * (60) H. Liu, S. Chandrasekharan, R. K. Kaul, Hamiltonian models of lattice fermions solvable by the meron-cluster algorithm, arXiv:2011.13208.
  * (61) D. Banerjee, F.-J. Jiang, T. Z. Olesen, P. Orland, U.-J. Wiese, From the SU(2) quantum link model on the honeycomb lattice to the quantum dimer model on the kagomé lattice: Phase transition and fractionalized flux strings, Phys. Rev. B 97, 20, 205108 (2018).
  * (62) D. Banerjee, F.-J. Jiang, M. Kon, U.-J. Wiese, Real-Time Simulation of Large Open Quantum Spin Systems driven by Measurements, Phys. Rev. B 90, 241104 (2014).
  * (63) F. Hebenstreit, D. Banerjee, M. Hornung, F.-J. Jiang, F. Schranz, U.-J. Wiese, Real-time dynamics of open quantum spin systems driven by dissipative processes, Phys. Rev. B 92, 035116 (2015).
  * (64) D. Banerjee, F. Hebenstreit, F.-J. Jiang, U.-J. Wiese, Real-time simulation of non-equilibrium transport of magnetization in large open quantum spin systems driven by dissipation, Phys. Rev. B 92, 121104 (2015).
  * (65) M. Lüscher, Peter Weisz, Locality and exponential error reduction in numerical lattice gauge theory, JHEP 0109, 010 (2001).
  * (66) P. Majumdar, Experiences with the multilevel algorithm, Nucl. Phys. B Proc. Suppl. 119, 1021 (2003).