# Toward global fits using Higgs STXS data with Lilith

Dang Bao Nhi Nguyena,b Duc Ninh Lec Sabine Kramld Quang Loc Trana,b Van Dung
Lea,b aDepartment of Theoretical Physics, University of Science, Ho Chi Minh
City 70000, Vietnam
bVietnam National University, Ho Chi Minh City 70000, Vietnam
cFaculty of Fundamental Sciences, PHENIKAA University, Hanoi 12116, Vietnam
dLaboratoire de Physique Subatomique et de Cosmologie, Université Grenoble-
Alpes,
CNRS/IN2P3, 53 Avenue des Martyrs, F-38026 Grenoble, France

###### Abstract

In this talk, we present the program Lilith, a python package for constraining
new physics from Higgs measurements. We discuss the usage of signal strength
results in the latest published version of Lilith, which allows for
constraining deviations from SM Higgs couplings through coupling modifiers.
Moreover, we discuss the on-going development to include Higgs STXS data and
SMEFT parametrizations in Lilith with the aim of performing global fits of the
ATLAS and CMS data. As we point out, detailed information on Standard Model
uncertainties and their correlations is important to enable the proper reuse
of the experimental results.

##  1 Introduction

The discovery of the Higgs boson by ATLAS and CMS in 2012 at CERN has opened a
new window to the universe, in particular for the search for physics beyond
the Standard Model (SM). Of particular interest are the couplings of the 125
GeV Higgs boson to heavy fermions (top, bottom quarks) and to massive gauge
bosons ($W$ and $Z$). The current ATLAS and CMS results from Run 2 of the LHC
already provide very stringent tests of these couplings. The data types
concerning us in this contribution are Higgs signal strengths (SS) and
Simplified Template Cross Sections (STXS), which are published in many Higgs-
analysis papers by ATLAS and CMS since the Higgs discovery; see, e.g., Ref.
[1] for a detailed discussion.

The SS are defined as

$\displaystyle\mu_{i}^{j}=\dfrac{(\sigma_{i}\times\mathcal{B}^{j})_{\text{exp}}}{(\sigma_{i}\times\mathcal{B}^{j})_{\text{SM}}},$
(1)

where $\sigma_{i}$ is the inclusive Higgs production cross section of the
channel $i$, which can be $ggF$ ($87\%$), $VBF$ ($7\%$), $VH$ ($4\%$),
$t\bar{t}H$ ($1\%$), $tH$ ($<0.1\%$). These cross sections are usually
measured with $|y_{H}|<2.5$. $\mathcal{B}^{j}$ is the Higgs branching fraction
in the decay mode $j$, which can be $b\bar{b}$, $WW^{*}$, $ZZ^{*}$,
$\gamma\gamma$, etc. The best-fit values for $\mu_{i}^{j}$ together with their
uncertainties and correlations, as reported in ATLAS and CMS publications, can
be used to constrain deviations of the Higgs couplings from their SM
expectations.

The SS data type is simple and has been used efficiently in many experimental
and theoretical papers to check the consistency of the SM with the ATLAS and
CMS data. Its power for constraining new physics is, however, limited for
(essentially) two reasons: 1) the dominant contribution to the inclusive cross
section comes from the low energy region, while new physics is expected to
affect the high energy region, e.g., when $p_{T}^{H}$ is large; and 2) the SS
framework explicitly assumes that the signal acceptances are the same as in
the SM. The latter means that the SS framework is valid only for overall
scalings, but not for changes in kinematic distributions. The STXS scheme was
introduced to overcome these limitations without going all the way to fiducial
differential distributions. Concretely, the STXS split the inclusive phase
space into different regions, called _bins_ , corresponding to different
ranges of $p_{T}^{H}$, jet multiplicity, etc. Similar to the SS results, ATLAS
and CMS have been providing best-fit values together with uncertainties and
correlations for the STXS.

##  2 Signal strength data and reduced couplings in Lilith

Lilith, Light Likelihood Fit for the Higgs, is a python package for
constraining new physics from Higgs measurements [2, 3]. The current stable
release Lilith-2.1 includes a database with the complete ATLAS and CMS SS data
from Run 1 and from Run 2 with 36/fb, plus some old Tevatron results. The core
of the program is the likelihood calculation based on multi-dimensional
Gaussian distributions of variable width to deal with correlated, asymmetric
uncertainties [3]. A global likelihood is built from this by trivially
combining the individual experimental results in the database.

Within Lilith, the SS are expressed in terms of coupling modifiers $C_{X}$
($X=t,b,c,W,Z,....$), which rescale the Higgs couplings to SM particles with
respect to their SM values [2]. Likelihood fits can thus be performed in terms
of SS or in terms of the reduced couplings $C_{X}$. An example for the latter
is shown in Fig. 1 (left), which compares the Lilith fit of universal scaling
factors for the couplings to fermions, $C_{F}$, and the coupling to massive
gauge bosons, $C_{V}$, to the corresponding results from the CMS measurement
in the $H\to WW^{*}\to 2\ell 2\nu$ decay channel. We note that such
comparisons are particularly useful for validations of the database and of our
statistical method.

This approach can be used for constraining the parameters of any new physics
model that has the same tensor structure as the SM, such that the couplings of
the SM-like Higgs simply rescale. The simplest case are two Higgs doublet
models (2HDMs), where, for decoupled $H^{\pm}$, the Higgs couplings depend on
only on the two angles $\alpha$ and $\beta$. As an example, Fig. 1 (right)
shows constraints in the $\tan\beta$ vs. $\cos(\beta-\alpha)$ plane for the
2HDM of Type I, where $C_{V}=\cos(\beta-\alpha)$ and
$C_{F}=\cos\alpha/\sin\beta$.

Figure 1: Comparisons of Lilith results (colored regions) using SS data
against ATLAS and CMS results (dotted contours). On the left, the red and
orange regions are $68\%$ and $95\%$ CL regions, respectively, from Lilith,
while the dotted contours are the corresponding results from CMS-HIG-20-013.
On the right, the light yellow, orange and red regions correspond to $68\%$,
$95\%$ and $99.7\%$ CL, respectively, from Lilith, while the dotted contour is
the $95\%$ CL result from ATLAS-HIGG-2018-57.

##  3 Extension of Lilith to STXS data and SMEFT fits

We aim to go beyond the SS database and include important new physics
parametrizations. The step taken in this work is to implement the STXS data
type and SM effective field theory (SMEFT) parametrizations in Lilith. A first
comparison between the SS and STXS data types is shown in Fig. 2, still using
the reduced coupling scheme. We can see clearly that the STXS data allow to
better reproduce the official ATLAS result.

Figure 2: Comparison of a fit of $C_{F}$ vs. $C_{V}$ using SS (left) and STXS
(right) data from the ATLAS combined measurements of Higgs boson production
and decay with up to 80/fb, ATLAS-HIGG-2018-57. The red, orange and yellow
regions are the $68\%$, $95\%$ and $99.7\%$ CL regions from Lilith, while the
blue dotted contours are the $68\%$ and $95\%$ CL results from ATLAS.

While the STXS result in Fig. 2 is encouraging, there is still a sizeable
discrepancy between Lilith and ATLAS. This discrepancy is even more pronounced
for the ATLAS measurement in the $H\to\gamma\gamma$ channel with 139/fb,
ATLAS-HIGG-2020-16, shown in Fig. 3. We think that the discrepancy comes to
large extent from the treatment of the theoretical uncertainties, and in
particular from unknown theory correlations. The problem here is that the
determination of the theory correlations is not unique, as it depends on the
uncertainty sources and assumptions on their correlations (see Section I.4.2.a
of Ref. [1] for a nice discussion on this issue). Indeed there exist various
approximations of theory correlations without a guidance which one is the best
choice. For the $ggF$ bins, these approximations are called “2017-scheme”,
“JVE-scheme”, “STXS-scheme”, “WG1-scheme”, provided in Ref. [4] (see also the
LHC HXSWG webpage [5]). For the ATLAS-HIGG-2020-16 data, we observe that the
above approximations give almost identical results, see Fig. 3, without
completely resolving the discrepancy with the ATLAS fit.

Figure 3: Comparisons between different treatments of theoretical
uncertainties, for the STXS data of ATLAS-HIGG-2020-16. Left: without SM
errors, middle: SM errors included assuming the “2017-scheme”, right: SM
errors included assuming the “WG1-scheme”. The results for the “JVE” and
“STXS” schemes are not shown but almost identical to those of the middle and
right plots.

The way forward, in our opinion, would be that the SM uncertainties _and_ at
least one approximation for their correlations are provided together with the
STXS measurements. As of now, ATLAS and CMS have published the best-fit values
together with measurement uncertainties and their correlations. Sometimes, the
SM predictions and their uncertainties are also provided, but information on
the assumed shapes and correlations are missing. Without this information, the
usage of the STXS measurements to constrain new physics is severely limited.

Another ongoing development in Lilith is the extension to SMEFT fits. To this
end, the parametrization of the STXS bins in terms of Wilson coefficients of
the SMEFT is necessary. This is available from the ATLAS measurement in the
$H\to ZZ^{*}\to 4\ell$ channel with 139/fb, ATLAS-HIGG-2018-28. In Fig. 4, we
show the SMEFT fit results obtained with Lilith using the data from ATLAS-
HIGG-2018-28; as can be seen, they agree very well with the official ATLAS
results.

Figure 4: Comparison between Lilith (using the “2017-scheme” $ggF$ SM
correlations) and ATLAS for a SMEFT parametrization. Details of this
parametrization and of the Wilson coefficients $cHB$, $cHW$, $cHG$ are
provided in the ATLAS-HIGG-2018-28 paper.

##  4 Conclusions

We have discussed the usage of Higgs SS and STXS data for constraining new
physics from Higgs measurements, as well as the program Lilith, which has been
developed for this purpose. The current stable release of Lilith makes use of
SS data, and likelihood fits are based on rescaling of the Higgs couplings
with respect to the SM. Ongoing developments include the extension of Lilith
with STXS data and likelihood fits within the SMEFT, for both of which first
results have been presented here. We have also shown that more information on
SM uncertainties and their correlations would be important for the proper
reuse of the experimental results.

## Acknowledgement

This research is funded by the Vietnam National Foundation for Science and
Technology Development (NAFOSTED) under grant number 103.01-2020.17. SK
acknowledges funding by the IN2P3 theory master project BSMGA and the ANR
project EFTatLHC, ANR-22-CE31-0022-02.

## References

## References

  * [1] D. de Florian et al. [LHC Higgs Cross Section Working Group], _Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector_ , [1610.07922].
  * [2] J. Bernon and B. Dumont, _Lilith: a tool for constraining new physics from Higgs measurements_ , _Eur. Phys. J._ C75 (2015) 440 [1502.04138].
  * [3] S. Kraml, T. Q. Loc, D. T. Nhung and L. D. Ninh, _Constraining new physics from Higgs measurements with Lilith: update to LHC Run 2 results_ , _SciPost Phys._ 7 (2019) no.4, 052 [1908.03952].
  * [4] Dag Gillberg, https://dgillber.web.cern.ch/dgillber/ggF_uncertainty_2017/
  * [5] https://twiki.cern.ch/twiki/bin/view/LHCPhysics/LHCHWGFiducialAndSTXS


# Phases of the spin-1/2 Heisenberg antiferromagnet on the diamond-decorated
square lattice in a magnetic field

Nils Caci  Institute for Theoretical Solid State Physics, JARA FIT, and JARA
CSD, RWTH Aachen University, 52056 Aachen, Germany  Katarína Karl’ová
Department of Theoretical Physics and Astrophysics, Faculty of Science, P. J.
Šafárik University, Park Angelinum 9, 04001 Košice, Slovakia  Taras Verkholyak
Institute for Condensed Matter Physics, National Academy of Sciences of
Ukraine, Svientsitskii Street 1, 790 11, L’viv, Ukraine  Jozef Strečka
Department of Theoretical Physics and Astrophysics, Faculty of Science, P. J.
Šafárik University, Park Angelinum 9, 04001 Košice, Slovakia  Stefan Wessel
Institute for Theoretical Solid State Physics, JARA FIT, and JARA CSD, RWTH
Aachen University, 52056 Aachen, Germany  Andreas Honecker  Laboratoire de
Physique Théorique et Modélisation, CNRS UMR 8089,
CY Cergy Paris Université, 95000 Cergy-Pontoise, France

###### Abstract

The spin-1/2 Heisenberg antiferromagnet on the frustrated diamond-decorated
square lattice is known to feature various zero-field ground-state phases,
consisting of extended monomer-dimer and dimer-tetramer ground states as well
as a ferrimagnetic regime. Using a combination of analytical arguments,
density matrix renormalization group (DMRG), exact diagonalization, as well as
sign-problem-free quantum Monte Carlo (QMC) calculations, we investigate the
properties of this system and the related Lieb lattice in the presence of a
finite magnetic field, addressing both the ground-state phase diagram as well
as several thermodynamic properties. In addition to the zero-field ground
states, we find at high magnetic field a spin-canted phase with a continuously
rising magnetization for increasing magnetic field strength, as well as the
fully polarized paramagnetic phase. At intermediate field strength, we
identify a first-order quantum phase transition line between the ferrimagnetic
and the monomer-dimer regime. This first-order line extends to finite
temperatures, terminating in a line of critical points that belong to the
universality class of the two-dimensional Ising model.

##  I Introduction

The study of strongly frustrated quantum magnets is a central topic in
contemporary condensed matter research. Indeed, magnetic frustration,
introduced, e.g., by competing antiferromagnetic exchange couplings, can lead
to the stabilization of non-classical ground states in quantum magnets Richter
_et al._ (2004); Balents (2010); Lacroix  _et al._ (2011); Diep (2013). In
most cases, these non-magnetic states are characterized by the formation of
strong local singlets among small sub-clusters of spins, as well as the
emergence of an extensive ground-state entropy. In the most favorable case, it
is possible to obtain exact analytical expressions for the ground-state
properties, such as for the Shastry-Sutherland model in the regime of strong
dimer coupling Sriram Shastry and Sutherland (1981); Albrecht and Mila (1996);
Miyahara and Ueda (1999, 2003). In this system, quantum spin degrees of
freedom are arranged on a two-dimensional lattice in an orthogonal manner to
form a frustrated array of coupled spin dimers. For strong intra-dimer
coupling (as compared to the inter-dimer coupling), an exact product state of
dimer singlets forms the system’s ground state. Later, it was furthermore
found that the spin-1/2 version of this quantum spin model finds an almost
perfect realization in the copper-based compound SrCu2(BO3)2 Kageyama  _et
al._ (1999); Miyahara and Ueda (1999). This system has since then been studied
extensively with respect to both the ground state and thermal properties
Kageyama  _et al._ (1999, 2000a, 2000b); Lemmens  _et al._ (2000); Gaulin  _et
al._ (2004); Zayed  _et al._ (2014, 2017); Guo  _et al._ (2020); Larrea
Jiménez  _et al._ (2021) as well as its rich physics in the additional
presence of a magnetic field, notably various plateaux in its magnetization
curve Kageyama  _et al._ (1999); Onizuka  _et al._ (2000); Kodama  _et al._
(2002); Sebastian  _et al._ (2008); Takigawa and Mila (2011); Takigawa  _et
al._ (2013); Jaime  _et al._ (2012); Matsuda  _et al._ (2013); Haravifard  _et
al._ (2016).

Figure 1: Illustration of the diamond-decorated square lattice, with a unit
cell indicated (dashed square), along with the labeling of the five different
sites (circles) within the unit cell and the two different exchange couplings
$J_{1}$ (thin black lines) and $J_{2}$ (thick red lines).

The dimerized nature of the low-energy states in the Shastry-Sutherland model
not only gives rise to interesting physics, but is actually also favorable for
a numerical treatment. Indeed, the Shastry-Sutherland model is not only a
showcase for tensor-network approaches Takigawa  _et al._ (2013); Corboz and
Mila (2013, 2014); Wessel  _et al._ (2018); Wietek  _et al._ (2019); Larrea
Jiménez  _et al._ (2021), but it also allows one to use efficient Quantum
Monte Carlo (QMC) simulations throughout a large part of the dimer phase
Wessel  _et al._ (2018); Wietek  _et al._ (2019); Honecker  _et al._ (2022).
Remarkably, the latter extends to a generalized version of the Shastry-
Sutherland model Weihong  _et al._ (1999); Müller-Hartmann  _et al._ (2000)
where in a certain limit, that is equivalent to a fully frustrated bilayer
model Hida (1992); Sandvik and Scalapino (1994); Sommer  _et al._ (2001); Wang
_et al._ (2006), the QMC sign problem disappears completely. The fully
frustrated bilayer model thus becomes accessible to detailed investigations at
finite temperature via QMC simulations Alet  _et al._ (2016); Ng and Yang
(2017); Stapmanns  _et al._ (2018). In fact, the identification of a first-
order line that terminates at a finite-temperature critical point Stapmanns
_et al._ (2018) in the fully frustrated bilayer model was an important guiding
element to identify similar physics in the Shastry-Sutherland model and
ultimately SrCu2(BO3)2 Larrea Jiménez  _et al._ (2021). Note, furthermore,
that the low-energy high-field region of the fully frustrated bilayer model
permits a mapping to a classical lattice gas, thus allowing for a rigorous
treatment of its low-energy thermodynamics, including a finite-temperature
ordering transition Richter  _et al._ (2006); Derzhko  _et al._ (2007).

Figure 2: Ground-state phase diagram of the spin-1/2 Heisenberg
antiferromagnet on the diamond-decorated square lattice in zero magnetic field
as obtained from Ref. Hirose  _et al._ (2017) containing the Lieb-Mattis (LM),
dimer-tetramer (DT) and the monomer-dimer (MD) phase. In the illustration of
the different ground states, blue (orange) ovals denote spin triplet (singlet)
states on the dimers. A tetramer singlet of the DT phase is illustrated by a
rhombus.

Another highly frustrated two-dimensional quantum spin system of coupled
orthogonal spin dimers is the Heisenberg antiferromagnet on the diamond-
decorated square lattice, shown in Fig. 1. This model contains, in addition to
the dimers (along the $J_{2}$ bonds), a further set of spins that are coupled
to other (dimer) spins only by the $J_{1}$ bonds. In the large $J_{2}$-limit,
for $J_{1}/J_{2}\rightarrow 0$, the spins coupled solely through the $J_{1}$
bonds thus lack a partner spin to form a singlet, and we therefore refer to
these spins as monomer spins. Hirose et al. have performed a detailed
investigation of its zero-field ground-state properties Hirose  _et al._
(2016, 2017, 2018, 2020), but little is known otherwise about this model. The
zero temperature zero-field phase diagram exhibits three distinct ground
states, as illustrated in Fig. 2, and promises interesting physics also in
finite fields, respectively at finite temperature. Here, we shortly introduce
these phases, with further details provided in the following sections. In the
case of large dimer coupling $J_{2}$, the ground state is an exact product
state formed by dimer singlet states on all the $J_{2}$ dimers, while the
remaining spins (referred to as monomer spins) are effectively decoupled. This
leads to an extensive ground-state entropy of $\ln(2)$ per unit cell in this
regime ($J_{2}/J_{1}>2$), which is denoted the monomer-dimer (MD) phase. On
the other hand, for weak $J_{2}$, the system prefers to form dimer triplet
states on all the $J_{2}$ dimers, while the monomer spins predominantly orient
themselves opposite to the polarization of the dimers. This leads to a
ferrimagnetic state, akin to the ferrimagnetic ground state of the mixed
spin-1 and spin-1/2 model on the Lieb lattice 111This lattice appeared in the
context of the CuO2 planes of the high-$T_{c}$ superconductors Mattis (1988)
and thanks to extensions in Ref. Lieb (1989) is now known as “Lieb” lattice.
Actually, an Ising version of this model had appeared in the literature
already significantly earlier Fisher (1960a, b), but we will stick to
established terminology. . Its ferrimagnetic polarization follows from the
Lieb-Mattis theorem Lieb and Mattis (1962) in terms of the two different
sublattices of the Lieb lattice. This phase is therefore also denoted by “LM”
in the following. These two phases, MD and LM, are separated by a further
gapped phase, the dimer-tetramer (DT) phase, cf. Fig. 2. In this phase, two
different kinds of local singlets form: besides the $J_{2}$-dimer singlets,
also singlets on larger clusters with four spins are formed: namely, among the
tetramers that are each composed of one $J_{2}$-dimer and its two neighboring
$J_{1}$-coupled monomer spins. In the DT phase, the ground-state manifold is
again highly degenerate and consists of all configurations of closed packings
of tetramers, with the remaining $J_{2}$-dimers forming dimer (two-site)
singlets.

We here examine the spin-1/2 Heisenberg antiferromagnet on the diamond-
decorated square lattice in a magnetic field. In particular, we explore the
ground-state phase diagram in the presence of a finite magnetic field as well
as the thermal properties. For this purpose, we use a combination of
analytical approaches and various computational methods, including exact
numerical diagonalization, density matrix renormalization group (DMRG)
calculations White (1992); White and Huse (1993); Schollwöck (2005) and
stochastic series expansion (SSE) QMC simulations Sandvik and Kurkijärvi
(1991); Sandvik (1999); Syljuåsen and Sandvik (2002), based on a dimer-
decoupling of the Hamiltonian Honecker  _et al._ (2016); Alet  _et al._
(2016), in order to render the QMC sign-problem free.

After introducing the model in more detail in the following section II, we
describe the analytical and computational approaches that we used in Secs. III
and IV, respectively. Our results for the ground-state properties are
presented in Sec. V, and those at finite temperatures in Sec. VI. In passing,
we provide reference data for the mixed spin-1/2 and spin-1 Heisenberg model
on the Lieb lattice, compare also App. A. Finally, we provide our conclusions
and future perspectives in Sec. VII.

##  II Model

In the following, we consider the spin-1/2 Heisenberg antiferromagnet on the
diamond-decorated square lattice in a magnetic field. The lattice is shown
schematically in Fig. 1 and the Hamiltonian of the model is given by

$\displaystyle H$ $\displaystyle=$ $\displaystyle
J_{1}\sum\limits_{i=1}^{N}\Bigl{[}\mathbf{S}_{i,1}\cdot\Bigl{(}\mathbf{S}_{i,2}+\mathbf{S}_{i,3}+\mathbf{S}_{i,4}+\mathbf{S}_{i,5}$
(1)
$\displaystyle+\mathbf{S}_{i-\hat{x},2}+\mathbf{S}_{i-\hat{x},3}+\mathbf{S}_{i-\hat{y},4}+\mathbf{S}_{i-\hat{y},5}\Bigr{)}\Bigr{]}$
$\displaystyle+J_{2}\sum\limits_{i=1}^{N}\Bigl{(}\mathbf{S}_{i,2}\cdot\mathbf{S}_{i,3}+\mathbf{S}_{i,4}\cdot\mathbf{S}_{i,5}\Bigr{)}$
$\displaystyle-h\sum\limits_{i=1}^{N}\sum\limits_{\mu=1}^{5}S_{i,\mu}^{z}\,,$

where $\mathbf{S}_{i,\mu}=({S}_{i,\mu}^{x},{S}_{i,\mu}^{y},{S}_{i,\mu}^{z})$
represents the spin-1/2 operators assigned to the $\mu$-th spin within the
$i$-th unit cell. We denote the corresponding lattice site by $(i,\mu)$.
Furthermore, the index $i-\hat{x}$ ($i-\hat{y}$) refers to the unit cell to
the left (below) the $i$-th unit cell. Here, we consider a finite lattice with
$N$ unit cells and $N_{s}=5N$ sites, imposing periodic boundary conditions,
and with $N\rightarrow\infty$ in the thermodynamic limit (TDL). Typically, we
use square lattices with $N=L^{2}$. Furthermore, $J_{1}$ and $J_{2}$ are the
two exchange interactions drawn in Fig. 1 by black and red lines,
respectively. The last term in $H$ accounts for the Zeeman coupling of the
spin-1/2 particles to an external magnetic field $h$.

The Hamiltonian (1) can also be expressed in terms of the composite spins on
the $2N$ dimers formed by the $J_{2}$ bonds: In each unit cell $i$, a vertical
dimer is formed by the spins $\mathbf{S}_{i,2}$ and $\mathbf{S}_{i,3}$, and
the total dimer spin for this dimer $d$ is then
$\mathbf{S}_{d}=\mathbf{S}_{i,2}+\mathbf{S}_{i,3}$. Similarly, the spins
$\mathbf{S}_{i,4}$ and $\mathbf{S}_{i,5}$ form a horizontal dimer, and in this
case $\mathbf{S}_{d}=\mathbf{S}_{i,4}+\mathbf{S}_{i,5}$. All total dimer spins
represent locally conserved quantities with well defined quantum spin numbers.
The remaining spins $\mathbf{S}_{i,1}$ are referred to as monomer spins. One
can then express the Hamiltonian (1) in a more compact form:

$\displaystyle H\\!$ $\displaystyle=$
$\displaystyle\\!J_{1}\sum_{d=1}^{2N}\sum_{(i,1)\in\mathcal{N}_{d}}\mathbf{S}_{d}\cdot\mathbf{S}_{i,1}+\frac{J_{2}}{2}\sum_{d=1}^{2N}\left(\mathbf{S}_{d}^{2}-\frac{3}{2}\right)$
(2) $\displaystyle-$ $\displaystyle
h\sum_{d=1}^{2N}S_{d}^{z}-h\sum_{i=1}^{N}S_{i,1}^{z},$

where summations over $d$ extend over all the $2N$ dimers, and the inner sum
of the first term extends over the two monomer spins $\mathbf{S}_{i,1}$ that
are nearest neighbors of the $d$-th dimer (cf. Fig. 1), i.e., the lattice site
$(i,1)$ is an element of the set of the two nearest-neighbor sites
$\mathcal{N}_{d}$ of the $d$-th dimer. More specifically, for a vertical
(horizontal) dimer, these are the two monomer spins to the left and right (top
and bottom) of that dimer.

The first term in the Hamiltonian (2) corresponds to the mixed spin-$S_{d}$
and spin-1/2 Heisenberg model on a Lieb lattice, whereas the second term
provides a trivial shift of the energy depending on the quantum spin numbers
$S_{d}$. Note that two different values of the quantum spin numbers $S_{d}=0$,
$1$ are available for the composite spin on each dimer, whereby the value
$S_{d}=0$ corresponds to a singlet-dimer state

$\displaystyle|s\rangle_{d}=\frac{1}{\sqrt{2}}\left(|\\!\uparrow\downarrow\rangle_{d}-|\\!\downarrow\uparrow\rangle_{d}\right).$
(3)

This leads to a fragmentation of the effective mixed-spin Heisenberg models
obtained from the Hamiltonian (2) upon considering all possible combinations
of quantum spin numbers $S_{d}$ for all the dimers. Hence, the ground state of
the Heisenberg antiferromagnet on the diamond-decorated square lattice can be
related to the lowest-energy eigenstates of the effective Heisenberg models
(2) taking into consideration all available combinations of the quantum spin
numbers $S_{d}$. In the following, we first introduce our methods and then
explore the rich ground-state phase diagram of the Hamiltonian $H$, shown
further below in Fig. 3.

##  III Exact analytical ground states

We first consider the parameter regime with a dominant dimer coupling $J_{2}$,
in which we can obtain exact analytical results for the ground state. More
specifically, for $J_{2}/J_{1}>2$ one can use the variational principle in
order to derive an exact ground state of $H$. The main idea of this approach
consists in decomposing the Hamiltonian into $4N$ commuting cell Hamiltonians:

$\displaystyle H=\sum_{d=1}^{2N}\sum_{(i,1)\in\mathcal{N}_{d}}H_{d,i},$ (4)

with each cell Hamiltonian $H_{d,i}$ corresponding to a single triangle
involving one dimer $d$ and one of its two nearest-neighbor monomer spins,
i.e.,

$\displaystyle
H_{d,i}=\frac{J_{2}}{4}\left(\mathbf{S}_{d}^{2}-\frac{3}{2}\right)+J_{1}\mathbf{S}_{i,1}\cdot\mathbf{S}_{d}-\frac{h}{2}S_{d}^{z}-hS_{i,1}^{z}.$
(5)

Note that each dimer $d$ is part of two triangles, leading to the additional
factor of $1/2$ for the intra-dimer term proportional to $J_{2}$ in $H_{d,i}$
as compared to Eq. (2).

According to the variational principle Sriram Shastry and Sutherland (1981);
Bose (1989, 1990, 1992), the ground-state energy of $H$ has a lower bound,
given by the sum of the lowest-energy eigenvalues $\varepsilon^{(0)}_{d,i}$ of
the cell Hamiltonians (5),

$\displaystyle E_{0}$ $\displaystyle=$
$\displaystyle\langle\Psi_{0}|H|\Psi_{0}\rangle$ (6) $\displaystyle=$
$\displaystyle\Bigl{\langle}\Psi_{0}\Bigr{|}\sum_{d=1}^{2N}\sum_{(i,1)\in\mathcal{N}_{d}}\\!\\!\\!H_{d,i}\,\Bigl{|}\Psi_{0}\Bigr{\rangle}\geq\sum_{d=1}^{2N}\sum_{(i,1)\in\mathcal{N}_{d}}\varepsilon_{d,i}^{(0)}\,.\
$

The energy-spectrum of each cell Hamiltonian $H_{d,i}$ can be expressed in
terms of quantum spin numbers $S_{t}$, $S_{d}$, $m_{t}$, and $m_{d}$ which are
assigned to the composite spin operators
$\mathbf{S}_{t}=\mathbf{S}_{d}+\mathbf{S}_{i,1}$ and $\mathbf{S}_{d}$,
$S^{z}_{t}$, and $S^{z}_{d}$, respectively, as follows,

$\displaystyle\varepsilon_{d,i}$ $\displaystyle=$
$\displaystyle-\frac{3}{8}(J_{1}+J_{2})+\frac{J_{1}}{2}S_{t}(S_{t}+1)$ (7)
$\displaystyle+\left(\frac{J_{2}}{4}-\frac{J_{1}}{2}\right)S_{d}(S_{d}+1)-h\left(m_{t}-\frac{1}{2}m_{d}\right)\,.\quad$

It is straightforward to show that for $h=0$ the eigenstate with quantum spin
numbers $S_{t}=1/2$ and $S_{d}=0$ represents the true ground state of
$H_{d,i}$ whenever $J_{2}/J_{1}>2$. Hence, in this regime
$\epsilon_{d,i}^{(0)}=-\frac{3}{8}J_{2}$. A (weak) finite field then simply
leads to a polarization of the monomer spins. Owing to this fact, the overall
ground state of $H$ for $J_{2}/J_{1}>2$ and in the monomer-dimer (MD) phase is

$\displaystyle|\mathrm{MD}\rangle=\left\\{\begin{array}[]{l}\prod_{i=1}^{N}|\sigma\rangle_{i,1}\otimes\prod_{d=1}^{2N}|s\rangle_{d},\,\sigma\in\\{\uparrow,\downarrow\\},\,h=0\\\
\prod_{i=1}^{N}|\\!\uparrow\rangle_{i,1}\otimes\prod_{d=1}^{2N}|s\rangle_{d},\,h>0\end{array}\right.\
$ (10)

which has the following energy

$\displaystyle E_{\mathrm{MD}}/N=-\frac{3}{2}J_{2}-\frac{h}{2}.$ (11)

Note, that for $h=0$ the MD phase has an extensive ground-state degeneracy,
$2^{N}$, as each of the $N$ monomer spins can be either in the up or down
state. We will examine in Sec. V up to which field strength the MD phase is
actually stable.

The stability condition $J_{2}/J_{1}>2$ of the MD phase at $h=0$ is in
agreement with the results reported previously by Hirose et al. Hirose  _et
al._ (2020). They also verified the presence of the other exact ground state,
referred to as the dimer-tetramer (DT) phase. The DT ground state of the
spin-1/2 Heisenberg antiferromagnet on a diamond square lattice involves the
singlet-dimer states $|s\rangle_{d}$ and singlet-tetramer states
$|t\rangle_{d}$, which are formed between a dimer $d$ and its two neighboring
monomer spins, denoted $(i,1)$ and $(i^{\prime},1)$ in the following:

$\displaystyle|t\rangle_{d}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{3}}(|\\!\uparrow\rangle_{i,1}|\\!\downarrow\uparrow\rangle_{d}|\\!\downarrow\rangle_{i^{\prime},1}+|\\!\downarrow\rangle_{i,1}|\\!\uparrow\downarrow\rangle_{d}|\\!\uparrow\rangle_{i^{\prime},1})$
(12)
$\displaystyle-\frac{1}{2}\left(|\\!\uparrow\rangle_{i,1}|\\!\uparrow\downarrow\rangle_{d}|\\!\downarrow\rangle_{i^{\prime},1}+|\\!\uparrow\rangle_{i,1}|\\!\downarrow\downarrow\rangle_{d}|\\!\uparrow\rangle_{i^{\prime},1}\right.$
$\displaystyle+\left.|\\!\downarrow\rangle_{i,1}|\\!\uparrow\uparrow\rangle_{d}|\\!\downarrow\rangle_{i^{\prime},1}+|\\!\downarrow\rangle_{i,1}|\\!\downarrow\uparrow\rangle_{d}|\\!\uparrow\rangle_{i^{\prime},1}\right)\,.$

In the DT phase, the highly degenerate ground-state manifold corresponds to
the most dense packing of the singlet-tetramer states (12) on the diamond-
decorated square lattice, whereby one cannot accommodate more than ${N}/{2}$
singlet tetramers $|t\rangle_{d}$ on the diamond-decorated square lattice (the
remaining dimers are in the singlet-dimer state $|s\rangle_{d}$). The ground-
state energy in the DT phase is thus given by

$\displaystyle E_{\rm
DT}/N=\frac{3}{2}\varepsilon_{s}+\frac{1}{2}\varepsilon_{t}.$ (13)

Here, $\varepsilon_{s}=-\frac{3}{4}J_{2}$ refers to the energy of the singlet-
dimer state $|s\rangle_{d}$, and $\varepsilon_{t}=-2J_{1}+\frac{J_{2}}{4}$
denotes the energy of the singlet-tetramer state $|t\rangle_{d}$. In order to
obtain the actual stability regions of these two phases for finite fields, we
turn to computational methods.

##  IV Computational Approaches

For our further analysis of the phase diagram of the spin-1/2 Heisenberg
diamond-decorated square lattice as well as its thermodynamic properties, we
have used a combination of various computational approaches. In this section,
we provide some details regarding the application of these different methods
to the model considered here.

###  IV.1 DMRG

The energy of the MD ground state (10) can be obtained from the Hamiltonian
(2) upon fixing all dimers into the singlet state, $S_{d}=0$. We thus expect
that other ground states are obtained from the effective mixed spin-1 and
spin-1/2 Heisenberg model on a Lieb lattice, given by the Hamiltonian (2) for
the particular choice $S_{d}=1$ for all $2N$ dimers $d$. For $h=0$, this
results in the effective Hamiltonian

$\displaystyle
H^{\mathrm{LM}}_{\mathrm{eff}}=J_{1}\sum_{d=1}^{2N}\sum_{(i,1)\in\mathcal{N}_{d}}\mathbf{S}_{d}\cdot\mathbf{S}_{i,1}+\frac{J_{2}}{2}N.$
(14)

In contrast to the case of fixed dimer-singlet states, the Hamiltonian (14)
cannot be solved analytically and we have therefore adopted the DMRG method
implemented in the Algorithms and Libraries for Physics Simulations (ALPS)
project Bauer  _et al._ (2011) in order to find its lowest-energy eigenstates.
For this purpose, we have performed DMRG calculations taking into account up
to 2000 kept states and up to 20 sweeps for lattices with up to $N=36$ unit
cells with periodic boundary conditions. The respective lowest-energy
eigenvalue of the spin-1/2 Heisenberg antiferromagnet on a diamond-decorated
square lattice is given for $h=0$ by the equation:

$\displaystyle E_{\mathrm{LM}}=E_{\rm L}+\frac{J_{2}}{2}N,$ (15)

where $E_{\rm L}$ denotes the ground state energy of the mixed spin-1 and
spin-1/2 Heisenberg model on the corresponding Lieb lattice with $N$ unit
cells at zero magnetic field. According to the Lieb-Mattis theorem Lieb and
Mattis (1962), the lowest-energy eigenstate of the mixed spin-1 and spin-1/2
Heisenberg model on a Lieb lattice in a zero field belongs to the sector with
the total spin given by the absolute value of the difference of the total spin
on two sublattices $S=|S_{A}-S_{B}|$. For the Lieb lattice composed of $N=36$
unit cells we have indeed obtained a ferrimagnetic ground state with total
spin $S=|S_{A}-S_{B}|=|72-18|=54$ and energy $E_{L}=-88.5600047J_{1}$, i.e.,
the ground-state energy $\varepsilon_{\mathrm{L}}=-2.46000J_{1}$ per unit
cell. We note that this value compares well to the value
$\varepsilon_{\mathrm{L}}=-2.46083J_{1}$ for the ground-state energy of the
mixed spin-1 and spin-1/2 Heisenberg antiferromagnet on the Lieb lattice in
the TDL, given in Ref. Hirose  _et al._ (2018).

In order to construct the ground-state phase diagram we have compared the
energies (11),(13) and (15), which were obtained either by analytical or by
numerical calculations of a lattice with $N=36$ unit cells. To study the
magnetization process and thermodynamic quantities in finite magnetic field,
all energies in zero field are shifted by the Zeeman term according to the
formula

$\displaystyle
E(m_{\mathrm{tot}},N,h)=E(m_{\mathrm{tot}},N,h=0)-h\>m_{\mathrm{tot}}\,,\ $
(16)

where $m_{\mathrm{tot}}$ are the eigenvalues of
$S_{\mathrm{tot}}^{z}=\sum_{i=1}^{N}\sum_{\mu=1}^{5}S^{z}_{i,\mu}$. Field-
driven changes of the lowest-energy eigenstates from the sectors with the
total spins $m_{\mathrm{tot}}$ and $m_{\mathrm{tot}}^{\prime}$ are obtained
from Honecker  _et al._ (2000)

$\displaystyle
h=\frac{E(m_{\mathrm{tot}},N,h=0)-E(m_{\mathrm{tot}}^{\prime},N,h=0)}{m_{\mathrm{tot}}-m_{\mathrm{tot}}^{\prime}}.$
(17)

We note that within the LM phase, the ground-state energy in the TDL and
finite $h$ is given by

$E_{\textrm{LM}}/N=\varepsilon_{\mathrm{L}}+\frac{J_{2}}{2}-\frac{3}{2}h,$
(18)

a result that will be useful further below.

###  IV.2 Exact Diagonalization

In addition to DMRG, we also performed exact diagonalizations of the
Hamiltonian $H$ on systems with up to $N_{s}=30$ spins. First, we exploit
conservation of the local spin of each dimer by expressing the Hamiltonian in
the form (2). Thus, the problem boils down to diagonalizing the Hamiltonian
for a given configuration of total dimer spins $S_{d}=0$, $1$ Derzhko  _et
al._ (2007); Honecker  _et al._ (2011); Ohanyan and Honecker (2012); Honecker
_et al._ (2016). In fact, it suffices to perform this computation for each
topologically inequivalent pattern. We have identified the inequivalent
patterns by computer enumeration. For example, for $N_{s}=30$ spins (6 unit
cells), we find 178 inequivalent configurations of the dimer spins $S_{d}$.
Then, we need to take the degeneracy of the corresponding configurations into
account. For example, the configurations with all $S_{d}=0$ or all $S_{d}=1$
are unique, and generally there are $2N_{s}/5$ configurations with exactly one
$S_{d_{i}}=1$. Furthermore, in the case of the $N_{s}=30$ system, there are up
to 120 different realizations of a given pattern with an intermediate number
of dimer triplets.

To diagonalize each of these cases, we first use conservation of
$S^{z}_{\mathrm{tot}}$, as well as spin inversion. We further use $SU(2)$
symmetry to reconstruct the sector with $m_{\mathrm{tot}}=1$ from the other
ones. Thermodynamic quantities such as the specific heat and magnetic
susceptibility can then be computed from the eigenvalues and the associated
quantum numbers. In sectors where we have so many $S_{d}=1$ dimers that the
Hilbert space dimensions become large, we also use the remaining spatial
symmetries of the configuration to further block-diagonalize the system. The
largest matrix to be diagonalized then occurs in the sector with
$m_{\mathrm{tot}}=2$ and for all dimers in the triplet configuration; for
$N_{s}=30$ the resulting maximal dimension is $257304\approx 2.6\cdot 10^{5}$,
which is considerably smaller than the total dimension $2^{30}\approx 10^{9}$
of this system. Still, this significantly exceeds the size of a previous
computation Heidrich-Meisner  _et al._ (2006) where we had used a custom
diagonalization routine Honecker and Schüle (2008), while the present
diagonalization is instead carried out with a recent version of the Intel®
Math Kernel Library.

At the end of this procedure, the full spectrum can be reconstructed for any
value of $J_{2}$, $J_{1}$ and $h$ thanks to the conservation of the total spin
on the dimers and $z$-component of the total spin $m_{\mathrm{tot}}$. Thus, we
can evaluate thermodynamic properties for all $(J_{2}/J_{1},h/J_{1},T/J_{2})$
by post-processing the results of a single diagonalization run for a given
system size $N_{s}$.

###  IV.3 Quantum Monte Carlo

In order to study the thermodynamic properties of the spin-1/2 Heisenberg
antiferromagnet on the diamond-decorated square lattice on system sizes that
extend beyond those accessible to exact diagonalization, we make use of QMC
simulations. In the following we comment on the QMC method that we used for
this purpose.

The SSE QMC method with directed loops Sandvik and Kurkijärvi (1991); Sandvik
(1999); Syljuåsen and Sandvik (2002) offers a highly efficient and unbiased
approach to study quantum spin models. However, introducing geometric
frustration while working in the conventional local spin-$S^{z}$ basis
generally leads to a sign problem, i.e., an exponential drop of the
statistical accuracy of the QMC simulations at low temperatures and large
system sizes Troyer and Wiese (2005); Hangleiter  _et al._ (2020); Hen (2021).
Fortunately, in certain cases, this issue can be eliminated when performing
the QMC simulations in a basis different from the local spin-$S^{z}$ basis.
More specifically, one considers instead appropriate basis states after
decomposing the Hamiltonian into separate terms of few-sites clusters, such as
dimers or trimers Honecker  _et al._ (2016); Alet  _et al._ (2016); Weber  _et
al._ (2022a). The case of dimers can be used to eliminate the sign problem
completely for, e.g., the fully frustrated bilayer model Alet  _et al._
(2016); Ng and Yang (2017); Stapmanns  _et al._ (2018), while a local spin-
trimer basis avoids the sign-problem for the fully frustrated trilayer Weber
_et al._ (2022a) $S=1/2$ antiferromagnet. For the diamond-decorated square
lattice considered here, a finite value of the coupling $J_{2}$ leads to
geometric frustration. We can avoid the associated sign problem that persists
when using the local spin-$S^{z}$ basis, by treating all $J_{2}$-dimer spins
in the spin-dimer basis, while leaving the local $S^{z}$-basis to the monomer
spins. In this combined 5-site basis for each unit cell, the Hamiltonian $H$
can be simulated free of a sign problem, using the SSE approach based on the
abstract operator loop update introduced in Ref. Weber  _et al._ (2022a).
During the operator-loop update of the SSE simulations, binary operators (such
as the bit-wise exclusive-or operation) are used in a binary representation of
the local cluster states. We refer to Refs. Weber  _et al._ (2022a, b) for
further details on this QMC approach. Here, we performed QMC simulations for
systems with values $L$ up to $24$.

##  V Ground-state phase diagram

In the following we describe in detail the ground-state phase diagram of the
spin-1/2 Heisenberg antiferromagnet on the diamond-decorated square lattice,
described by the Hamiltonian $H$, up to high magnetic fields, as obtained from
our combination of DMRG, exact diagonalization as well as analytical results.

Figure 3: Ground-state phase diagram of the spin-1/2 Heisenberg
antiferromagnet on the diamond-decorated square lattice in the
$J_{2}/J_{1}-h/J_{1}$ plane, containing the Lieb-Mattis (LM), dimer-tetramer
(DT), monomer-dimer (MD), spin-canted (SC) and the saturated paramagnetic (PM)
phase. Dashed (solid) lines denote continuous (discontinuous) quantum phase
transitions.

First, we review the phases that appear at zero field $h=0$. The ground-state
phase diagram in this limit has been obtained in previous works Hirose  _et
al._ (2016, 2017, 2018) and it features three distinct phases – LM, DT and MD,
depending on the coupling ratio $J_{2}/J_{1}$, cf. Fig. 2. In the LM regime,
the ground-state energy is minimized by all $J_{2}$-dimers being in the
triplet state, while the monomer spins are oriented predominantly opposite to
the dimer spins. Note that the quantum nature of the LM ferrimagnetic state is
well observable by the reduction of the local magnetizations of the monomer
and dimer spins as shown in Fig. 4. While quantum fluctuations reduce the
magnetization of the monomer spins by approximately 20%, the quantum reduction
of the magnetization of the dimer spins is much more subtle – only about of
5%. In agreement with the LM theorem, both local magnetizations are consistent
with the total magnetization per site of $M/M_{S}=3/5$ (where $M_{S}$ denotes
full saturation) that is inherent to the LM ferrimagnetic phase. For
$J_{2}/J_{1}$ between about 0.974 and 2, the zero-field ground state is the
highly degenerate DT phase, characterized by a dense packing of the singlet-
tetramer states $|t\rangle_{d}$ given by Eq. (12), while the remaining
$J_{2}$-dimers reside in a singlet-dimer state (3). In the highly frustrated
parameter region $J_{2}/J_{1}>2$ the MD phase (10) is realized in the ground
state, which was described in detail in Sec. III.

The ground-state phase diagram including the magnetic field $h$ is shown in
Fig. 3 in the $J_{2}/J_{1}-h/J_{1}$ plane. As one can see from Fig. 3, the LM
phase is stable up to about $h/J_{1}\approx 4$ and spreads out to larger
interaction ratios $J_{2}/J_{1}$ with increasing magnetic field. By contrast,
the MD phase (10) extends towards lower interaction ratio $J_{2}/J_{1}<2$ in
finite magnetic fields as compared to the parameter regime that is accessible
to the variational approach, cf. Sec. III. We also find that for
$J_{2}/J_{1}\geq 4$, the MD phase is stable all the way up to the saturation
field $h_{\mathrm{sat}}=J_{1}+J_{2}$, beyond which the fully polarized,
saturated paramagnetic (PM) regime is entered. On the other hand, the DT phase
narrows quickly for finite magnetic fields and it disappears completely at
$h/J_{1}\approx 0.5$.

In addition to the LM, DT and MD phases, the phase diagram in Fig. 3 exhibits
two high-field phases. Besides the fully polarized, saturated PM regime, we
identify a spin-canted (SC) phase, with a continuously rising magnetization
upon increasing the magnetic field. As shown in Fig. 4, inside the SC phase
the local monomer spins continuously align with the magnetic field upon
increasing the field strength. Initially, the local dimer magnetization
decreases slightly, before it eventually increases to full polarization as
well. Qualitatively, this behavior is well captured by the classical
Heisenberg model of the mixed spin-1 and spin-1/2 model Heisenberg model on
the underlying Lieb lattice, as detailed in App. A: Within the SC phase, the
spins are canted with respect to the magnetic field direction, displaying
biconical structures, resembling those found in, e.g., the classical
anisotropic Heisenberg model at finite magnetic fields Fisher (1975);
Holtschneider  _et al._ (2007). Fig. 3 furthermore shows that the SC phase is
separated from the PM and the LM phase by continuous field-driven quantum
phase transitions. In contrast, all other field-driven phase transitions
between the various ground-state phases are discontinuous.

Figure 4: Local dimer magnetization $\langle S_{d}^{z}\rangle$, local monomer
magnetization $\langle S_{i,1}^{z}\rangle$ and total magnetization $M$
(divided by the saturated magnetization $M_{S}$) at zero temperature, as
functions of magnetic field $h/J_{1}$ for $J_{2}=0$ as obtained from DMRG for
a $L=6$ system. Figure 5: Zero-temperature magnetization curves of the
spin-1/2 Heisenberg antiferromagnet on a diamond-decorated square lattice as
obtained from full exact diagonalization and DMRG method for system sizes
$N_{s}=20$, $N_{s}=30$ and $N_{s}=180$, respectively, for several values of
the interaction ratio $J_{2}/J_{1}$. In all plots, the total magnetization $M$
is scaled with respect to its saturated value $M_{S}$.

The phase diagram and the nature of the transitions in Fig. 3 can be directly
identified from the zero-temperature magnetization curves for the Hamiltonian
$H$. These are shown in Fig. 5 along several vertical cuts through the ground-
state phase diagram. The magnetization curves presented in Fig. 5 were
obtained by two different numerical methods: by exact diagonalization for
$N_{s}=20,30$ and by the DMRG method supplemented with exact analytical
results for $N_{s}=180$, respectively. Overall, the results obtained from both
methods are in excellent agreement, taking into account that within the SC
phase, the smaller size of the system either with $N_{s}=20$ or $N_{s}=30$
leads to more pronounced discretion steps in the stair-case profile of the
magnetization.

More specifically, for an interaction ratio of $J_{2}/J_{1}=0.5$, the
magnetization exhibits an extended 3/5-plateau at low fields, characteristic
of the ferrimagnetic LM phase, followed by a stair-case increase of the
magnetization, which evolves into a continuous magnetization increase in the
TDL, and eventually terminates in the fully saturated PM phase at a magnetic
field of $h/J_{1}=5$ [cf. 5(a)].

Similarly, in agreement with the ground-state phase diagram, there exists a
zero-magnetization plateau for $J_{2}/J_{1}=1.3$ and $J_{2}/J_{1}=1.7$ in
Figs. 5(b) and 5(c) inherent to the gapped DT state (closer inspection
resolves a tiny $1/5$-plateau on the $N_{s}=30$ system in Figs. 5(b) – a
finite-size effect on this particular cluster). For larger values of
$J_{2}/J_{1}>2$, shown in Figs. 5(d), 5(e) and 5(f), the zero-magnetization
plateau disappears. Instead, here, the monomer spins become fully polarized
within the MD phase already for an arbitrarily weak, finite magnetic field.
This results in the immediate onset of the intermediate 1/5-plateau,
characteristic of the MD phase. While for the lower value of $J_{2}/J_{1}=2.5$
the 3/5-plateau of the LM phase, and a subsequent steady increase of the
magnetization in the SC phase can be observed, the magnetization curve for the
higher value of $J_{2}/J_{1}=3.3$ exhibits a discontinuous field-driven
transition from the MD phase (1/5-plateau) into the SC phase. Finally, the
magnetization curve for sufficiently high values of $J_{2}/J_{1}\geq 4$ shows
a direct jump of the magnetization from the 1/5-plateau of the MD phase
towards the fully saturated PM regime.

##  VI Thermal Properties

In the following, we will investigate several aspects of the thermal
properties of the spin-1/2 Heisenberg model on the diamond-decorated square
lattice in the presence of a magnetic field, focusing on the different regions
of the ground-state phase diagram, which were detailed in the previous
section.

###  VI.1 Thermodynamics in the MD regime

We start by investigating the thermodynamic properties in the MD regime. Here,
a simple lattice-gas model Derzhko  _et al._ (2007); Strečka  _et al._ (2017,
2022) can be used to describe the relevant low-energy excitations in the
regime $J_{2}/J_{1}\geq 4$. In addition to a set of free $S=1/2$ spins,
corresponding to the monomer spins $S_{i,1}$ in a magnetic field, this model
contains a lattice gas of hard-core particles that correspond to the dimer-
singlet states on the $J_{2}$-dimers. These particles describe localized
magnons, i.e., a flat band of magnetic excitations, relative to the fully
polarized state Schulenburg  _et al._ (2002); Derzhko  _et al._ (2007);
Strečka  _et al._ (2017, 2022). The lattice-gas model is given, up to a
constant, by the effective Hamiltonian

$H^{\mathrm{MD}}_{\mathrm{eff}}=-h\sum\limits_{i=1}^{N}S_{i}^{z}-\mu\sum\limits_{i=1}^{2N}n_{d}\,,$
(19)

where $n_{d}\in{0,1}$ denotes the local occupation number of the hard-core
particles. A value of $n_{d}=1$ ($n_{d}=0$) corresponds to the presence of a
singlet (triplet) state on the $J_{2}$-dimer $d$. The chemical potential
$\mu=J_{2}-h$ is given by the energy difference between the singlet and lowest
energy triplet state. All thermodynamic properties then follow from the free
energy

$\frac{F}{N}=\frac{J_{2}}{2}-2h-T\ln\Bigl{[}2\cosh\Bigl{(}\frac{h}{2T}\Bigr{)}\Bigr{]}-2T\ln\bigl{[}1+\exp(\mu/T)\bigr{]}\,,$
(20)

where the constant ensures that, at $T=0$, the ground-state energy detailed in
Sec. III is recovered. The ground state corresponds to a fully occupied
(empty) lattice of singlets below (above) the saturation field
$h_{\mathrm{sat}}=J_{1}+J_{2}$. One can approximately describe the
thermodynamic properties in the MD regime also by a spin model that accounts
for all the dimer states (this model corresponds to the limit $J_{1}=0$ of
$H$). However, the above lattice-gas model already turns out to describe the
low-temperature thermodynamics remarkably well. Indeed, Figs. 6 and 7 show
that the lattice-gas model describes the thermodynamic properties at low
temperatures rather accurately, up to $T/J_{1}\approx 0.4$, where the specific
heat starts to show noticeable deviations to the data obtained using exact
diagonalization ($N_{s}=30$) and QMC ($N_{s}=80$), indicating that at higher
temperatures additional excitations become relevant (note that over the whole
regime there is an excellent agreement between the $N_{s}=30$ ED data and the
$N_{s}=80$ data obtained from QMC). From the phase diagram in Fig. 3 we expect
additional excitations to become most relevant in the regime near
$h/J_{1}\approx 5$ at $J_{2}/J_{1}=4$, due to the SC phase. Indeed, in Fig.
7(c), for $T/J_{1}=0.3$, small differences are already resolved in this
magnetic-field range.

Figure 6: Specific heat $C$ divided by the temperature $T/J_{1}$ as a function
of temperature in the MD phase at an interaction ratio of $J_{2}/J_{1}=4$ and
a magnetic field of $h/J_{1}=1$ as obtained from full exact diagonalization,
QMC as well as the effective lattice-gas model.

The magnetization is less susceptible to these additional states, and exhibits
an excellent agreement with the numerical data for all temperatures considered
here. It is noteworthy that in all cases, the results obtained from exact
diagonalization and QMC agree very well with each other, indicating that the
system sizes considered for the exact diagonalization are already
representative of the TDL in this regime. Our analysis therefore indicates
that in addition to the exact analytical results for the MD ground state, also
the low-temperature thermal properties can be understood analytically by means
of a simple effective lattice-gas model given by the Hamiltonian (19).

Figure 7: Thermodynamic properties in the MD phase at interaction ratio
$J_{2}/J_{1}=4$ as a function of magnetic field $h/J_{1}$ obtained from full
exact diagonalization, QMC as well as the effective lattice-gas model. Top
row: Specific heat $C$ for temperatures: (a) $T/J_{1}=0.2$, (b) $T/J_{1}=0.3$,
(c) $T/J_{1}=0.4$. Bottom row: Magnetization $M$, divided by the saturation
magnetization $M_{S}$ for temperatures: (d) $T/J_{1}=0.2$, (e) $T/J_{1}=0.3$,
(f) $T/J_{1}=0.4$

###  VI.2 Thermodynamics in the LM and SC regime

Next, we consider the thermodynamic properties in the LM and SC phase. In the
LM phase, we expect the low-temperature thermodynamics to be governed by the
underlying mixed-spin Lieb lattice. For a quantitative comparison we
considered the case of zero interaction $J_{2}=0$. In Fig. 8, we show the
specific heat $C$ for the diamond-decorated square lattice model as obtained
from ED (for $N_{s}=30$ sites) and QMC as well as for the corresponding mixed-
spin Lieb lattice model. We observe that in both cases, $h=0$ and $h=J_{1}$,
the mixed-spin Lieb lattice model captures the low-temperature asymptotic
behavior, while the behavior differs noticeably at intermediate temperatures.
This deviation is due to additional contributions with singlet configurations
for the diamond-decorated square lattice. These additional states are located
at higher energies, but they have a high density such that they lead to a
relevant contribution to $C$ in the temperature window of Fig. 8.

Another point concerns the strong finite-size effects in Fig. 8(a). These are
due to the $h=0$ ground state being a spin-$3N_{s}/10$ multiplet (compare Sec.
IV.1) such that on a finite-size system, part of the entropy is located at
$T=0$. The case $N_{s}=2000$ ($1200$) should be a good approximation to the
TDL, as is indicated by comparison with the $N_{s}=980$ ($588$) data. Indeed,
these data for $C/T$ in Fig. 8(a) approach a constant for $T\to 0$, as
expected for a ferro- or ferrimagnet in two dimensions, while the activated
low-temperature asymptotics of the $N_{s}=30$ ED data reflects a finite-size
gap of about $0.434J_{1}$.

Applying a magnetic field $h=J_{1}$ opens a gap of the same size in the
excitation spectrum. This leads to activated low-temperature behavior and
negligible finite-size, as the good agreement of the ED and QMC results for
$N_{s}=30$ and $80$ ($18$ and $48$) in Fig. 8(b) shows. As in the case $h=0$,
we again observe significant differences between the diamond-decorated square
lattice model and the mixed-spin Lieb lattice throughout most of Fig. 8(b)
with the exception of the region $T\lesssim 0.2J_{1}$, where we observe the
exponentially activated low-temperature asymptotics. This deviation can again
be attributed to the large number of additional contributions with singlet
configurations for the diamond-decorated square lattice that we already
mentioned in the context of Fig. 8(a).

Figure 8: Specific heat $C$ divided by temperature $T/J_{1}$ as a function of
temperature at magnetic field $h/J_{1}=0$ (top panel) and $h/J_{1}=1$ (bottom
panel). We compare results for the diamond-decorated square lattice model for
$J_{2}=0$, with results on the mixed-spin Lieb lattice. For comparison, the
Lieb lattice data are scaled by a factor of $3/5$ to account for the larger
unit cell on the diamond-decorated square lattice.

Figure 9 presents similar results for $C/T$ at $h/J_{1}=4.5$, i.e., a point in
the SC phase. There is a maximum at temperatures $T>J_{1}$ that, like for Fig.
8(b), is affected by only small finite-size effects. Furthermore, also like in
Fig. 8, the diamond-decorated square lattice model and the mixed-spin Lieb
lattice differ in this high-temperature region. There is a second low-
temperature maximum for $T$ just above $0.1J_{1}$ that is remarkably well
captured by the mixed-spin Lieb lattice. This low-temperature maximum in $C/T$
is affected by stronger finite-size effects, but we checked that the residual
finite-size effects on the $L=20$ data shown in Fig. 9 are not substantial.

We recall that the SC phase is characterized by (quasi-)long-range
antiferromagnetic (XY) order in the plane perpendicular to the magnetic-field
direction at zero (low) temperatures, as discussed in the previous section. At
non-zero temperatures, this leads to a (Berezinskii-)Kosterlitz-Thouless (KT)
transition Berezinskii (1971); Kosterlitz and Thouless (1973); Kosterlitz
(1974) at a finite temperature $T_{\mathrm{KT}}$. Beyond $T_{\mathrm{KT}}$,
the XY quasi-long range order is destroyed by the proliferation of vortex
excitation. We have estimated $T_{\mathrm{KT}}$ using a standard finite-size
scaling analysis of the spin stiffness $\rho_{S}$, as obtained from the spin
winding number fluctuations Weber and Minnhagen (1988); Harada and Kawashima
(1998) (cf. App. B for details). The inset of Fig. 9 shows the temperature
dependence of $\rho_{S}$ for $J_{2}/J_{1}=0$ and $h/J_{1}=4.5$ for increasing
system sizes, exhibiting a drop near $T/J_{1}\approx 0.08$. From a
quantitative analysis Weber and Minnhagen (1988); Harada and Kawashima (1998)
of the QMC data, along the line $h/J_{1}=4.5$, the KT transition temperature
is found to be $T_{\mathrm{KT}}/J_{1}\approx 0.0825$ across the SC regime (cf.
App. B for details). The specific heat $C$ displays a maximum at a temperature
slightly above the KT transition temperature, as typical for the KT
transition, associated to the entropy release from vortex unbinding Chaikin
and Lubensky (1995). This is shown explicitly in Fig. 9 for the interaction
$J_{2}=0$ and the magnetic field $h/J_{1}=4.5$.

Figure 9: Specific heat $C$ divided by the temperature $T/J_{1}$ as a function
of temperature in the SC phase for $J_{2}/J_{1}=0$ and $h/J_{1}=4.5$ as
obtained from QMC. The dashed vertical line denotes the KT transition
temperature $T_{\mathrm{KT}}$, given in the main text. The corresponding data
for the mixed-spin Lieb lattice model is shown for comparison as well, and has
been rescaled by a factor of $3/5$ to account for the larger unit cell on the
diamond-decorated square lattice. The inset shows the temperature dependence
of the spin stiffness of the diamond-decorated square lattice model for the
same parameters for different system sizes, along with the dashed line
indicating $T_{\mathrm{KT}}$.

###  VI.3 Thermal LM-MD phase boundary

Finally, we turn to consider the thermal properties within the parameter
regime where the transition between the LM and MD phases takes place at zero
temperature. As detailed in Sec. V, in the presence of a finite magnetic field
$h\gtrsim 0.5J_{1}$, the LM and MD states at $T=0$ are separated by a direct
discontinuous quantum phase transition line. Across this line, the
$J_{2}$-dimer states change from triplets in the LM phase to singlet states in
the MD phase (the monomer spins are fully polarized along the magnetic field
in the MD phase, while their mean value is reduced due to quantum fluctuations
in the LM phase).

Recently, such discontinuous quantum phase transitions in coupled spin-dimer
and spin-trimer systems were examined in other models Stapmanns  _et al._
(2018); Larrea Jiménez  _et al._ (2021); Weber  _et al._ (2022a), and it was
shown that first-order thermal phase transitions emerge out from the
discontinuous quantum phase transition line, terminating in a line of thermal
critical points. Moreover, these thermal critical points belong to the two-
dimensional Ising universality class, in accord with the binary variable
associated to the presence/absence of a singlet state on the spin dimers (such
as the variable $n_{d}$ introduced above). Here, the LM-MD transition line
offers another realization for such a scenario. We thus examine it in more
detail.

As an example, Fig. 10 shows the mean singlet occupation

$n_{s}=\left\langle\frac{1}{2N}\sum_{d=1}^{2N}n_{d}\right\rangle$ (21)

of the $J_{2}$-dimers as a function of $J_{2}$ along a cut at constant
$h/J_{1}=2.5$ across the LM-MD transition region. At low temperatures, this
quantity exhibits a jump from a value of 0 to a value of 1 upon increasing
$J_{2}/J_{1}$ across the quantum phase transition near $J_{2}/J_{1}\approx
2.5$. For temperatures beyond about $0.3J_{1}$, we instead observe a smooth
variation of $n_{s}$ with increasing $J_{2}$. This already provides indication
for the existence of a low-$T$ discontinuous thermal phase transition line and
its termination in a critical point. While the precise position of the
critical point needs to be extracted from QMC simulations (as detailed below),
the first-order transition line at finite temperature can be estimated by
simply comparing the free energies of both phases, following the approach used
in Refs. Weber  _et al._ (2022a); Stapmanns  _et al._ (2018).

Figure 10: Mean singlet occupation in the vicinity of the LM-MD transition for
$h/J_{1}=2.5$ as obtained from QMC for the system size $L=4$. The red line
shows the first-order transition line obtained from comparing the free
energies of both phases, extended up to the location of the critical point
(symbol), as extracted from a finite-size analysis of the QMC data (see text
for details). In the low-$T$ region (white), the QMC data exhibit large
statistical fluctuations and have been cut off.

Since both the LM and MD phase have a finite excitation gap atop their
respective ground states, we can use the generic estimate of the free energy
for a gapped system at low $T$, given by

$\displaystyle
F=-T\,\ln(Z)=-T\,\ln\left(e^{-E_{0}/T}+e^{-(E_{0}+\Delta)/T}\right)\,,$ (22)

where $E_{0}$ is the ground-state energy and $\Delta$ the excitation gap. At a
fixed magnetic field, the transition line is then obtained from the points
$J_{2}^{c}(T_{c})$, for which the coexistence condition
$F_{\mathrm{LM}}(J_{2}^{c},T^{c})=F_{\mathrm{MD}}(J_{2}^{c},T^{c})$ holds.
Based on Eq. (22), we expect the first-order line not to be vertical, but to
bend towards the phase with the larger excitation gap.

For a quantitative evaluation of the transition line, we require the values of
$E_{0}$ and $\Delta$ in both phases upon approaching the transition point. The
ground-state energies are given by Eq. (18) and Eq. (11) for the LM and MD
phase, respectively. We used exact diagonalization for a system of $N_{s}=30$
sites to extract an estimate for the excitation gaps. The excitation gap
$\Delta$ as a function of the interaction ratio $J_{2}/J_{1}$ for different
magnetic fields $h$ is shown in Fig. 11. We find that upon going from the
magnetic field $h/J_{1}=2$ to $h/J_{1}=2.5$, the excitation gap in the LM
phase becomes larger than that in the MD phase. We thus expect the bending of
the first-order line to change upon increasing the magnetic field. In
particular, for the case of $h/J_{1}=2.5$, considered already in Fig. 10, the
line bends to the left. This is however hardly seen on the scale of the main
panel of Fig. 10. The bending is better seen in the inset, which also shows
the location of the critical point as extracted from further QMC simulations
(as detailed below). Note that based on the free-energy argument, we cannot
determine the location of the critical point, but from the inset of Fig. 10,
we find that its location matches well with the estimated first-order
transition line. The main panel of Fig. 12 shows the specific heat $C$ in the
transition regime. Here, we observe two well pronounced lines of maxima that
expand out from the location of the critical point, very similar to the
behavior observed previously in related systems Stapmanns  _et al._ (2018);
Larrea Jiménez  _et al._ (2021); Weber  _et al._ (2022a).

Figure 11: Excitation gap $\Delta/J_{1}$ as function of the interaction ratio
$J_{2}/J_{1}$ for various magnetic fields $h/J_{1}$ as obtained from exact
diagonalization for the $N_{s}=30$ system. Solid (dashed) lines are used to
denote the excitation gap in the LM (MD) regime.

To accurately locate the critical point, we performed a finite-size scaling
analysis for the fluctuations of the mean singlet occupancy. More
specifically, we consider the corresponding singlet susceptibility Stapmanns
_et al._ (2018),

$\chi_{s}=\frac{\beta}{4N}\left(\left\langle\left(\sum_{d=1}^{2N}{n}_{d}\right)^{2}\right\rangle-\left\langle\sum_{d=1}^{2N}{n}_{d}\right\rangle^{2}\right).$
(23)

In the left panel of Fig. 13, we show this quantity for different system sizes
at a fixed temperature $T/J_{1}=0.32$ across the transition region. The data
exhibit pronounced maxima. Within the two-dimensional Ising universality of
the critical point, the maximum value scales as
$\chi_{s}^{\textrm{max}}\propto L^{7/4}$ at criticality Stapmanns  _et al._
(2018). This property can be used to extract the value of $T_{c}$ from
performing a finite-size scaling analysis of the peak position, as shown in
the upper right panel of Fig. 13, giving $T_{c}/J_{1}=0.315(5)$. From
analyzing the corresponding values of $J_{2}/J_{1}$ of the peak position, cf.
the lower right panel of Fig. 13, we can extract the critical coupling ratio
$(J_{2}/J_{1})_{c}=2.4745(5)$ as well. Together they give the estimated
location of the critical point already shown in Figs. 10 and 12. We also
performed a corresponding analysis at $h/J_{1}=2$. Here, according to the
excitation gaps shown in Fig. 11, we expect the first-order line to bend to
the right instead. This is indeed confirmed by our analysis, cf. the
corresponding data for the specific heat shown in Fig. 14.

Figure 12: Specific heat in the vicinity of the LM-MD transition for
$h/J_{1}=2.5$ as obtained from QMC for the system size with $L=4$. The red
line shows the first-order transition line obtained from comparing the free
energies of both phases, extended up to the location of the critical point
(symbol), as extracted from a finite-size analysis of the QMC data (see text
for details). In the low-$T$ region (white), the QMC data exhibit large
statistical fluctuations and have been cut off. Figure 13: Left panel: Singlet
susceptibility $\chi_{s}$ for a fixed temperature $T/J_{1}=0.32$ across the
transition region for $h/J_{1}=2.5$ and for various system sizes as obtained
from QMC simulations. Right panel: finite-size scaling analysis of the peak
value and its positions to extract the location of the critical point at
$h/J_{1}=2.5$. Figure 14: Specific heat in the vicinity of the LM-MD
transition for $h/J_{1}=2$ as obtained from QMC for the system size with
$L=4$. The red line shows the first-order transition line obtained from
comparing the free energies of both phases, extended up to the location of the
critical point (symbol), as extracted from a finite-size analysis of the QMC
data (see text for details). In the low-$T$ region (white), the QMC data
exhibit large statistical fluctuations and have been cut off.

For the future, it would be interesting to investigate the thermal properties
of this model with respect to several other aspects, such as (i) how the KT
transition lines of the SC phase merge with the discontinuous quantum phase
transition between the SC and MD phases, and (ii) how the thermal properties
of the DT phase can be quantitatively described by effective models of low-
energy excitations, similar to the lattice-gas model for the MD phase. We hope
that our investigations motivate further research on these challenging topics
in the future.

##  VII Conclusions

In this article we considered the spin-1/2 Heisenberg antiferromagnet on the
diamond-decorated square lattice in the presence of a finite magnetic field,
using a combination of analytical arguments and exact diagonalization, density
matrix renormalization group, as well as sign-problem free stochastic series
expansion quantum Monte Carlo simulations.

We identified the ground-state properties at finite magnetic field and mention
here several aspects: (i) the previously identified zero-field Lieb-Mattis
(LM), dimer-tetramer (DT) and monomer-dimer (MD) phases all extend to finite
magnetic fields, with a magnetization 3/5, 0, and 1/5-plateau characterizing
the LM, DT and MD regime, respectively, (ii) at intermediate fields, the DT
phase vanishes and beyond this magnetic-field range, a direct discontinuous
quantum phase transition takes place between the LM and MD phases, (iii) at
high magnetic fields, in addition to the fully saturated paramagnetic phase
(PM), a spin-canted (SC) phase with (quasi)-long-range order emerges.

Additionally, we showed that in the MD regime for $J_{2}/J_{1}\geq 4$, the
low-temperature thermodynamic properties can be well described in terms of a
simple effective lattice-gas model.

Motivated by related results in other quantum spin models, we showed that the
direct, discontinuous quantum phase transition line between the LM and MD
phase extends up to finite temperature, ending in a line of critical points
that belong to the two-dimensional Ising universality class. We furthermore
demonstrated that the slope of the transition line changes sign upon
increasing the strength of the magnetic field.

For the future it would certainly be interesting to extend this analysis to
other phases such as the DT phase, and furthermore to investigate its
excitations in more detail.

###### Acknowledgements.

We acknowledge discussions with Lukas Weber, and support by the Deutsche
Forschungsgemeinschaft (DFG) through Grant No. WE/3649/4-2 of the FOR 1807 and
through RTG 1995, and thank the IT Center at RWTH Aachen University and JSC
Jülich for access to computing time through the JARA Center for Simulation and
Data Science. Part of the ED computations were carried out on the “osaka”
cluster at the Centre de Calcul (CDC) of CY Cergy Paris Université. We
acknowledge support under the Stefanik program for Slovak-France bilateral
projects SK-FR-19-0013 / 45125RC. KK and JS were financially supported under
the grants VEGA 1/0105/20 and APVV-20-0150.

##  Appendix A Ground states of the classical mixed spin-1 and spin-1/2 model
on the Lieb lattice

Figure 15: Angles $\theta/\pi$ of the spins of the classical mixed-spin
Heisenberg model on the Lieb lattice as a function of the magnetic field $h/J$
at zero temperature. Figure 16: Local magnetizations $M_{1}$ and $M_{d}$ of
the two inequivalent sites of the classical mixed-spin Heisenberg model on the
Lieb lattice as well as the total magnetization $M_{\mathrm{tot}}$ (normalized
by the saturated magnetization $M_{S}$) as a function of the magnetic field
$h/J$ at zero temperature.

In this appendix we consider a classical version of the mixed spin-1 and
spin-1/2 Heisenberg model in a magnetic field on the Lieb lattice. To this end
we first note that the quantum Hamiltonian can be written as a sum of bond
operators such that

$\displaystyle
H=\sum\limits_{b=1}^{N_{b}}H_{b}=\sum\limits_{b=1}^{N_{b}}J\,\mathbf{S}_{1}\cdot\mathbf{S}_{d}-\frac{h}{4}S_{1}^{z}-\frac{h}{2}S_{d}^{z},$
(24)

where $\mathbf{S}_{1}$ and $\mathbf{S}_{d}$ denote the spin-1/2 and spin-1
sites respectively and the sum runs over all bonds on the lattice. For a
classical description, we replace the quantum spins by three-dimensional
vectors of length 1/2 and 1, i.e.,

$\mathbf{S}_{1}\rightarrow\frac{1}{2}\begin{pmatrix}\sin\theta_{1}\cos\phi_{1}\\\
\sin\theta_{1}\sin\phi_{1}\\\
\cos\theta_{1}\end{pmatrix},\quad\mathbf{S}_{d}\rightarrow\begin{pmatrix}\sin\theta_{d}\cos\phi_{d}\\\
\sin\theta_{d}\sin\phi_{d}\\\ \cos\theta_{d}\end{pmatrix}.$ (25)

The bond terms in the classical model are then given by

$\displaystyle H_{b}=$
$\displaystyle\frac{J}{2}\left(\sin\theta_{1}\sin\theta_{d}\cos(\phi_{1}-\phi_{d})+\cos\theta_{1}\cos\theta_{d}\right)$
$\displaystyle-\frac{h}{2}\left(\frac{1}{4}\cos(\theta_{1})+\cos(\theta_{d})\right).$
(26)

The ground state of the total classical model is obtained upon minimizing
$H_{b}$ with respect to all four angles
$\theta_{1},\theta_{d},\phi_{1},\phi_{d}$. Since $\theta_{1}$ and $\theta_{d}$
are restricted between $0$ and $\pi$, minimizing with respect to both
$\phi_{d}$ and $\phi_{1}$ yields $\phi_{1}-\phi_{d}=\pi$, i.e., an
antiferromagnetic alignment of neighboring spins transverse to the field
direction. This yields

$\displaystyle H_{b}=$
$\displaystyle\frac{J}{2}\left(-\sin\theta_{1}\sin\theta_{d}+\cos\theta_{1}\cos\theta_{d}\right)$
$\displaystyle-\frac{h}{2}\left(\frac{1}{4}\cos(\theta_{1})+\cos(\theta_{d})\right).$
(27)

Differentiating with respect to $\theta_{1}$ and $\theta_{d}$ gives the
conditions

$\displaystyle\frac{\partial H_{b}}{\partial\theta_{1}}$
$\displaystyle=\frac{J}{2}(-\cos\theta_{1}\sin\theta_{d}-\sin\theta_{1}\cos\theta_{d})+\frac{h}{8}\sin\theta_{1}=0$
(28) $\displaystyle\frac{\partial H_{b}}{\partial\theta_{d}}$
$\displaystyle=\frac{J}{2}(-\sin\theta_{1}\cos\theta_{d}-\cos\theta_{1}\sin\theta_{d})+\frac{h}{2}\sin\theta_{d}=0.$
(29)

Subtracting Eq. (28) from Eq. (29) yields the relation
$\sin\theta_{1}=4\sin\theta_{d}$ which, when reinserted into Eq. (28), finally
leads to

$\displaystyle\left(\frac{h}{8}-\frac{J}{8}\cos\theta_{1}-\frac{J}{2}\sqrt{1-\frac{1}{16}\sin^{2}\theta_{1}}\right)\sin\theta_{1}=0.$
(30)

Based on Eq. (30), one can identify three different regimes, by requiring
either factor to be zero, shown in Fig. 15. We identify first a ferrimagnetic
(FI) regime for $0\leq h/J\leq{}3$, where the spins align in opposite
directions with $\theta_{1}=\pi$ and $\theta_{d}=0$, with a ground-state
energy of

$\displaystyle E_{\mathrm{FI}}/N_{b}=-\frac{J}{2}-\frac{3}{8}h.$ (31)

Next, we identify in the regime $3\leq h/J\leq 5$ a phase in which the spin
directions change continuously upon varying the magnetic field, given by

$\displaystyle\cos\theta_{d}{=}\frac{J}{8h}\left[\left(\frac{h}{J}\right)^{2}{+}15\right],\>\cos\theta_{1}{=}\frac{J}{2h}\left[\left(\frac{h}{J}\right)^{2}{-}15\right].$
(32)

Here, the spins are canted with respect to the direction of the magnetic
field, forming biconical structures with the total magnetization
$M_{tot}/M_{s}=(1/2\cos\theta_{1}+2\cos\theta_{d})/(5/2)=h/(5J)$. Note that at
small fields, $\mathbf{S}_{1}$ is aligned in the opposite direction of the
magnetic field. Upon increasing the magnetic-field strength, however, both
spins align with the magnetic-field direction. Finally, we identify a fully
saturated paramagnetic (PM) phase, where all spins align in direction of the
magnetic field with $\theta_{1}=\theta_{d}=0$ and a ground-state energy

$\displaystyle E_{\mathrm{PM}}/N_{b}=\frac{J}{2}-\frac{5}{8}h.$ (33)

Figure 16 shows the local magnetization of both sites $\mathbf{S}_{1}$ and
$\mathbf{S}_{d}$ as well as the (normalized) total magnetization as a function
of the magnetic field. We note that this classical result is in good
qualitative agreement with the results obtained for the quantum model, where
the local magnetizations within the LM phase are slightly suppressed compared
to their saturated values due to quantum fluctuations.

##  Appendix B Determination of the KT transition

Figure 17: Spin stiffness $\rho_{S}$ for different system sizes $L$ as a
function of temperature $T/J_{1}$ for the diamond-decorated square lattice
model at $J_{2}/J_{1}=0$ and $h/J_{1}=4.5$. The dashed line denotes the
scaling form of the universal jump. The inset shows the quantity $A(T)$ from
the finite-size scaling analysis. The KT transition temperature is denoted by
the dashed vertical line, where $A(T)=1$ holds, obtained using a linear fit
(solid line).

In this appendix, we detail the determination of the KT transition temperature
$T_{\mathrm{KT}}$ within the SC phase. A standard means of identifying
$T_{\mathrm{KT}}$ in $O(2)$-symmetric systems is based on the behavior of the
spin stiffness $\rho_{S}$, which is predicted to exhibit a universal jump of
$\rho_{S}=2\,T_{\mathrm{KT}}/\pi$ at $T_{\mathrm{KT}}$ Nelson and Kosterlitz
(1977). Within the SSE QMC approach, $\rho_{S}$ can be calculated from the
spin winding number fluctuations Pollock and Ceperley (1987); Sandvik (1997);
Caci  _et al._ (2021)

$\rho_{S}=\frac{T}{2A_{\mathrm{uc}}}\bigl{(}\langle W_{x}^{2}\rangle+\langle
W_{y}^{2}\rangle\bigr{)},$ (34)

where $W_{x}$ and $W_{y}$ are the total winding numbers in the orthogonal $x$
and $y$ direction, respectively. Here $A_{\mathrm{uc}}$ is the unit cell area
of the underlying Bravais lattice. For the diamond-decorated square lattice,
$A_{\mathrm{uc}}=1$ holds. To extract $T_{\mathrm{KT}}$ from finite-size QMC
data, we follow the standard approach of Ref. Harada and Kawashima (1998),
which is based on the finite-size scaling form Weber and Minnhagen (1988)

$\frac{\rho_{S}\,\pi}{2\,T}=A(T)\,\biggl{(}1+\frac{1}{2\,\log(L/L_{0}(T))}\biggr{)}$
(35)

that holds exactly at the transition point with $A(T_{\mathrm{KT}})=1$. We
fitted this finite-size dependence to the data for different temperatures,
using $A(T)$ and $L_{0}(T)$ as fit parameters. This allows us to accurately
estimate $T_{\mathrm{KT}}$, where $A(T_{\mathrm{KT}})=1$ holds. Our results
from this approach are shown in Fig. 17, and we obtain from this analysis an
estimate of $T_{\mathrm{KT}}/J=0.08248(3)$ at $J_{2}/J_{1}=0$ and
$h/J_{1}=4.5$. Performing the same analysis at different ratios of
$J_{2}/J_{1}$ within the SC phase for $h/J_{1}=4.5$, we obtain similar values
for the KT transition temperature.

## References

  * Richter  _et al._ (2004) J. Richter, J. Schulenburg, and A. Honecker, “Quantum magnetism in two dimensions: From semi-classical Néel order to magnetic disorder,” in  _Quantum Magnetism_, edited by U. Schollwöck, J. Richter, D. J. J. Farnell, and R. F. Bishop (Springer Berlin Heidelberg, Berlin, Heidelberg, 2004) pp. 85–153.
  * Balents (2010) L. Balents, Nature 464, 199 (2010).
  * Lacroix  _et al._ (2011) C. Lacroix, P. Mendels, and F. Mila,  _Introduction to Frustrated Magnetism: Materials, Experiments, Theory_, Springer Series in Solid-State Sciences, Vol. 164 (Springer Berlin Heidelberg, Berlin, Heidelberg, 2011).
  * Diep (2013) H. T. Diep,  _Frustrated Spin Systems, 2nd edition_ (World Scientific, Singapore, 2013).
  * Sriram Shastry and Sutherland (1981) B. Sriram Shastry and B. Sutherland, Physica B+C 108, 1069 (1981).
  * Albrecht and Mila (1996) M. Albrecht and F. Mila, Europhys. Lett. 34, 145 (1996).
  * Miyahara and Ueda (1999) S. Miyahara and K. Ueda, Phys. Rev. Lett. 82, 3701 (1999).
  * Miyahara and Ueda (2003) S. Miyahara and K. Ueda, J. Phys.: Condens. Matter 15, R327 (2003).
  * Kageyama  _et al._ (1999) H. Kageyama, K. Yoshimura, R. Stern, N. V. Mushnikov, K. Onizuka, M. Kato, K. Kosuge, C. P. Slichter, T. Goto, and Y. Ueda, Phys. Rev. Lett. 82, 3168 (1999).
  * Kageyama  _et al._ (2000a) H. Kageyama, K. Onizuka, Y. Ueda, M. Nohara, H. Suzuki, and H. Takagi, J. Exp. Theor. Phys. 90, 129 (2000a).
  * Kageyama  _et al._ (2000b) H. Kageyama, M. Nishi, N. Aso, K. Onizuka, T. Yosihama, K. Nukui, K. Kodama, K. Kakurai, and Y. Ueda, Phys. Rev. Lett. 84, 5876 (2000b).
  * Lemmens  _et al._ (2000) P. Lemmens, M. Grove, M. Fischer, G. Güntherodt, V. Kotov, H. Kageyama, K. Onizuka, and Y. Ueda, Phys. Rev. Lett. 85, 2605 (2000).
  * Gaulin  _et al._ (2004) B. D. Gaulin, S. H. Lee, S. Haravifard, J. P. Castellan, A. J. Berlinsky, H. A. Dabkowska, Y. Qiu, and J. R. D. Copley, Phys. Rev. Lett. 93, 267202 (2004).
  * Zayed  _et al._ (2014) M. E. Zayed, C. Rüegg, T. Strässle, U. Stuhr, B. Roessli, M. Ay, J. Mesot, P. Link, E. Pomjakushina, M. Stingaciu, K. Conder, and H. M. Rønnow, Phys. Rev. Lett. 113, 067201 (2014).
  * Zayed  _et al._ (2017) M. E. Zayed, C. Rüegg, J. Larrea J., A. M. Läuchli, C. Panagopoulos, S. S. Saxena, M. Ellerby, D. F. McMorrow, T. Strässle, S. Klotz, G. Hamel, R. A. Sadykov, V. Pomjakushin, M. Boehm, M. Jiménez-Ruiz, A. Schneidewind, E. Pomjakushina, M. Stingaciu, K. Conder, and H. M. Rønnow, Nature Phys. 13, 962 (2017).
  * Guo  _et al._ (2020) J. Guo, G. Sun, B. Zhao, L. Wang, W. Hong, V. A. Sidorov, N. Ma, Q. Wu, S. Li, Z. Y. Meng, A. W. Sandvik, and L. Sun, Phys. Rev. Lett. 124, 206602 (2020).
  * Larrea Jiménez  _et al._ (2021) 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, C. Rüegg, P. Corboz, H. M. Rønnow, and F. Mila, Nature 592, 370 (2021).
  * Onizuka  _et al._ (2000) K. Onizuka, H. Kageyama, Y. Narumi, K. Kindo, Y. Ueda, and T. Goto, J. Phys. Soc. Jpn. 69, 1016 (2000).
  * Kodama  _et al._ (2002) K. Kodama, M. Takigawa, M. Horvatić, C. Berthier, H. Kageyama, Y. Ueda, S. Miyahara, F. Becca, and F. Mila, Science 298, 395 (2002).
  * Sebastian  _et al._ (2008) S. E. Sebastian, N. Harrison, P. Sengupta, C. D. Batista, S. Francoual, E. Palm, T. Murphy, N. Marcano, H. A. Dabkowska, and B. D. Gaulin, Proc. Natl. Acad. Sci. 105, 20157 (2008).
  * Takigawa and Mila (2011) M. Takigawa and F. Mila, “Magnetization plateaus,” in  _Introduction to Frustrated Magnetism: Materials, Experiments, Theory_, edited by C. Lacroix, P. Mendels, and F. Mila (Springer Berlin Heidelberg, Berlin, Heidelberg, 2011) pp. 241–267.
  * Takigawa  _et al._ (2013) M. Takigawa, M. Horvatić, T. Waki, S. Krämer, C. Berthier, F. Lévy-Bertrand, I. Sheikin, H. Kageyama, Y. Ueda, and F. Mila, Phys. Rev. Lett. 110, 067210 (2013).
  * Jaime  _et al._ (2012) M. Jaime, R. Daou, S. A. Crooker, F. Weickert, A. Uchida, A. E. Feiguin, C. D. Batista, H. A. Dabkowska, and B. D. Gaulin, Proc. Natl. Acad. Sci. 109, 12404 (2012).
  * Matsuda  _et al._ (2013) Y. H. Matsuda, N. Abe, S. Takeyama, H. Kageyama, P. Corboz, A. Honecker, S. R. Manmana, G. R. Foltin, K. P. Schmidt, and F. Mila, Phys. Rev. Lett. 111, 137204 (2013).
  * Haravifard  _et al._ (2016) S. Haravifard, D. Graf, A. E. Feiguin, C. D. Batista, J. C. Lang, D. M. Silevitch, G. Srajer, B. D. Gaulin, H. A. Dabkowska, and T. F. Rosenbaum, Nat. Comm. 7, 11956 (2016).
  * Corboz and Mila (2013) P. Corboz and F. Mila, Phys. Rev. B 87, 115144 (2013).
  * Corboz and Mila (2014) P. Corboz and F. Mila, Phys. Rev. Lett. 112, 147203 (2014).
  * Wessel  _et al._ (2018) S. Wessel, I. Niesen, J. Stapmanns, B. Normand, F. Mila, P. Corboz, and A. Honecker, Phys. Rev. B 98, 174432 (2018).
  * Wietek  _et al._ (2019) A. Wietek, P. Corboz, S. Wessel, B. Normand, F. Mila, and A. Honecker, Phys. Rev. Research 1, 033038 (2019).
  * Honecker  _et al._ (2022) A. Honecker, L. Weber, P. Corboz, F. Mila, and S. Wessel, J. Phys.: Conf. Ser. 2207, 012032 (2022).
  * Weihong  _et al._ (1999) Z. Weihong, C. J. Hamer, and J. Oitmaa, Phys. Rev. B 60, 6608 (1999).
  * Müller-Hartmann  _et al._ (2000) E. Müller-Hartmann, R. R. P. Singh, C. Knetter, and G. S. Uhrig, Phys. Rev. Lett. 84, 1808 (2000).
  * Hida (1992) K. Hida, J. Phys. Soc. Jpn. 61, 1013 (1992).
  * Sandvik and Scalapino (1994) A. W. Sandvik and D. J. Scalapino, Phys. Rev. Lett. 72, 2777 (1994).
  * Sommer  _et al._ (2001) T. Sommer, M. Vojta, and K. W. Becker, Eur. Phys. J. B 23, 329 (2001).
  * Wang  _et al._ (2006) L. Wang, K. S. D. Beach, and A. W. Sandvik, Phys. Rev. B 73, 014431 (2006).
  * Alet  _et al._ (2016) F. Alet, K. Damle, and S. Pujari, Phys. Rev. Lett. 117, 197203 (2016).
  * Ng and Yang (2017) K.-K. Ng and M.-F. Yang, Phys. Rev. B 95, 064431 (2017).
  * Stapmanns  _et al._ (2018) J. Stapmanns, P. Corboz, F. Mila, A. Honecker, B. Normand, and S. Wessel, Phys. Rev. Lett. 121, 127201 (2018).
  * Richter  _et al._ (2006) J. Richter, O. Derzhko, and T. Krokhmalskii, Phys. Rev. B 74, 144430 (2006).
  * Derzhko  _et al._ (2007) O. Derzhko, J. Richter, A. Honecker, and H.-J. Schmidt, Low Temp. Phys. 33, 745 (2007).
  * Hirose  _et al._ (2017) Y. Hirose, A. Oguchi, and Y. Fukumoto, J. Phys. Soc. Jpn. 86, 014002 (2017).
  * Hirose  _et al._ (2016) Y. Hirose, A. Oguchi, and Y. Fukumoto, J. Phys. Soc. Jpn. 85, 094002 (2016).
  * Hirose  _et al._ (2018) Y. Hirose, S. Miura, C. Yasuda, and Y. Fukumoto, AIP Advances 8, 101427 (2018).
  * Hirose  _et al._ (2020) Y. Hirose, A. Oguchi, and Y. Fukumoto, Phys. Rev. B 101, 174440 (2020).
  * Note (1) This lattice appeared in the context of the CuO2 planes of the high-$T_{c}$ superconductors Mattis (1988) and thanks to extensions in Ref. Lieb (1989) is now known as “Lieb” lattice. Actually, an Ising version of this model had appeared in the literature already significantly earlier Fisher (1960a, b), but we will stick to established terminology.
  * Lieb and Mattis (1962) E. Lieb and D. Mattis, J. Math. Phys. 3, 749 (1962).
  * White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
  * White and Huse (1993) S. R. White and D. A. Huse, Phys. Rev. B 48, 3844 (1993).
  * Schollwöck (2005) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
  * Sandvik and Kurkijärvi (1991) A. W. Sandvik and J. Kurkijärvi, Phys. Rev. B 43, 5950 (1991).
  * Sandvik (1999) A. W. Sandvik, Phys. Rev. B 59, R14157 (1999).
  * Syljuåsen and Sandvik (2002) O. F. Syljuåsen and A. W. Sandvik, Phys. Rev. E 66, 046701 (2002).
  * Honecker  _et al._ (2016) A. Honecker, S. Wessel, R. Kerkdyk, T. Pruschke, F. Mila, and B. Normand, Phys. Rev. B 93, 054408 (2016).
  * Bose (1989) I. Bose, J. Phys.: Condens. Matter 1, 9267 (1989).
  * Bose (1990) I. Bose, J. Phys.: Condens. Matter 2, 5479 (1990).
  * Bose (1992) I. Bose, Phys. Rev. B 45, 13072 (1992).
  * Bauer  _et al._ (2011) B. Bauer, L. D. Carr, H. G. Evertz, A. Feiguin, J. Freire, S. Fuchs, L. Gamper, J. Gukelberger, E. Gull, S. Guertler, A. Hehn, R. Igarashi, S. V. Isakov, D. Koop, P. N. Ma, P. Mates, H. Matsuo, O. Parcollet, G. Pawłowski, J. D. Picon, L. Pollet, E. Santos, V. W. Scarola, U. Schollwöck, C. Silva, B. Surer, S. Todo, S. Trebst, M. Troyer, M. L. Wall, P. Werner, and S. Wessel, J. Stat. Mech.: Theor. Exp. 2011, P05001 (2011).
  * Honecker  _et al._ (2000) A. Honecker, F. Mila, and M. Troyer, Eur. Phys. J. B 15, 227 (2000).
  * Honecker  _et al._ (2011) A. Honecker, S. Hu, R. Peters, and J. Richter, J. Phys.: Condens. Matter 23, 164211 (2011).
  * Ohanyan and Honecker (2012) V. Ohanyan and A. Honecker, Phys. Rev. B 86, 054412 (2012).
  * Heidrich-Meisner  _et al._ (2006) F. Heidrich-Meisner, A. Honecker, and T. Vekua, Phys. Rev. B 74, 020403 (2006).
  * Honecker and Schüle (2008) A. Honecker and J. Schüle, Parallel Computing: Architectures, Algorithms and Applications. NIC Series 38, 271 (2008).
  * Troyer and Wiese (2005) M. Troyer and U.-J. Wiese, Phys. Rev. Lett. 94, 170201 (2005).
  * Hangleiter  _et al._ (2020) D. Hangleiter, I. Roth, D. Nagaj, and J. Eisert, Science Advances 6, eabb8341 (2020).
  * Hen (2021) I. Hen, Phys. Rev. Research 3, 023080 (2021).
  * Weber  _et al._ (2022a) L. Weber, A. Honecker, B. Normand, P. Corboz, F. Mila, and S. Wessel, SciPost Phys. 12, 054 (2022a).
  * Weber  _et al._ (2022b) L. Weber, N. Caci, and S. Wessel, Phys. Rev. B 106, 035141 (2022b).
  * Fisher (1975) M. E. Fisher, AIP Conference Proceedings 24, 273 (1975).
  * Holtschneider  _et al._ (2007) M. Holtschneider, S. Wessel, and W. Selke, Phys. Rev. B 75, 224417 (2007).
  * Strečka  _et al._ (2017) J. Strečka, J. Richter, O. Derzhko, T. Verkholyak, and K. Karľová, Phys. Rev. B 95, 224415 (2017).
  * Strečka  _et al._ (2022) J. Strečka, T. Verkholyak, J. Richter, K. Karl’ová, O. Derzhko, and J. Schnack, Phys. Rev. B 105, 064420 (2022).
  * Schulenburg  _et al._ (2002) J. Schulenburg, A. Honecker, J. Schnack, J. Richter, and H.-J. Schmidt, Phys. Rev. Lett. 88, 167207 (2002).
  * Berezinskii (1971) V. Berezinskii, Sov. Phys. JETP 32, 493 (1971).
  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1181 (1973).
  * Kosterlitz (1974) J. Kosterlitz, J. Phys. C 7, 1046 (1974).
  * Weber and Minnhagen (1988) H. Weber and P. Minnhagen, Phys. Rev. B 37, 5986 (1988).
  * Harada and Kawashima (1998) K. Harada and N. Kawashima, J. Phys. Soc. Jpn. 67, 2768 (1998).
  * Chaikin and Lubensky (1995) P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, 1995).
  * Nelson and Kosterlitz (1977) D. R. Nelson and J. M. Kosterlitz, Phys. Rev. Lett. 39, 1201 (1977).
  * Pollock and Ceperley (1987) E. L. Pollock and D. M. Ceperley, Phys. Rev. B 36, 8343 (1987).
  * Sandvik (1997) A. W. Sandvik, Phys. Rev. B 56, 11678 (1997).
  * Caci  _et al._ (2021) N. Caci, L. Weber, and S. Wessel, Phys. Rev. B 104, 155139 (2021).
  * Mattis (1988) D. C. Mattis, Phys. Rev. B 38, 7061 (1988).
  * Lieb (1989) E. H. Lieb, Phys. Rev. Lett. 62, 1201 (1989).
  * Fisher (1960a) M. E. Fisher, Proc. R. Soc. London, Ser. A 254, 66 (1960a).
  * Fisher (1960b) M. E. Fisher, Proc. R. Soc. London, Ser. A 256, 502 (1960b).