Quantum Monte Carlo simulations in the trimer basis: first-order transitions
and thermal critical points in frustrated trilayer magnets

L. Weber,1* A. Honecker,2 B. Normand,3,4 P. Corboz,5 F. Mila,4 S. Wessel1

1 Institute for Theoretical Solid State Physics, JARA-FIT, and JARA-HPC, RWTH
Aachen University, 52056 Aachen, Germany

2 Laboratoire de Physique Théorique et Modélisation, CNRS UMR 8089,

CY Cergy Paris Université, 95302 Cergy-Pontoise, France

3 Paul Scherrer Institute, CH-5232 Villigen-PSI, Switzerland

4 Institute of Physics, Ecole Polytechnique Fédérale de Lausanne (EPFL),
CH-1015 Lausanne, Switzerland

5 Institute for Theoretical Physics and Delta Institute for Theoretical
Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The
Netherlands

*<EMAIL_ADDRESS>
## Abstract

The phase diagrams of highly frustrated quantum spin systems can exhibit
first-order quantum phase transitions and thermal critical points even in the
absence of any long-ranged magnetic order. However, all unbiased numerical
techniques for investigating frustrated quantum magnets face significant
challenges, and for generic quantum Monte Carlo methods the challenge is the
sign problem. Here we report on a quantum Monte Carlo approach that
nevertheless allows us to study a frustrated model of coupled spin-1/2
trimers: simulations of the trilayer Heisenberg antiferromagnet in the spin-
trimer basis are sign-problem-free when the intertrimer couplings are fully
frustrated. This model features a first-order quantum phase transition, from
which a line of first-order transitions emerges at finite temperatures and
terminates in a thermal critical point. The trimer unit cell hosts an internal
degree of freedom that can be controlled to induce an extensive entropy jump
at the quantum transition, which alters the shape of the first-order line. We
explore the consequences for the thermal properties in the vicinity of the
critical point, which include profound changes in the lines of maxima defined
by the specific heat. Our findings reveal trimer quantum magnets as
fundamental systems capturing in full the complex thermal physics of the
strongly frustrated regime.

###### Contents

  1. 1 Introduction
  2. 2 Model
  3. 3 Trimer-basis QMC
    1. 3.1 Operator decomposition
    2. 3.2 Abstract directed-loop update
  4. 4 Thermodynamics of the FFTL
    1. 4.1 First-order transition
    2. 4.2 Thermal critical point
    3. 4.3 Specific heat
    4. 4.4 Correlation length and characteristic lines
  5. 5 Conclusion
  6. A Ground-state energy of the spin-3/2 Heisenberg model

##  1 Introduction

Frustrated magnets typically possess a highly degenerate low-energy subspace,
because no spin configuration can minimize every interaction term
simultaneously. This degeneracy can be lifted by collective mechanisms leading
to (quantum) order-by-disorder [1] or to a kaleidoscope of exotic disordered
phases, including quantum spin liquids [2]. While some of this physics is
captured by certain exactly solvable models [3, 4, 5], and the perturbative
regimes close to them [6], many realistic models require numerical methods for
their solution. For two-dimensional (2D) systems, these include exact
diagonalization [7, 8, 9], the density matrix renormalization group [10, 11,
12, 13], and recently developed tensor-network approaches [14, 15, 16, 17];
although both are powerful methods, each has specific drawbacks. Another
approach is the use of quantum Monte Carlo (QMC) methods, which using the
stochastic series expansion (SSE) with directed loop updates are very
efficient for quantum spin systems [18, 19, 20]. However, in the presence of
geometric frustration, QMC typically suffers severely from the negative-sign
problem [21, 22].

Although the sign problem causes an exponential decrease in the efficiency of
QMC, it is essentially basis-dependent, and for certain highly frustrated
models it has been shown that a basis exists in which QMC simulations are
exactly sign-problem-free (henceforth “sign-free”). Examples include fully
frustrated two-leg $S=1/2$ ladders in one dimension and the fully-frustrated
bilayer (FFB) [Fig. 1(a)] in two [23, 24, 25, 26, 27]. In the FFB, the basic
components are spin dimers, which are arranged on a square lattice and
connected symmetrically by nearest-neighbor Heisenberg interactions, and the
QMC sampling is performed in the local spin-dimer basis [23] (related QMC
approaches were reported in Refs. [28, 24]). Similar sign-free bases exist not
only for spin dimers but for arbitrary spin simplices, rendering
generalizations of these models accessible by QMC [23, 24]. Further, each
sign-free model is surrounded by an island in parameter space with a tolerable
sign problem, which makes QMC useful also for understanding related frustrated
quantum magnets, such as the Shastry-Sutherland lattice (SSL) model [29] by
its proximity to the FFB. In the SSL model, if the dimer coupling is not too
large, then using the dimer basis removes the sign problem completely in the
zero-temperature limit. A similar effect occurs for frustrated quantum magnets
in applied magnetic fields sufficiently far above the saturation field [30].

These recent advances in QMC methodology allow us to explore a number of
important aspects of the thermodynamic properties in specific frustrated
quantum magnets, and to trace their link to the ground-state properties. In
the example of the FFB, the spin dimers shape the ground-state phase diagram,
which is divided into a dimer-singlet quantum disordered and a dimer-triplet
antiferromagnetic (AFM) phase, separated by a first-order quantum phase
transition [31]. At finite temperature, it was observed that this transition
persists, forming a first-order line, which terminates at a critical point
[27], a scenario that has immediate parallels [32] to the liquid-gas [33],
ferromagnetic [34], and Mott transitions [35].

Figure 1: The fully frustrated bilayer (a) and trilayer models (b). In the
trilayer, the layers are labelled $1$, $2$, $3$ from bottom to top. The spins
of a trimer unit cell (inset) are coupled by intra-trimer interactions
$J_{1}$, $J_{2}$, and $J_{3}$ (red) and to all spins in the nearest-neighbor
unit cells by the same inter-trimer interaction, $J$ (black). The dashed red
line denoting the coupling ($J_{2}$) between spins in the bottom (1) and top
layers (3) of each trimer may differ from the other two couplings
($J_{1}=J_{3}$).

In this work, we study a generalization of the FFB, the spin-1/2 Heisenberg
trilayer with fully frustrated inter-trimer couplings (FFTL), illustrated in
Fig. 1(b). The FFTL constitutes an extension of the fully frustrated three-leg
ladder [36] to two dimensions, similar to the way in which the FFB is an
extension of the fully frustrated two-leg ladder. In place of the dimers
forming the basic components of the FFB, the FFTL is composed of trimer unit
cells, whose internal frustration can be varied from zero at $J_{2}=0$ to
maximal at $J_{1}=J_{2}=J_{3}$, and we will demonstate that rich physics
emerges from this internal degree of freedom. We apply the SSE QMC method,
whose sign-free basis is set by the local eigenbasis of the trimers. This
basis consists of 8 states per unit cell, in place of the 4 states in the
spin-dimer unit cell, and simulating such large local bases usually comes at
the cost of significant algorithmic complexity. To facilitate this task we
therefore present a formulation of the SSE directed-loop update scheme that
generalizes readily to arbitrary bases.

These advances allow us to explore the thermal physics of the FFTL model based
on sign-free QMC simulations. In particular, we are able to identify a first-
order transition line emerging from a first-order quantum phase transition,
which in the FFTL separates two distinct AFM ground states. Similar to the
physics of the FFB and SSL [27, 32], we show that the first-order line of the
FFTL is also terminated by a critical point, which we identify as belonging to
the 2D Ising universality class. This phenomenology shares several
similarities to the well known phase diagram of water, in which the line of
first-order transitions in the pressure-temperature plane is terminated by a
critical point belonging to the 3D Ising universality class. However, it was
noted in Ref. [32] that in the 2D quantum magnets the specific heat exhibits
two pronounced lines of maxima near the thermal critical point, in contrast to
the single line of specific-heat maxima found in water. This behavior can be
traced to the special case of the Ising model, whose $Z_{2}$ symmetry in fact
enforces an exact symmetry between the two specific-heat lines in the
magnetic-field-temperature phase diagram [32]. As we report below, varying the
intra-trimer frustration in the unit cell allows us to (i) control the
relative strengths of the two lines of maxima in the specific heat of the
FFTL, and thus (ii) connect the generic “Ising-like” behavior observed to date
in 2D quantum magnets to a more “water-like” type of behavior, featuring a
single line of pronounced specific-heat maxima. Further, we connect the
strongest water-like behavior, which occurs when the intra-trimer frustration
is maximal, to an extensive entropy jump across the first-order quantum
transition and hence an enhanced “slanting” of the first-order transition line
in the plane of “pressure” (coupling ratio) and temperature. Our observations
thus establish the FFTL as a generic quantum spin model in which to explore
the full complexity of the thermal physics uncovered in Refs. [32, 27], by the
unbiased numerical technique of sign-free QMC simulations. Taken together,
these provide us with the basis on which to understand the differences between
the superficially distinct forms of thermodynamic behavior in certain quantum
magnets and in water.

This article is structured as follows. In Sec. 2, we provide an overview of
the FFTL model. In Sec. 3 we outline the QMC method for simulating the FFTL in
the trimer basis. This allows us to present in Sec. 4 our findings for the
thermal physics of the FFTL, discussing in particular the first-order line,
the thermal critical point, the behavior of the specific heat and of other
quantities characterizing the critical fluctuations, and the effect of intra-
trimer frustration on these properties. We conclude in Sec. 5 by discussing
the implications of our results and their relation to other recent studies of
thermal critical points in frustrated quantum magnets.

##  2 Model

In analogy to the well known FFB [37, 38, 39, 40, 41, 24, 26, 27] and the
fully frustrated three-leg ladder [36], the FFTL is a square lattice of
$S=1/2$ spin trimers where all the sites in nearest-neighbor unit cells are
connected by the same interaction, $J$, giving the Hamiltonian

$\displaystyle H$
$\displaystyle=\sum_{\Delta}H_{\Delta}+H_{\Delta,\Delta+\hat{x}}+H_{\Delta,\Delta+\hat{y}},$
$\displaystyle H_{\Delta}$
$\displaystyle=J_{1}\,\mathbf{S}_{\Delta,2}\cdot\mathbf{S}_{\Delta,3}+J_{2}\,\mathbf{S}_{\Delta,3}\cdot\mathbf{S}_{\Delta,1}+J_{3}\,\mathbf{S}_{\Delta,1}\cdot\mathbf{S}_{\Delta,2},$
$\displaystyle H_{\Delta,\Delta^{\prime}}$
$\displaystyle=J\sum_{i,j=1}^{3}\mathbf{S}_{\Delta,i}\cdot\mathbf{S}_{\Delta^{\prime},j}=J\,\mathbf{S}_{\Delta}\cdot\mathbf{S}_{\Delta^{\prime}}.$
(1)

The sum over $\Delta$ enumerates the trimer unit cells,
$\mathbf{S}_{\Delta,i}$ is the $i$th spin of trimer $\Delta$, which belongs to
the $i$th layer in Fig. 1(b), and
$\mathbf{S}_{\Delta}=\mathbf{S}_{\Delta,1}+\mathbf{S}_{\Delta,2}+\mathbf{S}_{\Delta,3}$
is the total-spin operator of trimer $\Delta$. The unit cells are arranged as
a square superlattice of $L\times L$ trimers.

Figure 2: Spectrum of a single trimer with spin rotational symmetry. The 8
states belong to one quartet, Q, and two doublets, D$0$ and D$1$. The doublets
are degenerate if and only if $J_{1}=J_{2}=J_{3}$.

The spectrum of a single trimer consists of 8 states that belong, due to the
spin rotational symmetry of $H_{\Delta}$, to one quartet (denoted
$\mathrm{Q}$) or two doublets ($\mathrm{D}0$ and $\mathrm{D}1$), with
respective eigenenergies $\varepsilon_{\mathrm{Q}}$,
$\varepsilon_{\mathrm{D}0}$, or $\varepsilon_{\mathrm{D}1}$ (Fig. 2). If
$J_{1}=J_{2}=J_{3}$, the two doublets are degenerate, resulting for AFM
interactions in a fourfold ground-state degeneracy. This is a direct
consequence of the symmetry group $G=\mathrm{C}_{3\mathrm{v}}\times$SU(2) at
the degeneracy point: the irreducible representations of $G$ are tensor
products of the irreducible representations of the two factors, and one
obtains in this case a decomposition of the trimer Hilbert space into
$\mathrm{A}_{1}\otimes\mathbf{4}\>\oplus\>\mathrm{E}\otimes\mathbf{2}$. Away
from this point, the doublet levels split and the ground state of a single
trimer is only doubly degenerate. For $J_{1}=J_{3}$, the eigenstates can be
classified completely by the expectation values of the commuting operators

$\displaystyle\mathbf{S}^{2}_{\Delta}$
$\displaystyle=l_{\Delta}(l_{\Delta}+1),$ (2)
$\displaystyle(\mathbf{S}_{\Delta,1}+\mathbf{S}_{\Delta,3})^{2}$
$\displaystyle=l_{\Delta,13}(l_{\Delta,13}+1),$ (3) $\displaystyle
S^{z}_{\Delta}$ $\displaystyle=m_{\Delta},$ (4)

and for convenience we restrict ourselves henceforth to this symmetric case
(Fig. 1).

In the full Hamiltonian, the fully frustrated inter-trimer interaction, $J$,
couples the total spins, $\mathbf{S}_{\Delta}$, of adjacent trimers. Because

$\displaystyle[\mathbf{S}^{2}_{\Delta},S^{\alpha}_{\Delta}]$
$\displaystyle=0,$ (5)
$\displaystyle[(\mathbf{S}_{\Delta,1}+\mathbf{S}_{\Delta,3})^{2},S^{\alpha}_{\Delta}]$
$\displaystyle=0,$ (6) $\displaystyle[S^{z}_{\Delta},S^{\alpha}_{\Delta}]$
$\displaystyle=i\epsilon^{z\alpha\beta}S^{\beta}_{\Delta}\neq 0,$ (7)

$l_{\Delta}$ and $l_{\Delta,13}$ remain good (local) quantum numbers of the
full model. However, $m_{\Delta}$ is no longer conserved.

A notable consequence of Eqs. (6) and (7) is that the matrix elements of the
fully frustrated inter-trimer interactions in the “trimer” basis labelled by
$l_{\Delta}$, $l_{\Delta,13}$, and $m_{\Delta}$ take the form

$\braket{l_{\Delta},l_{\Delta,13},m_{\Delta}}{\mathbf{S}_{\Delta}}{l^{\prime}_{\Delta},l^{\prime}_{\Delta,13},m^{\prime}_{\Delta}}=\delta_{l_{\Delta},l_{\Delta}^{\prime}}\delta_{l_{\Delta,13},l_{\Delta,13}^{\prime}}\braket{m_{\Delta}}{\mathbf{S}^{S=l_{\Delta}}_{\Delta}}{m_{\Delta}^{\prime}},$
(8)

where $\mathbf{S}^{S=l_{\Delta}}_{\Delta}$ is a spin-$l_{\Delta}$ spin
operator. The matrix elements therefore take the same values in the
$S_{\Delta}^{z}$ basis as in an effective square-lattice AFM of mixed $S=1/2$
and $S=3/2$. The spin-$S$ Heisenberg AFM on a bipartite lattice can, however,
be rendered sign-free [42], and consequently the inter-trimer Hamiltonian,
$H_{\Delta,\Delta^{\prime}}$, of the FFTL in the trimer basis,
$\ket{l_{\Delta},l_{\Delta,13},m_{\Delta}}$, is also sign-free. The remaining
intra-trimer part, $H_{\Delta}$, is diagonal in this basis and is thus sign-
free as well. The trimer basis thus enables an efficient QMC algorithm to
simulate the FFTL model, which we present next.

##  3 Trimer-basis QMC

An efficient means of simulating quantum spin systems is offered by the SSE
QMC method, which is based on expressing the quantum Boltzmann operator as an
infinite-order high-temperature expansion with easily evaluated, positive
coefficients [18]. This representation, which can be interpreted as a
classical probability distribution, is then supplemented by the global
directed-loop [19] class of updates that allow an efficient sampling of the
probability distribution using the Markov-chain Monte Carlo approach.

###  3.1 Operator decomposition

The starting point for the SSE is the decomposition of the Hamiltonian into
non-branching operators (those which map each basis vector to a single basis
vector, and not to a superposition of different basis vectors). We express the
FFTL model Hamiltonian in the form

$H=\sum_{\braket{\Delta,\Delta^{\prime}}}\tilde{H}_{\Delta,\Delta^{\prime}},\quad\tilde{H}_{\Delta,\Delta^{\prime}}=H_{\Delta,\Delta^{\prime}}+\frac{1}{4}\left(H_{\Delta}+H_{\Delta^{\prime}}\right),$

where the intra-trimer terms have been absorbed into the inter-trimer part,
and decompose $H$ into the operators

$\displaystyle H$
$\displaystyle=\sum_{\begin{subarray}{c}\braket{\Delta,\Delta^{\prime}}\\\
x_{\Delta},x_{\Delta^{\prime}}\\\
y_{\Delta},y_{\Delta^{\prime}}\end{subarray}}\braket{x_{\Delta},x_{\Delta^{\prime}}}{\tilde{H}_{\Delta,\Delta^{\prime}}}{y_{\Delta},y_{\Delta^{\prime}}}\ket{x_{\Delta},x_{\Delta^{\prime}}}\\!\bra{y_{\Delta},y_{\Delta^{\prime}}}=:\sum_{\alpha}h_{\alpha},$
(9)

denoting the sets of trimer-basis states
$\ket{l_{\Delta},l_{\Delta,13},m_{\Delta}}$ for brevity by $\ket{x_{\Delta}}$
and $\ket{y_{\Delta}}$. The operators $h_{\alpha}$ of this decomposition are
labelled by the bond $\braket{\Delta,\Delta^{\prime}}$ and for notational
simplicity we combine the summation indices into a single multi-index,
$\alpha$. Using this decomposition, the partition function can be expanded
into a form that is simple to evaluate because it contains only sums over non-
branching operator strings (details may be found in Ref. [42]),

$\displaystyle Z=\operatorname{Tr}e^{-\beta H}$
$\displaystyle=\sum_{n=0}^{\infty}\sum_{\psi_{0}}\frac{\beta^{n}}{n!}\braket{\psi_{0}}{(-H)^{n}}{\psi_{0}}$
(10)
$\displaystyle=\sum_{n=0}^{\infty}\sum_{\psi_{0}}\sum_{\\{\alpha_{p}\\}_{p=1}^{n}}\frac{\beta^{n}}{n!}\braket{\psi_{n}}{\psi_{0}}\prod_{p=1}^{n}\braket{\psi_{p-1}}{(-h_{\alpha_{p}})}{\psi_{p}},$
(11)

where $\ket{\psi}=\ket{x_{1},x_{2},x_{3},\cdots}$ is the trimer-product basis
of the full Hilbert space and the sum is over all $n$-sequences of $\alpha$.
The sequence of intermediate states $\\{\psi_{1},\dots,\psi_{n}\\}$ is
generated by the successive application of the non-branching operators
$h_{\alpha}$ to the trace state $\ket{\psi_{0}}$.

To increase numerical efficiency, the sequences of $h_{\alpha}$ are in
practice padded by identity operators so that they always have a fixed length,
$\mathcal{L}$. This length is chosen adaptively to be sufficiently large that
it has no measurable effect on the result. However, the padding procedure
introduces a combinatorial factor that must be accounted for in the padded
expansion (so that $n$ now denotes the number of non-identity operators
contained in the operator string $\\{\alpha_{p}\\}$), leading to

$\displaystyle\tilde{Z}$
$\displaystyle=\sum_{\psi_{0},\\{\alpha_{p}\\}_{p=1}^{\mathcal{L}}}\frac{\beta^{n}(\mathcal{L}-n)!}{\mathcal{L}!}\braket{\psi_{\mathcal{L}}}{\psi_{0}}\prod_{p=1}^{\mathcal{L}}\braket{\psi_{p-1}}{(-h_{\alpha_{p}})}{\psi_{p}}=:\sum_{\mathcal{C}}W(\mathcal{C}).$
(12)

Because the weights $W(\mathcal{C})$ are non-negative for the FFTL in the
trimer basis, the quantum partition function can be interpreted as a classical
probability distribution and different observables can be calculated using
Markov-chain Monte Carlo methods.

###  3.2 Abstract directed-loop update

To sample from a Markov chain, we need an ergodic update for the configuration
$\mathcal{C}=\\{\psi_{0},\\{\alpha_{p}\\}\\}$ that fulfills the detailed-
balance criterion. A combination of two types of update is used widely in the
context of the SSE. First, the diagonal update inserts or removes diagonal
operators (in the computational basis) locally in the string
$\\{\alpha_{p}\\}$ [18]. Second, for the off-diagonal operators and the state
$\psi_{0}$, global loop updates respecting the trace boundary condition on
$\braket{\psi_{\mathcal{L}}}{\psi_{0}}\neq 0$ are employed [19, 43]. In
general, directed-loop updates function as follows: starting from a given SSE
configuration, two discontinuities are introduced at lattice site $i$ with
operator index $p$. Referring to these discontinuities as the loop head and
tail, the head is propagated stochastically through the SSE operator string
until it meets with the tail. At this point the two discontinuities can
annihilate each other again, leaving a loop along which the spin configuration
has changed.

The diagonal updates generalize readily to arbitrary bases, and so can be
applied directly. For the directed-loop updates, however, the situation is
more complicated. Here, the behavior of the loops is typically hand-crafted to
a specific model by associating the possible discontinuities with physically
intuitive operators. As an example, the operators $S^{+}$ and $S^{-}$ are used
for spin systems in the conventional single-site $S^{z}$ basis. For local
bases containing more than a single site, making this identification becomes
increasingly difficult, because the number of basis states grows faster than
the number of intuitive operators. An example from the dimer basis is the
inclusion of the spin-difference operators,
$\mathbf{D}=\mathbf{S}_{1}-\mathbf{S}_{2}$, in addition to normal spin flips
[23, 25]. For the trimer basis and other large local bases, following this
route would lead to even more complexity.

As an alternative we outline an abstract directed-loop update, which can be
constructed canonically for any basis or model. As a general starting point,
we associate the discontinuities at the loop head and tail not with operators
but with bijective functions or “actions,”

$a:X\rightarrow Y;\quad X,Y\subset\mathcal{B},$ (13)

which map the state before the discontinuity, $x$, to the state after the
discontinuity, $a(x)$, both states being elements of the local (trimer or
other arbitrary) computational basis, $\mathcal{B}$. In this picture, the
equivalent of $S^{+}$ and $S^{-}$ for the single-spin $S^{z}$ basis would be
the two actions $a^{\pm}:\ket{m}\mapsto\ket{m\pm 1}$, where $m$ denotes the
$S^{z}$ eigenvalue. Note that the domain $X$ and image $Y$ of an action $a$
may not include all of the basis states, as is the case here for $a^{+}$ and
the state $\ket{m=S}$.

Figure 3: Abstract directed-loop updates. Whenever the loop head arrives at
leg $l_{\text{in}}$ of an operator $h_{\alpha}$ carrying the action
$a_{\text{in}}$, it may continue on a different leg, $l_{\text{out}}$, with a
different action, $a_{\text{out}}$. Along the path of the head, the four
states $x_{i}$, $x_{j}$, $y_{i}$, and $y_{j}$ defining $h_{\alpha}$ are
changed by the application of $a_{\text{in}}$ and $a_{\text{out}}$.

The loop head moves through the configuration carrying one of these actions,
and whenever it hits an operator $h_{\alpha}$ there are different ways to
proceed, as represented in Fig. 3. The head can change its real-space position
and direction by hopping to a different operator leg, $l$, or it can change
the type of its discontinuity, meaning its action $a$. Depending on these
choices, the operator is changed, $h_{\alpha}\rightarrow h_{\alpha^{\prime}}$,
after the head has passed. A way to bring this change of configuration into
agreement with detailed balance is to demand that each junction, or scattering
event, within the propagation process itself fulfills detailed balance.

Adopting the notation of actions $a$ and legs $l$, each operator scattering
event can be described as a transition

$\alpha,a_{\text{in}},l_{\text{in}}\rightarrow\alpha^{\prime},a_{\text{out}},l_{\text{out}}.$
(14)

If the operator before scattering was
$h_{\alpha}=-w_{\alpha}\ket{x_{1},x_{2}}\\!\bra{x_{3},x_{4}}$, the new
operator, $h_{\alpha^{\prime}}$, is the one where first $x_{l_{\text{in}}}$ is
replaced by $x^{\prime}_{l_{\text{in}}}=a_{\text{in}}(x_{l_{\text{in}}})$, and
in a second step the state at $l_{\text{out}}$ is set to
$x^{\prime\prime}_{l_{\text{out}}}=a_{\text{out}}(x^{\prime}_{l_{\text{out}}})$.
The detailed-balance condition for this process is then

$w_{\alpha}\,p(\alpha,a_{\text{in}},l_{\text{in}}\rightarrow\alpha^{\prime},a_{\text{out}},l_{\text{out}})=w_{\alpha^{\prime}}\,p(\alpha^{\prime},a_{\text{out}}^{-1},l_{\text{out}}\rightarrow\alpha,a_{\text{in}}^{-1},l_{\text{in}}).$
(15)

Taken together, this condition and the normalization conditions for the
transition probabilities produce a system of equations from which not all
probabilities are determined uniquely. In practice, one typically chooses an
optimal solution according to some heuristic criteria, such as the
minimization of “bounce” processes (those for which
$l_{\text{out}}=l_{\text{in}},~{}a_{\text{out}}=a_{\text{in}}^{-1}$, Fig.
3(a)) [19]. Finally, it is necessary to specify how the loop ends. If the
head, carrying the action $a_{h}$, meets the tail, carrying the action
$a_{t}$, and $a_{h}=a_{t}^{-1}$, the loop terminates. Otherwise, $a_{t}$ is
replaced by the composition $a_{h}\circ a_{t}$.

So far, everything we have done is equivalent to the hand-crafted approach,
with a slight change of nomenclature. The question for a given model is always
how to choose the set of actions $\\{a_{k}\\}$ that the loop head can carry.
In a sense, the most complete set of actions, $\\{a_{k},~{}k=1,...,d-1\\}$,
that can be applied to a local basis, $\mathcal{B}=\\{x_{1},\dots,x_{d}\\}$
fulfills the property

$\\{a_{k}(x),~{}k=1,\dots,d-1\\}=\mathcal{B}\backslash\\{x\\}\qquad\forall
x\in\mathcal{B},$ (16)

i.e., for all basis states $x\neq y$ there exists an action $a_{k}$ so that
$a_{k}(x)=y$. In other words, the discontinuity at a loop head can change
every state to every other state, given the right loop action. The identity
transitions $a_{k}(x)=x$ can be excluded because they are in general not
useful in changing the configuration.

Adding more actions to a complete set does not change the behavior of the
loops, because all state discontinuities are already possible. Using a
minimally complete set means that all hand-crafted implementations with fewer
operators are included automatically as a solution in Eq. (15). Such an
“incomplete” set of actions would not allow all discontinuities but does
fulfill detailed balance; a solution equivalent to the incomplete set can be
obtained within a complete set of actions by setting the probabilities leading
to the “forbidden” discontinuities to zero. Thus for the abstract directed-
loop update we use a complete set of actions.

There are clearly many choices for a set of functions fulfilling the
completeness property from Eq. (16). However, because all possible state
discontinuities are already allowed, the exact implementation of the set does
not matter. Changing from one implementation to another merely constitutes a
relabelling of the transition indices
$\alpha,a_{\text{in}},l_{\text{in}}\rightarrow\alpha^{\prime},a_{\text{out}},l_{\text{out}}$.
As a result, we can choose a computationally convenient set of functions, such
as

$a^{\text{xor}}_{k}:x_{l}\mapsto x_{\text{xor}(k,l)},\qquad k=1,\dots,d-1,$
(17)

where $\text{xor}(k,l)$ is the bitwise exclusive-or function of the integers
$k$ and $l$.

One may ask how allowing all discontinuities, as opposed to allowing only a
selected subset, impacts computational performance. In the example of an $S=1$
Heisenberg model, the complete set of actions allows the head to interchange
the states $\ket{+1}$ and $\ket{-1}$, which cannot be expressed using just the
actions $a^{+}$ and $a^{-}$. It is well known that including $(S^{\pm})^{2}$
loops, which would correspond to such discontinuities, does not increase the
efficiency because of the $S^{+}_{i}S^{-}_{j}$ structure of the bond
Hamiltonian. Starting from a complete set of actions, one would therefore need
to go beyond the bounce-minimization heuristic for optimal performance and
search for solutions of Eq. (15) that discourage magnetization changes $\Delta
m\neq\pm 1$ explicitly. Using such a solution, the behavior of the abstract
directed-loop update is again equivalent to the standard $S^{\pm}$ one, other
than some performance cost associated with larger data structures.

The advantage of the abstract directed-loop update is in models and bases
where knowledge about the optimal set of allowed discontinuities does not yet
exist, such as the FFTL in the trimer basis. The abstract directed-loop update
instead provides a canonical starting point, namely a complete set of actions
with the bounce-minimization heuristic. This starting point is equivalent to
any set of hand-crafted operators that allow all discontinuities and, in
principle, can be optimized systematically.

In our implementation, we use the $a^{\text{xor}}$ actions and linear
programming to minimize the probability of bounces in the solution of Eq.
(15). While it is not always possible to avoid bounces completely, this
heuristic alone was sufficient to achieve acceptable performance for the FFTL.
More specifically, we were able to perform efficient simulations for systems
of $L\times L$ trimers up to $L=48$, at temperatures down to $T=0.1J_{1}$. For
scans of the full $T$-$J$ plane, we restricted our simulations to systems with
$L=12$, this being sufficient to access the physics away from the critical
point; close to the critical point, we performed finite-size scaling to
capture the effect of the diverging correlation length.

##  4 Thermodynamics of the FFTL

We now present the results obtained by applying the abstract directed-loop QMC
approach to examine the thermal properties of the FFTL. Phase diagrams in the
plane of temperature and the ratio of inter- to intra-trimer coupling are
presented in Figs. 4(a-d) for $J_{1}=J_{3}$ with 4 different values of
$J_{2}/J_{1}$. The phase diagram is obtained from the expectation value of the
total spin, $\mathbf{S}^{2}_{\Delta}=l_{\Delta}(l_{\Delta}+1)$, on each trimer
and the color scale indicates the mean of its associated quantum number,
$\braket{l_{\Delta}}$.

Figure 4: (a)-(d) Triangle total-spin quantum number, $\langle
l_{\Delta}\rangle$, computed for the $L=12$ FFTL model at different values of
the intra-trimer coupling ratio $J_{2}/J_{1}$, with $J_{1}=J_{3}$. Red lines
show the lines of first-order transitions at finite temperature, calculated
from the approximation developed in Subsec. 4.1, and red circles mark the
critical point determined numerically in Subsec. 4.2. (e) Comparison of the
shapes of the first-order lines obtained for different $J_{2}/J_{1}$ values
from the approximation of Subsec. 4.1.

At low temperatures, $\braket{l_{\Delta}}$ approaches the value $1/2$ if the
AFM intra-trimer couplings dominate, but sufficiently strong $J$ overrides
this preference, optimizing the inter-trimer interaction terms by creating a
ferromagnetic configuration on each trimer, causing $\braket{l_{\Delta}}$ to
approach $3/2$. These effective $S=1/2$ and $S=3/2$ AFM ground states are
lower in energy than any states of mixed $l_{\Delta}$, because the fully
frustrated inter-trimer interactions favor neighboring trimers with the same
$l_{\Delta}$. Thus they are connected directly by a level-crossing, or first-
order quantum phase transition, occurring at one $J_{2}$-independent value of
the inter-to-intra-trimer coupling ratio, $J_{c}/J_{1}$.

At low finite temperatures, this first-order transition persists, forming a
line in the phase diagram. However, sufficiently high temperatures cause the
transition to turn into a crossover, and hence each line in Figs. 4(a-d)
terminates at a critical point. This behavior can be associated with the
effects of thermal fluctuations into mixed-$l_{\Delta}$ states and we defer a
detailed determination of the critical points to Sec. 4.2.

At a superficial level these results resemble the phase diagram of the FFB,
but there are two major distinctions. First, the first-order transition takes
place between two phases that both have a broken symmetry at $T=0$. This is in
contrast to the quantum disordered dimer-singlet ground state on one side of
the transition in the FFB. Second, in the FFTL there is another independent
control parameter in addition to the inter-trimer coupling, and it is clear
from Figs. 4(a-d) that this $J_{2}/J_{1}$ parameter influences both the shape
of the first-order line and the location of the critical point.

To investigate the physics of the FFTL phase diagram, we will first introduce
some low-temperature approximations appropriate for gaining an analytical
understanding of the first-order transition line (Subsec. 4.1). Because this
approach does not allow us to extract the position of the critical point, in
Subsec. 4.2 we perform detailed finite-size scaling to identify the critical
point numerically and to discuss the influence of $J_{2}$ on the critical
temperature. In Subsec. 4.3, we examine the specific heat and discuss its
evolution with the shape of the first-order line. In Subsec. 4.4, we
supplement this analysis by presenting results for the correlation length
around the critical point, which allow us to compare its line of maxima with
that of the specific heat.

###  4.1 First-order transition

With the goal of finding a straightforward means of explaining the shape of
the first-order transition line in the low-temperature region, we use the
locally conserved quantities of the FFTL model to obtain low-temperature
approximations to the free energy on both sides of the transition, an approach
used for the FFB in Ref. [27]. We start by rewriting the partition function
using the local quantum numbers,
$\\{n_{\Delta}\\}:=\\{(l_{\Delta},l_{\Delta,13})\\}$, to obtain

$\displaystyle Z$
$\displaystyle=\operatorname{Tr}\exp\Big{[}-\beta\big{(}\sum_{\Delta}H_{\Delta}+\underbrace{\sum_{\Delta}H_{\Delta,\Delta+\hat{x}}+H_{\Delta,\Delta+\hat{y}}}_{=:H_{J}}\big{)}\Big{]}$
$\displaystyle=\sum_{\\{n_{\Delta}\\}}e^{-\beta\sum_{\Delta}\varepsilon({n_{\Delta}})}\operatorname{Tr}P_{\\{n_{\Delta}\\}}e^{-\beta
H_{J}}=\sum_{\\{n_{\Delta}\\}}\operatorname{Tr}e^{-\beta
H(\\{n_{\Delta}\\})}=:\sum_{\\{n_{\Delta}\\}}e^{-\beta F(\\{n_{\Delta}\\})},$
(18)

where $P_{\\{n_{\Delta}\\}}$ is a projection onto the subspace of the set of
local trimer quantum numbers, $\\{n_{\Delta}\\}$, and

$H(\\{n_{\Delta}\\})=\sum_{\Delta}\varepsilon(n_{\Delta})+P_{\\{n_{\Delta}\\}}H_{J}P_{\\{n_{\Delta}\\}}.$
(19)

Thus, $F(\\{n_{\Delta}\\})$ is the free energy of a system with a frozen
configuration of trimer quantum numbers belonging to the partition function

$Z(\\{n_{\Delta}\\})=\operatorname{Tr}e^{-\beta H(\\{n_{\Delta}\\})}.$ (20)

In this way, the FFTL model can be regarded as a mixture of all possible ways
to combine $S=1/2$ and $S=3/2$ trimer spins into a square-lattice AFM. Because
$F(\\{n_{\Delta}\\})$ is an extensive quantity, the distribution in Eq. (18)
will generally become sharp in the thermodynamic limit and at low temperature,
and the system will choose the configuration $\\{n_{\Delta}\\}$ that minimizes
$F(\\{n_{\Delta}\\})$.

Our first approximation is that $F(\\{n_{\Delta}\\})$ is minimized by a
uniform configuration, $n_{\Delta}=n$. At low temperatures, this is expected
because a system containing multiple domains of different $n_{\Delta}$ can
decrease its energy by growing the single domain with the lowest energy per
site. At high temperatures, however, “domain walls” become favorable due to
their entropy, and lead to the appearance of the critical point. By neglecting
domain walls, our approximation loses the ability to predict the end of the
first-order line. Further, it is conceivable that domain walls between
$l_{\Delta}=1/2$ and $l_{\Delta}=3/2$ could become energetically favorable
even at $T=0$, due to quantum fluctuation effects, although we will exclude
this possibility a posteriori for the bare FFTL by comparing our approximate
results to the QMC data.

For simplicity, we first consider $T=0$, where $F(\\{n_{\Delta}\\})$ reduces
to the ground-state energy, $E_{0}(\\{n_{\Delta}\\})$, of
$H(\\{n_{\Delta}\\})$. As explained above, we assume

$\min_{\\{n_{\Delta}\\}}E_{0}(\\{n_{\Delta}\\})=E_{0}(n)$ (21)

and restrict our considerations to three possible ground-state configurations
with uniform $n_{\Delta}$, having total energies per trimer ($N_{t}$ denotes
the number of trimers)

$\displaystyle E_{\mathrm{Q}}/N_{t}$
$\displaystyle=\varepsilon_{\mathrm{Q}}+J\varepsilon_{\text{AFM}}^{S=3/2},$
(22) $\displaystyle E_{\mathrm{D}1}/N_{t}$
$\displaystyle=\varepsilon_{\mathrm{D}1}+J\varepsilon_{\text{AFM}}^{S=1/2},$
(23) $\displaystyle E_{\mathrm{D}0}/N_{t}$
$\displaystyle=\varepsilon_{\mathrm{D}0}+J\varepsilon_{\text{AFM}}^{S=1/2},$
(24)

where $\varepsilon_{\text{AFM}}^{S=1/2}$ ($\varepsilon_{\text{AFM}}^{S=3/2}$)
is the ground-state energy per site of the spin-1/2 (spin-3/2) square-lattice
Heisenberg AFM with unit couplings. We denote by $\mathrm{D0}$ the doublet
with lower energy, $\varepsilon_{\mathrm{D}0}\leq\varepsilon_{\mathrm{D}1}$,
to obtain for the ground state

$\displaystyle n=\begin{cases}\mathrm{D}0&\text{if
}\varepsilon_{\mathrm{D}0}+J\varepsilon^{S=1/2}_{\text{AFM}}\leq\varepsilon_{\mathrm{Q}}+J\varepsilon^{S=3/2}_{\text{AFM}},\\\
\mathrm{Q}&\text{otherwise.}\end{cases}$ (25)

The critical inter-trimer coupling where the two levels cross is

$J_{c}=\frac{\varepsilon_{\mathrm{Q}}-\varepsilon_{\mathrm{D}0}}{\varepsilon^{S=1/2}_{\text{AFM}}-\varepsilon^{S=3/2}_{\text{AFM}}},$
(26)

and because $J_{c}>0$ a first-order quantum phase transition takes place
between a spin-1/2 and a spin-3/2 AFM phase. The value of $J_{c}$ does not
depend on $J_{2}$ because $l_{\Delta,13}=1$ in both the Q and D0 levels, and
the contribution from the $J_{2}$ dimer cancels in the numerator of Eq. (26).

The reasoning to this point is analogous to the FFB case, with the exception
of the additional quantum number, $l_{\Delta,13}$, that labels the doublets
$\mathrm{D}0$ and $\mathrm{D}1$. Although the $\mathrm{D}1$ level has no
direct influence at $T=0$, this changes at finite temperatures. When $T>0$,
the first-order transition is expected to continue along a line, $J_{c}(T)$,
up to the critical point at $(T_{c},J_{c}(T_{c}))$. For $T\ll T_{c}$, we
extend the approximation to finite temperatures by replacing energies with
free energies. For non-degenerate $\mathrm{D}0$ and $\mathrm{D}1$, the system
will adopt the state with the lower free energy of

$\displaystyle F_{\mathrm{Q}}/N_{t}$
$\displaystyle=\varepsilon_{\mathrm{Q}}+J\varepsilon_{\text{AFM}}^{S=3/2}(T/J)-Ts_{\text{AFM}}^{S=3/2}(T/J),$
(27) $\displaystyle F_{\mathrm{D}0}/N_{t}$
$\displaystyle=\varepsilon_{\mathrm{D}0}+J\varepsilon_{\text{AFM}}^{S=1/2}(T/J)-Ts_{\text{AFM}}^{S=1/2}(T/J).$
(28)

However, if $\epsilon_{\mathrm{D}1}-\epsilon_{\mathrm{D}0}\lesssim T$, it is
not meaningful to consider the configuration $n_{\Delta}=\mathrm{D}0$ as a
macroscopic state, because thermal fluctuations allow admixtures of
$\mathrm{D}0$ and $\mathrm{D}1$ to acquire a finite weight. Instead we subsume
all configurations with $n_{\Delta}\in\\{\mathrm{D}0,\mathrm{D}1\\}$ into a
single macrostate,

$\displaystyle e^{-\beta F_{\mathrm{D}}}$
$\displaystyle=\sum_{\begin{subarray}{c}\\{n_{\Delta}\\}\\\
n_{\Delta}=\mathrm{D}0,\mathrm{D}1\end{subarray}}e^{-\beta[\varepsilon_{n}+F(\\{n_{\Delta}\\})]}=\left(1+e^{-\beta(\varepsilon_{\mathrm{D}1}-\varepsilon_{\mathrm{D}0})}\right)^{N_{t}}e^{-\beta
F_{\mathrm{D}0}},$ (29)

with free energy

$\displaystyle F_{\mathrm{D}}$
$\displaystyle=F_{\mathrm{D}0}-TS_{\mathrm{D}},$ (30) $\displaystyle
S_{\mathrm{D}}$
$\displaystyle=N_{t}\log\left(1+e^{-\beta(\varepsilon_{\mathrm{D}1}-\varepsilon_{\mathrm{D}0})}\right).$
(31)

As expected, the additional entropy contribution due to the doublets vanishes
when $\varepsilon_{\mathrm{D}1}-\varepsilon_{\mathrm{D}0}\gg T$ and assumes
the maximal value of $N_{t}\log 2$ in the degenerate case
$\varepsilon_{\mathrm{D}1}=\varepsilon_{\mathrm{D}0}$ ($J_{1}=J_{2}=J_{3}$).

These expressions allow us to estimate the first-order line, $J_{c}(T)$, for
$T\ll T_{c}$ in a manner similar to Ref. [27]. Because $F_{\mathrm{Q}}$ and
$F_{\mathrm{D}0}$ describe square-lattice Heisenberg AFMs, we employ the
standard low-temperature approximation for the specific heat, $C(T,J)=\partial
E(T,J)/\partial T/N_{t}$, to capture the behavior of the free energy at low
temperatures,

$\displaystyle F(T,J)$
$\displaystyle=E(0,J)-TS(0,J)+\int_{0}^{T}dT^{\prime}~{}\left(1-\frac{T}{T^{\prime}}\right)N_{t}C(T^{\prime},J).$
(32)

For the spin-$S$ AFM, the spin-wave excitations (magnons) give a power-law
contribution to the specific heat [44, 45, 46, 47],

$\displaystyle C(T,J)$ $\displaystyle=\frac{\partial
J\varepsilon_{\text{AFM}}(T/J)}{\partial
T}=a(S)\left(\frac{T}{J}\right)^{2}+b(S)\left(\frac{T}{J}\right)^{4}+\dots,$
(33)

where we have also included the next-to-leading-order correction. The free
energy then becomes

$F(T)=E(0)-TS(0)-\frac{a(S)}{6J^{2}}T^{3}-\frac{b(S)}{20J^{4}}T^{5}+\dots,$
(34)

or for the FFTL model

$\displaystyle F_{\mathrm{Q}}(T,J)/N_{t}$
$\displaystyle=\epsilon_{\mathrm{Q}}+J\varepsilon_{\text{AFM}}^{S=3/2}-\frac{a(3/2)}{6J^{2}}T^{3}-\frac{b(3/2)}{20J^{4}}T^{5},$
(35) $\displaystyle F_{\mathrm{D}0}(T,J)/N_{t}$
$\displaystyle=\epsilon_{\mathrm{D}0}+J\varepsilon_{\text{AFM}}^{S=1/2}-\frac{a(1/2)}{6J^{2}}T^{3}-\frac{b(1/2)}{20J^{4}}T^{5},$
(36) $\displaystyle F_{\mathrm{D}}(T,J)/N_{t}$
$\displaystyle=F_{\mathrm{D}0}(T)/N_{t}-TS_{\mathrm{D}}/N_{t}.$ (37)

Values for $a(1/2)=0.8252$, $a(3/2)=0.1146$, $b(1/2)=7.2$, and $b(3/2)=0.067$
are provided by Ref. [47] and the value
$\varepsilon_{\text{AFM}}^{S=1/2}=-0.669437(5)$ by Ref. [48]. We did not find
a value for $\varepsilon_{\text{AFM}}^{S=3/2}$ in the literature, and so we
extracted the result $\varepsilon_{\text{AFM}}^{S=3/2}=-4.98603(3)$ from our
own SSE QMC simulations, as detailed in App. A.

Solving the equation

$\displaystyle F_{\mathrm{Q}}(T,J_{c})$
$\displaystyle=F_{\mathrm{D}}(T,J_{c})$ (38)

for different, fixed values of $J_{2}/J_{1}$ then yields the transition lines,
$J_{c}(T)$, shown in Fig. 4(e). Direct comparison with the data shown for a
finite-sized system in Figs. 4(a-d) confirms that the approximate treatment
holds remarkably well in describing the first-order transition line. In fact,
within the finite grid of our $(J,T)$ coverage, we do not detect any
deviations from it, even up to $T_{c}$. However, as discussed earlier, the
existence of the critical point is not visible by this approach.

The approximate analysis provides valuable insight into the shapes of the
transition lines in Fig. 4(e). In the limit of low temperatures and when
$J_{2}\neq J_{1}$, $S_{D}$ is suppressed exponentially and $J_{c}(T)$ is
dominated by the $T^{3}$ magnon contribution, leading to the asymptotic
behavior

$J_{c}(T)\sim
J_{c}+\frac{\left(\varepsilon_{\text{AFM}}^{S=3/2}-\varepsilon_{\text{AFM}}^{S=1/2}\right)}{6(\epsilon_{\mathrm{Q}}-\epsilon_{\mathrm{D0}})}\left(a(3/2)-a(1/2)\right)\,T^{3},$
(39)

which describes an asymptotically vertical transition line. However, in the
special case $J_{1}=J_{2}=J_{3}$, $S_{D}=N_{t}\log 2$ never becomes irrelevant
due to the extensive ground-state degeneracy in the doublet phase. Because
$-TS_{D}$ enters linearly, the result is an asymptotically sloped first-order
line,

$J_{c}(T)\sim J_{c}+\frac{\log
2}{\varepsilon_{\text{AFM}}^{S=1/2}-\varepsilon_{\text{AFM}}^{S=3/2}}\,T,$
(40)

accounting for the clear qualitative difference in lineshapes in Fig. 4(e). A
closer inspection of Eq. (34) reveals that obtaining a $T$-linear free-energy
contribution, and thus a first-order line with finite slope at $T=0$, is
possible only in the presence such a macroscopic jump in the ground-state
entropy.

For intermediate $J_{2}/J_{1}$ and at finite temperatures, the doublet entropy
term, $-TS_{D}$, competes with the magnon terms, leading to a smooth crossover
between the regime of well separated doublets and the case where they are
degenerate (Fig. 4(e)). Even close to $T_{c}$, the slope of the first-order
line is influenced by this competition, allowing us to alter the slope at
$T_{c}$ to a certain extent by changing $J_{2}/J_{1}$. We will explore the
effects of this parameter in Subsecs. 4.3 and 4.4. However, first we focus on
locating the critical point, which is not possible within the approximate
analytical approach, and for this we return to QMC simulations and systematic
finite-size scaling.

###  4.2 Thermal critical point

Figure 5: Determination of the critical point of the FFTL model for different
intra-trimer couplings, $J_{2}/J_{1}$, by finite-size scaling of the maximum
in the quartet susceptibility, $\chi_{Q}(T,J)$, computed along several
fixed-$T$ cuts. (left column) $\chi_{Q}(T\approx T_{c},J)$ for different $L$.
The maxima are fitted using a Lorentzian (black). (center column) Scaling of
the maximum value, $\chi^{\text{max}}_{Q}$, normalized by the power-law form,
$L^{\gamma/\nu}=L^{7/4}$, expected for the 2D Ising universality class. (right
column) Scaling of the position, $J^{\text{max}}$, of this maximum.

As we saw in Fig. 4, the first-order transition line, $J_{c}(T)$, continues
only up to a critical temperature, $T_{c}$. Because the fluctuations that
proliferate at this critical point are those of the scalar $l_{\Delta}$ degree
of freedom, one expects by analogy with the phase diagram of water, and of the
FFB [27] and SSL [32], a 2D Ising-type criticality in the FFTL. We now extract
$(T_{c},J_{c}(T_{c}))$ for the different values of $J_{2}/J_{1}$ shown in Fig.
4.

Because both the critical temperature, $T_{c}$, and the critical coupling,
$J_{c}(T_{c})$, are unknown, locating the critical point numerically requires
the adjustment of two parameters. To accomplish this, we calculated the
“quartet susceptibility,”

$\displaystyle\chi_{Q}=\frac{\partial^{2}f}{\partial\varepsilon_{Q}^{2}}$
$\displaystyle=\frac{\beta}{N_{t}}\left(\Braket{N_{Q}^{2}}-\Braket{N_{Q}}^{2}\right),$
(41)
$\displaystyle=\frac{\beta}{N_{t}}\left(\Braket{\left(\sum_{\Delta}l_{\Delta}\right)^{2}}-\Braket{\sum_{\Delta}l_{\Delta}}^{2}\right),$
(42)

along lines of fixed $T$. $\chi_{Q}$ reflects the fluctuations in the number
of $l_{\Delta}=3/2$ sites, $N_{Q}$, or equivalently fluctuations in the sum
over all $l_{\Delta}$, and diverges along the first-order line, $J_{c}(T)$
with $T\leq T_{c}$, whereas for $T>T_{c}$ the divergence turns into a finite
peak. In a finite system, any divergence at a phase transition is regularized
by the system size and also appears as a finite peak. We therefore search for
the critical point along the line of susceptibility maxima,
$(T,J^{\text{max}}(T))$, where the peak height of the susceptibility will
follow a (divergent) power-law form as a function of the system size,

$\chi^{\text{max}}_{Q}(T_{c})\sim L^{\gamma/\nu},$ (43)

with the critical exponents $\gamma$ and $\nu$. If we assume a 2D Ising
critical point, we can find its position by ensuring the best match between
the finite-size scaling of the maximum and the exactly known exponents
$\gamma=7/4$ and $\nu=1$ [49].

In Fig. 5, this program is carried out for different values of $J_{2}/J_{1}$.
In practice, the simulations become less efficient close to the first-order
transition, making the extraction of the maximum, $\chi_{Q}^{\text{max}}$,
intractable for larger system sizes and $T\lesssim T_{c}$. Nevertheless, by
using the known 2D Ising exponents, we are able to extract accurate estimates.
Their error is dominated by the uncertainty in $T_{c}$, which is mostly a
consequence of the maximum system sizes we can simulate. We find that
$J_{c}(T_{c})$ rises significantly, whereas $T_{c}$ itself is enhanced only
weakly, as $J_{2}/J_{1}$ approaches $1$. Marking the positions of these
critical points in the phase diagrams of Figs. 4(a-d), we find again
remarkable agreement between the QMC data and the analytical estimate for the
first-order line: the fact that the critical points fall on top of the
estimated lines shows that the approximations of Subsec. 4.1 perform well even
up to criticality. The rise of $J_{c}(T_{c})$ with $J_{2}/J_{1}$ reflects the
increasing slope of the first-order lines with $J_{2}/J_{1}$, whose
thermodynamic consequences we consider next.

###  4.3 Specific heat

The critical point of the FFTL, and of similar models such as the FFB and SSL
[27, 32], can be viewed up to the differing dimensionality ($d=2$ vs. $d=3$)
as a relative to the critical point that terminates the liquid-gas transition
line of water. Indeed, these spin models share the critical exponents of the
$d=2$ Ising universality class, in the same way that the critical point of
water shows $d=3$ Ising universality. Further, the critical point in these
quantum spin models comes about without the presence of an explicit $Z_{2}$
symmetry, making them even closer relatives of water than is the $d=3$ Ising
model itself.

At a second glance, however, this correspondence is not trivial. An immediate
example arises in a comparison of the specific heat that was performed
recently [32]. In the quantum spin systems, the specific-heat data contain two
lines of maxima that both approach the critical point sideways relative to the
first-order line. While one line of maxima is well known in the specific heat
of water, appearing as an extension of the first-order line beyond the
critical point and into the supercritical regime, there is no equivalent to
the second line of maxima seen in the spin systems. Another distinction is the
slope of the first-order line at the critical point with respect to the
temperature axis, the FFB and SSL models featuring rather vertical lines quite
opposed to the strongly slanted one in water (explicit visualizations may be
found in Ref. [32]). These differences raise the general question of the
extent to which the two lines of maxima can be a universal feature of the
critical point, rather than a model-specific property, and we would like in
particular to understand their connection with the slope of the first-order
line.

Figure 6: Specific heat, $C(J,T)$, of the $L=12$ FFTL model for different
values of $J_{2}/J_{1}$. (a)–(d) $C$ as a function of $J$ and $T$. For each
fixed $J$, the temperatures, $T^{\text{max}}$, of the maxima in $C(T)$ are
shown as white dots. Red lines and points mark again the approximate first-
order line and the numerically determined critical point. Arrows indicate the
positions of the fixed-$J$ cuts shown in panels (e)–(h). In panel (h), the
inset shows a magnification of the cut at $J/J_{1}=0.26$, highlighting the
two-maxima structure.

The FFTL, in which $J_{2}/J_{1}$ produces quite explicitly the differently
slanted first-order transition lines shown in Fig. 4(e), thus provides a
valuable model for elucidating this connection. In Fig. 6 we show the specific
heat, $C=\partial\langle H\rangle/\partial T$, in the $T$-$J$ plane for 4
different values of $J_{2}$. In the vicinity of the critical point, we indeed
identify two lines of maxima, whose shapes depend on the strength of $J_{2}$.
Although neither the lines nor the actual $C(T)$ curves are ever symmetrical,
for $J_{2}=0$ (Figs. 6(a,e)) the lines of maxima do appear at comparable
temperatures. This is similar to the specific-heat results for the FFB and SSL
models [32], in which we stress that there is no exact symmetry, whereas this
symmetry is an exact property of the Ising model (enforced by the $Z_{2}$
symmetry of the Hamiltonian and displayed in Fig. 7(a)).

As $J_{2}$ is increased and the first-order line slants further to the right,
both the temperatures (Figs. 6(b-d)) and peak heights (Figs. 6(f-h)) of the
specific heat are suppressed on the left side of the phase diagram ($S=1/2$
AFM). At $J_{2}/J_{1}=1$ (Fig. 6(d,h)), the maxima to the left of $J_{c}$ are
barely discernible, least of all when compared to the robust peak in their
right-side counterparts, and the situation resembles the specific heat of
water, where only one line of maxima is visible. In this way, the FFTL
provides a continuous progression from the “Ising-like” to the “water-like”
form of the specific heat within a single model, making it a useful example to
understand their differences. At the microscopic level, it is the doublet-
degeneracy of the FFTL model that allows us to control this progression, and
to proceed all the way to the asymmetric, “water-like” form in a manner that
was not possible in the FFB and SSL models (which are not far from “Ising-
like” in this respect).

Focusing in detail on the positions of the maxima, we observe that only the
right-hand line of maxima approaches the critical point directly. The left-
hand line of maxima instead meets the first-order line below the critical
temperature, $T_{c}$. This behavior is most obvious around $J_{2}/J_{1}=1$ and
becomes harder to distinguish for smaller $J_{2}$. However, its presence does
point to the fact that the left-hand maxima are not a universal feature of the
critical point, and thus it is not surprising that there are models where only
one line of maxima appears.

Returning to the analysis of Subsec. 4.1, the specific heat of the FFTL for
$J<J_{c}$ can be understood as the sum of an $S=1/2$ AFM magnon contribution
and the specific heat of decoupled trimers ($J=0$). The magnon contribution
consists of a peak at the temperature scale of $J$, which is independent of
the value of $J_{2}$. Conversely, the trimer contribution is independent of
$J$ but depends strongly on $J_{2}$. When the degeneracy between
$\epsilon_{\mathrm{D}0}$ and $\epsilon_{\mathrm{D}1}$ is broken, the
$\mathrm{D}1$ level moves upwards into the gap between
$\epsilon_{\mathrm{D}0}$ and $\epsilon_{\mathrm{Q}}$, leading for small
$J_{2}/J_{1}$ to a “doublet peak” that obscures the magnon peak. As
$J_{2}/J_{1}\rightarrow 1$, the gap between the doublet levels closes and the
associated specific-heat peak moves to lower temperatures, revealing more of
the magnon peak. In the degenerate case ($J_{2}/J_{1}=1$), the doublet entropy
is not released at all, leading to a very weak specific-heat signal. In this
situation, the magnon and trimer contributions are sufficiently well separated
at weak inter-trimer couplings ($J/J_{1}\lesssim 0.27$) that they produce two
clearly distinguishable, if very weak, local maxima (inset, Fig. 6(h)).

From another perspective, one may ask how much of the specific heat can be
understood by considering the free energy close to the critical point, where
its singular part will become similar to the part associated with the 2D Ising
universality class. This similarity, however, holds only up to a coordinate
transformation in the scaling fields, i.e. the quantities known as the
temperature and magnetic field in the bare Ising model need not correspond to
the physical temperature and couplings of the FFTL.

This type of coordinate transformation not only controls the angle and
location of the first-order line in the $(J,T)$ phase diagram, but also has a
direct effect on the thermodynamic quantities by transforming and mixing their
free-energy derivatives. To elucidate this point, we performed such a
transformation explicitly for the specific heat of the bare Ising model on the
square lattice,

$H=-\sum_{\braket{i,j}}\sigma^{z}_{i}\sigma^{z}_{j}-\tilde{h}\sum_{i}\sigma^{z}_{i},$
(44)

at “bare” temperature $\tilde{T}$ and magnetic field $\tilde{h}$. To compare
with the FFTL, which has a variable slanting of the first-order line, we
introduce a coordinate transformation that is essentially a rotation around
the critical point,

$\begin{pmatrix}T-T_{c}\\\
h\end{pmatrix}=\begin{pmatrix}\cos(\phi)&-\sin(\phi)\\\
\sin(\phi)&\phantom{-}\cos(\phi)\end{pmatrix}\begin{pmatrix}\tilde{T}-T_{c}\\\
\tilde{h}\end{pmatrix};$ (45)

(a related transformation was considered in the context of liquid-liquid
transitions in Ref. [50]). Under this transformation, the specific heat mixes
with the magnetocaloric correlations and the magnetic susceptibility to give

$-C/T=\frac{\partial^{2}f}{\partial
T^{2}}=\cos^{2}(\phi)\,\frac{\partial^{2}f}{\partial\tilde{T}^{2}}+\cos(\phi)\sin(\phi)\,\frac{\partial^{2}f}{\partial\tilde{T}\partial\tilde{h}}+\sin^{2}(\phi)\,\frac{\partial^{2}f}{\partial\tilde{h}^{2}}.$
(46)

Figure 7: Specific heat, $C(T,h)$, of the bare $L=12$ Ising model in rotated
coordinates, shown for different values of the mixing angle, $\phi$, between
the bare coordinates $\tilde{T}$ and $\tilde{h}$.

We evaluate these three free-energy derivatives for the Ising model by
computing the correlation functions $\braket{H^{2}}-\braket{H}^{2}$,
$\braket{HM}-\braket{H}\\!\braket{M}$, and $\braket{M^{2}}-\braket{M}^{2}$,
where $M=\sum_{i}\sigma^{z}$ is the Ising magnetization. For this purpose we
employed classical Monte Carlo simulations with Metropolis [51] and Wolff [52]
updates using the “ghost spin” method [53] to simulate efficiently in a
magnetic field. Figure 7 shows the “rotated” specific-heat results we obtained
for different values of the coordinate mixing angle, $\phi$. Looking again at
the maxima, we observe that singular contributions appear along the first-
order line as soon as $\phi>0$, and that the right-hand maxima above it are
enhanced. The more slanted the first-order line becomes, the more strongly is
the right branch enhanced, until even at comparatively small $\phi$ values it
dominates the signal completely.

This picture is in good agreement with the specific heat of the FFTL close to
the critical point in Fig. 6, where all cases feature slanted first-order
lines whose angle increases with their proximity to the degenerate case.
Similar behavior is also present in the specific-heat data for the FFB and SSL
models [32]. At this point, we stress that the construction of the “rotated”
free energy has a meaning only in the vicinity of the critical point: it
cannot be written in terms of a partition function of a microscopic model and
away from criticality it is the microscopic degrees of freedom that are the
most important. Specifically for our discussion of the full first-order lines,
their low-temperature shapes are those given by the free-energy arguments of
Subsec. 4.1 and shown in Fig. 4(e).

Figure 8: Finite-size scaling of the specific heat, $C(J,T)$, close to the
critical point for (a) a non-degenerate FFTL and (b) the degenerate FFTL. (c)
Scaling of the specific-heat maxima, $C^{\text{max}}$, shown together with the
quartet susceptibility maxima, $\chi_{Q}^{\text{max}}$, normalized for
comparison to their $L=12$ values. The dashed black line shows a power law
with the critical exponent $\gamma/\nu=7/4$. The temperatures are those of
panels (a) and (b).

The physics underlying these results is the fact that the most divergent
quantity in the critical scaling of the Ising model is the susceptibility,
given by the third term of Eq. (46), whereas the Ising specific heat (first
term) has only a weak critical scaling ($\log L$ in 2D). The critical
properties of any system described by a rotated Ising model are a mixture of
the pure Ising terms, and can therefore be expected to change their form. This
behavior is seen most clearly in the critical specific heat of the FFTL, which
in contrast to the 2D Ising form grows strongly with system size, as we show
in Figs. 8(a-b). In the degenerate case ($J_{2}/J_{1}=1$), the scaling of the
specific heat coincides with that of the quartet susceptibility up to a
normalization factor (Fig. 8(c)). For $J_{2}/J_{1}=0$, the specific heat does
grow more slowly than the quartet susceptiblity, in agreement with the weaker
scaling-field mixing of the non-degenerate case, although the susceptibility
contribution is still expected to dominate at large system sizes. A purely
Ising-type $\log L$ scaling of the specific heat has been observed at the
thermal critical point in the FFB [27], indicating that the mixing of scaling
fields is weaker in that model than in the maximally non-degenerate FFTL. For
the degenerate FFTL, the match in the $L^{\gamma/\nu}$ scaling behavior of the
specific heat and the quartet susceptibility close to critical point can be
interpreted thermodynamically as an effect of scaling-field mixing and
microscopically by the fact that the energy fluctuations (reflected in $C$)
are dominated by quartet fluctuations ($\chi_{Q}$) in the absence of a gap
between the doublets (Fig. 6(h)).

Finally, in principle one may question the effect of this scaling-field mixing
on our estimation of the critical point in Subsec. 4.2, where we compared the
growth of the quartet susceptibility directly to the critical scaling of the
Ising susceptibility. However, these two quantities no longer have the same
direct relationship when the scaling fields become mixed: similar to Eq. (46),
the scaling of the singular part of the quartet susceptiblity is given by a
combination of Ising observables. However, once again the most divergent
contribution is that of the Ising susceptibility, which therefore dominates in
the thermodynamic limit, while the other mixed-in contributions behave as
subleading corrections to scaling and enter our estimate of the critical point
in the same way as additional finite-size effects. The contributions from
weaker divergences can be significant in some quantities for the system sizes
of our study, most notably for $C(T)$ at smaller values of $J_{2}/J_{1}$.
However, the results of Fig. 8(c) make clear that the quartet susceptibility,
$\chi_{Q}$, of the FFTL is affected only weakly by the mixing of scaling
fields for all values of $J_{2}/J_{1}$, making it eminently suitable for a
robust extraction of $T_{c}$ and $J_{c}(T_{c})$. In fact the minor deviation
from $L^{\gamma/\nu}$ behavior of $\chi_{Q}$ is most likely the consequence of
a small discrepancy from the exact critical point, which would allow our
estimate to be further improved.

Figure 9: Quartet correlation length, $\xi_{Q}$, shown as a function of $T$
and $J$ for the $L=12$ FFTL model at (a) a representative non-degenerate
intra-trimer coupling, $J_{2}/J_{1}=0$, and (b) the degenerate point,
$J_{2}/J_{1}=1$. Superimposed are the lines of constant-$J$ maxima in the
specific heat, $T^{\text{max}}$, and the “critical isochore” line where
$\braket{l_{\Delta}}=1$. In both panels, the critical point and the
approximate first-order line are shown in red.

###  4.4 Correlation length and characteristic lines

Additional insight into the physics above $T_{c}$ can be obtained by
considering the correlation length associated with the critical point. This we
estimate [42] from the quartet structure factor,

$S_{Q}(\mathbf{q})=\frac{1}{N_{t}}\sum_{\Delta,\Delta^{\prime}}e^{-i\mathbf{q}\cdot(\mathbf{R}_{\Delta}-\mathbf{R}_{\Delta^{\prime}})}\left(\braket{l_{\Delta}l_{\Delta^{\prime}}}-\braket{l_{\Delta}}\\!\braket{l_{\Delta^{\prime}}}\right),$
(47)

as

$\xi_{Q}=\frac{1}{|\delta\mathbf{q}|}\sqrt{\frac{S_{Q}(\delta\mathbf{q})/S_{Q}(2\delta\mathbf{q})-1}{4-S_{Q}(\delta\mathbf{q})/S_{Q}(2\delta\mathbf{q})}},\quad\delta\mathbf{q}=\frac{2\pi}{L}\hat{\mathbf{x}},$
(48)

where $\mathbf{R}_{\Delta}$ is the position of trimer $\Delta$. Based on this
standard Ornstein-Zernike form, we extract the length scale of the critical
fluctuations in terms of the corresponding correlation function.

Figure 9 shows scans of $\xi_{Q}$ across the phase diagram for the two cases
$J_{2}=0$ and $J_{2}/J_{1}=1$. These plots look rather similar to the
specific-heat data of Fig. 6, with the notable exception that there is no
trace of a second set of maxima to the left of the first-order line. From
$\xi_{Q}$ we therefore identify a single characteristic line within the
supercritical regime, which in essence continues the first-order line beyond
$T_{c}$. Along this continuation line, fluctuations in the local trimer states
proliferate, anticipating the critical point and then the first-order line
that control the low-temperature physics. For comparison, we include in Fig. 9
two other sets of characteristic quantities. One set is the lines of specific-
heat maxima from Fig. 6, although for further comparison we neglect the line
of weak maxima to the left of $J_{c}$. The other is the line along which
$\braket{l_{\Delta}}=1$, a value we base on our estimate that
$\braket{l_{\Delta}}\approx 1$ at the critical point. By analogy with the
density of water at the critical point, we refer to this line as the “critical
isochore.”

In the supercritical regime, the lines of maxima from these three different
thermodynamic quantities converge upon approaching the critical point for any
value of $J_{2}$. Away from criticality, they begin to diverge, each
apparently most sensitive to slightly different features of the non-divergent
thermal fluctuations. It is notable that both the convergence near the
critical point and the proximity of the characteristic lines further from
criticality are best in the degenerate limit, $J_{2}/J_{1}=1$, when the first-
order line is most “water-like.” Analogous observations regarding the loci of
extrema of various response functions were also made in Ref. [50] with a view
to determining the Widom line (strictly, the line of zero ordering field) in
the critical phenomena of liquid-liquid transitions with differently slanted
first-order lines. In this sense, the characteristic lines we have identified
may all serve as quantum magnetic analogs of the Widom line, particularly for
the doublet-degenerate FFTL.

##  5 Conclusion

We have presented a QMC approach to study highly frustrated quantum spin
models in the trimer basis. The abstract directed-loop approach that we
introduce requires no system-specific knowledge or bias. It generalizes easily
to other computational bases, and thus should be of use for QMC simulations of
many other classes of frustrated quantum spin system (for example those
composed of four-spin clusters). Here we have applied our algorithm to the
$S=1/2$ Heisenberg antiferromagnet on the FFTL, where it performs efficient
and sign-problem-free simulations up to large system sizes, allowing us to
investigate the full finite-temperature physics of this model.

A key property of the FFTL model is the appearance of a first-order quantum
phase transition between two AFM ground states that break the same symmetry.
From the QMC simulations, we found that a line of first-order transitions
emerges from this point at finite temperatures, and terminates at a thermal
critical point in the 2D Ising universality class. We used a straightforward
free-energy argument to understand the shape of the first-order line in terms
of the low-lying excitations and the doublet entropy intrinsic to each trimer.
What sets the FFTL apart from the quantum spin models studied to date in this
context is the curvature of the first-order line induced by the doublet
entropy contribution, which develops in the limit of degenerate doublets into
a fully slanting line in the plane of temperature and coupling ratio.

We demonstrated further that this increased slanting of the first-order line
is associated with an evolution of the specific heat from the Ising limit,
which features two symmetrical lines of maxima, to a single prominent line of
maxima that appears as an extension of the first-order line into the
supercritical regime. This effect can be explained within the Ising model by a
coordinate transformation in the vicinity of the critical point. We computed
different characteristic thermodynamic quantities for the FFTL, notably the
line of maxima in the correlation length and the critical isochore, to show
that they align with the line of primary specific-heat maxima upon approaching
the critical point from the supercritical regime, whereas the line of
secondary maxima is not related to the critical-point physics. In this way we
identified the quantum magnetic analog of the Widom line, serving as a
continuation of the first-order line in the supercritical regime.

Our results show that the FFTL model is very valuable for exploring the
physics of first-order (quantum) phase transitions and thermal criticality in
the absence of symmetry-breaking in frustrated quantum magnets. The model
possesses a number of attributes that warrant further investigation. In
particular, the first-order quantum transition separates two AFM ground states
that break the same symmetry. In principle, this allows quantum fluctuations
in an extended phase diagram to terminate the first-order transition at $T=0$,
leading to a quantum critical point. Beyond such a point, one may anticipate a
smooth crossover between the $S=1/2$ and $S=3/2$ representations or an exotic
intermediate phase that mixes the two. It is conceivable that this regime
could be reached within the sign-free parameter subspace, for example by
introducing a finite magnetic field or an anisotropy in the fully frustrated
interactions, both of which would be accessible by the same abstract directed-
loop QMC approach that we have developed here.

## Acknowledgements

We thank Fabien Alet for helpful discussions. We acknowledge the support of
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through
Grant No. WE/3649/4-2 of the program FOR 1807 and through project RTG 1995. We
are grateful to the Swiss National Science Foundation, and to the European
Research Council (ERC) for funding under the European Union Horizon 2020
research and innovation program (Grant No. 677061). We thank the IT Center at
RWTH Aachen University and the JSC Jülich for access to computing time through
JARA-HPC.

##  Appendix A Ground-state energy of the spin-3/2 Heisenberg model

The free-energy analysis of Subsec. 4.1 requires the ground-state energy per
site of the spin-3/2 Heisenberg model on the square lattice. To compute this
energy, we performed SSE QMC simulations of the model for several linear
system sizes up to $L=32$ and for inverse temperatures up to $\beta=2L$. The
AFM ground-state energy in the finite system converges to the thermodynamic
limit with the form [44]

$\varepsilon_{\text{AFM}}^{S=3/2}(L)=\varepsilon_{\text{AFM}}^{S=3/2}+c/\beta{}L^{2}+\dots,$
(49)

where $c$ is a non-universal constant. We use this relation to extract the
value $\varepsilon_{\text{AFM}}^{S=3/2}=-4.98603(3)$ (Fig. 10).

Figure 10: Ground-state energy per site,
$\varepsilon^{S=3/2}_{\mathrm{AFM}}(L)$ of the spin-3/2 Heisenberg model. The
dashed line shows a fit to Eq. (49), which we use to extract the value in the
thermodynamic limit.

## References

  * [1] O. A. Starykh,  _Unusual ordered phases of highly frustrated magnets: a review_ ,  Rep. Prog. Phys. 78, 052502 (2015),  10.1088/0034-4885/78/5/052502.
  * [2] L. Balents,  _Spin liquids in frustrated magnets_ ,  Nature 464, 199 (2010),  10.1038/nature08917.
  * [3] D. S. Rokhsar and S. A. Kivelson,  _Superconductivity and the quantum hard-core dimer gas_ ,  Phys. Rev. Lett. 61, 2376 (1988),  10.1103/PhysRevLett.61.2376.
  * [4] A. Kitaev,  _Fault-tolerant quantum computation by anyons_ ,  Ann. Phys. 303, 2 (2003),  10.1016/S0003-4916(02)00018-0.
  * [5] A. Kitaev,  _Anyons in an exactly solved model and beyond_ ,  Ann. Phys. 321, 2 (2006),  10.1016/j.aop.2005.10.005.
  * [6] M. Hermele, M. P. A. Fisher and L. Balents,  _Pyrochlore photons: The ${U}(1)$ spin liquid in a ${S}=\frac{1}{2}$ three-dimensional frustrated magnet_,  Phys. Rev. B 69, 064404 (2004),  10.1103/PhysRevB.69.064404.
  * [7] F. Heidrich-Meisner, A. Honecker and T. Vekua,  _Frustrated ferromagnetic spin- $\frac{1}{2}$ chain in a magnetic field: The phase diagram and thermodynamic properties_,  Phys. Rev. B 74, 020403 (2006),  10.1103/PhysRevB.74.020403.
  * [8] A. M. Läuchli, J. Sudan and R. Moessner,  _${S}=\frac{1}{2}$ kagome Heisenberg antiferromagnet revisited_,  Phys. Rev. B 100, 155142 (2019),  10.1103/PhysRevB.100.155142.
  * [9] A. Wietek and A. M. Läuchli,  _Valence bond solid and possible deconfined quantum criticality in an extended kagome lattice Heisenberg antiferromagnet_ ,  Phys. Rev. B 102, 020411 (2020),  10.1103/PhysRevB.102.020411.
  * [10] S. Yan, D. A. Huse and S. R. White,  _Spin-liquid ground state of the S = 1/2 kagome Heisenberg antiferromagnet_ ,  Science 332, 1173 (2011),  10.1126/science.1201080.
  * [11] S. Depenbrock, I. P. McCulloch and U. Schollwöck,  _Nature of the spin-liquid ground state of the $S=1/2$ Heisenberg model on the kagome lattice_,  Phys. Rev. Lett. 109, 067201 (2012),  10.1103/PhysRevLett.109.067201.
  * [12] Y.-C. He, M. P. Zaletel, M. Oshikawa and F. Pollmann,  _Signatures of Dirac cones in a DMRG study of the kagome Heisenberg model_ ,  Phys. Rev. X 7, 031020 (2017),  10.1103/PhysRevX.7.031020.
  * [13] L. Wang and A. W. Sandvik,  _Critical level crossings and gapless spin liquid in the square-lattice spin- $1/2$ ${J}_{1}-{J}_{2}$ Heisenberg antiferromagnet_,  Phys. Rev. Lett. 121, 107202 (2018),  10.1103/PhysRevLett.121.107202.
  * [14] P. Corboz and F. Mila,  _Crystals of bound states in the magnetization plateaus of the Shastry-Sutherland model_ ,  Phys. Rev. Lett. 112, 147203 (2014),  10.1103/PhysRevLett.112.147203.
  * [15] H. J. Liao, Z. Y. Xie, J. Chen, Z. Y. Liu, H. D. Xie, R. Z. Huang, B. Normand and T. Xiang,  _Gapless spin-liquid ground state in the ${S}=1/2$ kagome antiferromagnet_,  Phys. Rev. Lett. 118, 137202 (2017),  10.1103/PhysRevLett.118.137202.
  * [16] Z. Shi, W. Steinhardt, D. Graf, P. Corboz, F. Weickert, N. Harrison, M. Jaime, C. Marjerrison, H. A. Dabkowska, F. Mila and S. Haravifard,  _Emergent bound states and impurity pairs in chemically doped Shastry-Sutherland system_ ,  Nature Commun. 10, 2439 (2019),  10.1038/s41467-019-10410-x.
  * [17] A. Kshetrimayum, C. Balz, B. Lake and J. Eisert,  _Tensor network investigation of the double layer kagome compound Ca 10Cr7O28_,  Ann. Phys. 421, 168292 (2020),  10.1016/j.aop.2020.168292.
  * [18] A. W. Sandvik,  _A generalization of Handscomb’s quantum Monte Carlo scheme-application to the 1D Hubbard model_ ,  J. Phys. A 25, 3667 (1992),  10.1088/0305-4470/25/13/017.
  * [19] O. F. Syljuåsen and A. W. Sandvik,  _Quantum Monte Carlo with directed loops_ ,  Phys. Rev. E 66, 046701 (2002),  10.1103/PhysRevE.66.046701.
  * [20] O. F. Syljuåsen,  _Directed loop updates for quantum lattice models_ ,  Phys. Rev. E 67, 046701 (2003),  10.1103/PhysRevE.67.046701.
  * [21] M. Troyer and U.-J. Wiese,  _Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations_ ,  Phys. Rev. Lett. 94, 170201 (2005),  10.1103/PhysRevLett.94.170201.
  * [22] D. Hangleiter, I. Roth, D. Nagaj and J. Eisert,  _Easing the Monte Carlo sign problem_ ,  Sci. Adv. 6, eabb8341 (2020),  10.1126/sciadv.abb8341.
  * [23] A. Honecker, S. Wessel, R. Kerkdyk, T. Pruschke, F. Mila and B. Normand,  _Thermodynamic properties of highly frustrated quantum spin ladders: Influence of many-particle bound states_ ,  Phys. Rev. B 93, 054408 (2016),  10.1103/PhysRevB.93.054408.
  * [24] F. Alet, K. Damle and S. Pujari,  _Sign-problem-free Monte Carlo simulation of certain frustrated quantum magnets_ ,  Phys. Rev. Lett. 117, 197203 (2016),  10.1103/PhysRevLett.117.197203.
  * [25] S. Wessel, B. Normand, F. Mila and A. Honecker,  _Efficient Quantum Monte Carlo simulations of highly frustrated magnets: the frustrated spin-1/2 ladder_ ,  SciPost Phys. 3, 005 (2017),  10.21468/SciPostPhys.3.1.005.
  * [26] K.-K. Ng and M.-F. Yang,  _Field-induced quantum phases in a frustrated spin-dimer model: A sign-problem-free quantum Monte Carlo study_ ,  Phys. Rev. B 95, 064431 (2017),  10.1103/PhysRevB.95.064431.
  * [27] J. Stapmanns, P. Corboz, F. Mila, A. Honecker, B. Normand and S. Wessel,  _Thermal critical points and quantum critical end point in the frustrated bilayer Heisenberg antiferromagnet_ ,  Phys. Rev. Lett. 121, 127201 (2018),  10.1103/PhysRevLett.121.127201.
  * [28] T. Nakamura,  _Vanishing of the negative-sign problem of quantum Monte Carlo simulations in one-dimensional frustrated spin systems_ ,  Phys. Rev. B 57, R3197 (1998),  10.1103/PhysRevB.57.R3197.
  * [29] S. Wessel, I. Niesen, J. Stapmanns, B. Normand, F. Mila, P. Corboz and A. Honecker,  _Thermodynamic properties of the Shastry-Sutherland model from quantum Monte Carlo simulations_ ,  Phys. Rev. B 98, 174432 (2018),  10.1103/PhysRevB.98.174432.
  * [30] J. D’Emidio, S. Wessel and F. Mila,  _Reduction of the sign problem near $T=0$ in quantum Monte Carlo simulations_,  Phys. Rev. B 102, 064420 (2020),  10.1103/PhysRevB.102.064420.
  * [31] E. Müller-Hartmann, R. R. P. Singh, C. Knetter and G. S. Uhrig,  _Exact demonstration of magnetization plateaus and first-order dimer-Néel phase transitions in a modified Shastry-Sutherland model for SrCu 2(BO3)2_,  Phys. Rev. Lett. 84, 1808 (2000),  10.1103/PhysRevLett.84.1808.
  * [32] J. Larrea Jiménez, S. P. G. Crone, E. Fogh, M. E. Zayed, R. Lortz, E. Pomjakushina, K. Conder, A. M. Läuchli, L. Weber, S. Wessel, A. Honecker, B. Normand _et al._ ,  _A quantum magnetic analogue to the critical point of water_ ,  Nature 592, 370 (2021),  10.1038/s41586-021-03411-8.
  * [33] P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ ,  Cambridge University Press,  10.1017/CBO9780511813467 (1995).
  * [34] T. R. Kirkpatrick and D. Belitz,  _Exponent relations at quantum phase transitions with applications to metallic quantum ferromagnets_ ,  Phys. Rev. B 91, 214407 (2015),  10.1103/PhysRevB.91.214407.
  * [35] A. Georges, G. Kotliar, W. Krauth and M. J. Rozenberg,  _Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions_ ,  Rev. Mod. Phys. 68, 13 (1996),  10.1103/RevModPhys.68.13.
  * [36] A. Honecker, F. Mila and M. Troyer,  _Magnetization plateaux and jumps in a class of frustrated ladders: A simple route to a complex behaviour_ ,  Eur. Phys. J. B 15, 227 (2000),  10.1007/s100510051120.
  * [37] H.-Q. Lin and J. Shen,  _Exact ground states and excited states of net spin models_ ,  J. Phys. Soc. Jpn. 69, 878 (2000),  10.1143/JPSJ.69.878.
  * [38] H. Q. Lin, J. L. Shen and H. Y. Shik,  _Exactly soluble quantum spin models on a double layer: The net spin model_ ,  Phys. Rev. B 66, 184402 (2002),  10.1103/PhysRevB.66.184402.
  * [39] J. Richter, O. Derzhko and T. Krokhmalskii,  _Finite-temperature order-disorder phase transition in a frustrated bilayer quantum Heisenberg antiferromagnet in strong magnetic fields_ ,  Phys. Rev. B 74, 144430 (2006),  10.1103/PhysRevB.74.144430.
  * [40] O. Derzhko, J. Richter, A. Honecker and H.-J. Schmidt,  _Universal properties of highly frustrated quantum magnets in strong magnetic fields_ ,  Low Temp. Phys. 33, 745 (2007),  10.1063/1.2780166.
  * [41] O. Derzhko, T. Krokhmalskii and J. Richter,  _Emergent Ising degrees of freedom in frustrated two-leg ladder and bilayer $s=\frac{1}{2}$ Heisenberg antiferromagnets_,  Phys. Rev. B 82, 214412 (2010),  10.1103/PhysRevB.82.214412.
  * [42] A. W. Sandvik,  _Computational studies of quantum spin systems_ ,  AIP Conf. Proc 1297, 135 (2010),  10.1063/1.3518900.
  * [43] F. Alet, S. Wessel and M. Troyer,  _Generalized directed loop method for quantum Monte Carlo simulations_ ,  Phys. Rev. E 71, 036706 (2005),  10.1103/PhysRevE.71.036706.
  * [44] P. Hasenfratz and F. Niedermayer,  _Finite size and temperature effects in the AF Heisenberg model_ ,  Z. Phys. B 92, 91 (1993),  10.1007/bf01309171.
  * [45] Z. Weihong and C. J. Hamer,  _Spin-wave theory and finite-size scaling for the Heisenberg antiferromagnet_ ,  Phys. Rev. B 47, 7961 (1993),  10.1103/PhysRevB.47.7961.
  * [46] C. J. Hamer, Z. Weihong and J. Oitmaa,  _Spin-wave stiffness of the Heisenberg antiferromagnet at zero temperature_ ,  Phys. Rev. B 50, 6877 (1994),  10.1103/PhysRevB.50.6877.
  * [47] D. C. Johnston, R. J. McQueeney, B. Lake, A. Honecker, M. E. Zhitomirsky, R. Nath, Y. Furukawa, V. P. Antropov and Y. Singh,  _Magnetic exchange interactions in BaMn 2As2: A case study of the ${J}_{1}$-${J}_{2}$-${J}_{c}$ Heisenberg model_,  Phys. Rev. B 84, 094445 (2011),  10.1103/PhysRevB.84.094445.
  * [48] A. W. Sandvik,  _Finite-size scaling of the ground-state parameters of the two-dimensional Heisenberg model_ ,  Phys. Rev. B 56, 11678 (1997),  10.1103/PhysRevB.56.11678.
  * [49] M. E. Fisher and R. J. Burford,  _Theory of critical-point scattering and correlations. I. the Ising model_ ,  Phys. Rev. 156, 583 (1967),  10.1103/PhysRev.156.583.
  * [50] J. Luo, L. Xu, E. Lascaris, H. E. Stanley and S. V. Buldyrev,  _Behavior of the Widom line in critical phenomena_ ,  Phys. Rev. Lett. 112, 135701 (2014),  10.1103/PhysRevLett.112.135701.
  * [51] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller,  _Equation of state calculations by fast computing machines_ ,  J. Chem. Phys 21(6), 1087 (1953),  10.1063/1.1699114.
  * [52] U. Wolff,  _Collective Monte Carlo updating for spin systems_ ,  Phys. Rev. Lett. 62, 361 (1989),  10.1103/PhysRevLett.62.361.
  * [53] J. Kent-Dobias and J. P. Sethna,  _Cluster representations and the Wolff algorithm in arbitrary external fields_ ,  Phys. Rev. E 98, 063306 (2018),  10.1103/PhysRevE.98.063306.