# Unveiling nonmagnetic phase and many-body entanglement in
two-dimensional random quantum magnets Sr2CuTe1-xWxO6

Dian Wu<EMAIL_ADDRESS>Institute of Physics, École Polytechnique Fédérale de
Lausanne (EPFL), CH-1015 Lausanne, Switzerland  Fan Yang<EMAIL_ADDRESS>Institute of Physics, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015
Lausanne, Switzerland  Giuseppe Carleo<EMAIL_ADDRESS>Institute of
Physics, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne,
Switzerland

###### Abstract

We apply a random-plaquette $J_{1}$-$J_{2}$ model on the square lattice to
capture the physics of a series of spin-$1/2$ Heisenberg antiferromagnet
compounds Sr2CuTe1-xWxO6. With the input of experimentally relevant coupling
strengths, our exact diagonalization (ED) study probes the ground state
properties beyond previous linear spin-wave approach. An intermediate range of
$x\in[0.08,0.55]$ is identified for a nonmagnetic phase without the long-range
Néel or stripe order. The absence of both valence-bond-glass order and spin-
glass non-ergodic dynamics renders its nature intriguing. Deep inside this
phase around $x=0.3$, we observe signatures potentially linked to randomness-
induced short-range spin-liquid-like (SLL) states, including close to zero
spin-freezing parameter, vanishing spin-spin correlation beyond nearest
neighbors, almost uniform static spin structure factor, as well as a broad
tail in the dynamical spin structure factor. The nonmagnetic phase also
features multipartite entanglement in the ground state witnessed by quantum
Fisher information (QFI), which exhibits universal scaling behaviors at
quantum critical points.

##  I Introduction

Quantum spin liquids (QSLs) [1, 2, 3], as ground states of disordered quantum
spins characteristic of vanishing long-range magnetic order and massive many-
body entanglement, have renewed theoretical and experimental interests [4, 5,
6, 7, 8, 9] in recent years. From a few exactly solvable examples, such as the
Kitaev $\mathbb{Z}_{2}$ spin liquids in the spin-$1/2$ honeycomb model [10]
with generalizations to ladder geometries [11, 12], the triangle-honeycomb
[13] and three-dimensional lattices [14], as well as extensions to higher-spin
analogues [15], it is well known that certain types of QSL feature fractional
excitations of quasiparticles, including anyons [16] and Majorana zero modes
[17, 18], which can be adapted into braiding and fusion mechanisms [10, 19]
leading to potential fault-tolerant quantum computation [20, 21].

Pursuing QSLs in experimental platforms proves to be challenging. Quantum
simulators [22, 23, 24, 25, 26, 27] offer precise control of spin exchange
interactions, yet remain to be scaled up to large system sizes. Quantum
materials of frustrated magnets [28, 4, 29, 8, 9, 30], on the other hand,
deviate from the ideal frustration scenarios set up in simple models, and more
often than not involve intrinsic inhomogeneity or quenched randomness, which
results in disorder on sites or bonds. Meanwhile, motivated by theoretical
concerns on finding more methods to directly observe many-body entanglement as
independent evidence for quantum critical phenomena in strongly correlated
systems [31, 32, 33, 34], there arise increasing interests and efforts in
probing entanglement properties in quantum materials [35, 36, 37, 38, 39, 40],
and establishing links with quantum information science [41, 42].

Figure 1:  (a) The random-plaquette $J_{1}$-$J_{2}$ model of Sr2CuTe1-xWxO6.
The Heisenberg AFM coupling strengths are obtained from quantum chemistry
calculations [43]: $J_{1}=7.4$, $J_{2}=8.3$, $J^{\prime}_{1}=0.7$,
$J^{\prime}_{2}=0.1$, $J^{\prime\prime}_{1}=0.3$, in the unit of meV. Each red
(blue) plaquette is dominated by $J_{1}$ ($J_{2}$) interactions with a Te (W)
ion at the center. (b) Phase diagram of the ground state via ED as a function
of $x$, the mixing ratio of W or the probability of $J_{2}$-dominated
plaquette. At the two ends, two magnetically ordered phases exhibit the Néel
order and the strip order respectively, forming sublattices of spins shown in
different colors. A nonmagnetic phase appears in the middle region
$x\in[0.08,0.55]$.

In this work, we focus on the spin-$1/2$ square-lattice Heisenberg
antiferromagnet: the double perovskite Sr2CuTe1-xWxO6 [44, 45, 46, 47, 48, 49,
50, 43, 51, 52], whose degrees of frustration and randomness can be tuned by
the mixing ratio $x$ of Te and W ions. At $x=0$ and $1$, the compounds
Sr2CuTeO6 and Sr2CuWO6 are illustrated by the first two panels of Fig. 1 (a).
They are dominated by nearest-neighbor $J_{1}$ and next-nearest-neighbor
$J_{2}$ Heisenberg antiferromagnetic (AFM) interactions and also have weak
frustrated interactions $J^{\prime}_{2}$ and $J^{\prime}_{1}$, thus exhibit
the Néel order and the stripe order respectively, as in Fig. 1 (b). While
nuclear magnetic resonance (NMR) experiments [44, 45] suggest a SLL state in
between, numerical attempts to reveal its nature become challenging both in
the choice of effective spin models, and in the absence of lattice symmetry
and translational symmetry due to disorder. Preliminary ED computations on the
random-bond $J_{1}$-$J_{2}$ model point to a random-singlet phase and a spin-
glass phase [53], while density matrix renormalization group (DMRG) studies
find evidence only of the random-singlet state [54]. The latter is in
agreement with quantum Monte Carlo (QMC) simulations on the random $J$-$Q$
model [55] with $J$ the Heisenberg exchange and $Q$ the multispin interaction.
The two effective models studied so far are remotely linked to the material
and the role of the mixing ratio $x$ remains obscure.

Here, we focus on a different and more realistic microscopic description [43],
the random-plaquette $J_{1}$-$J_{2}$ model depicted in Fig. 1 (a). In the two-
dimensional square lattice, each plaquette can be $J_{2}$-dominated (W
counterion) with probability $x$, and otherwise $J_{1}$-dominated (Te
counterion). When two different kinds of plaquette are adjacent, i.e., when a
Te and a W counterions encounter, the strong $J_{1}$ bond on the shared edge
is replaced by a weaker $J^{\prime\prime}_{1}$ bond. From the linear spin-wave
approach [43, 52], the random-plaquette model is found to faithfully reproduce
the inelastic neutron scattering spectra, the magnetic excitations, as well as
the neutron diffraction patterns of Sr2CuTe1-xWxO6. Moreover, Monte Carlo
simulations show that even with classical spins, this model can explain the
extreme suppression of the Néel order at small $x=0.025$, consistent with muon
spin resonance ($\mu$SR) measurements [51]. Yet, a computational simulation of
the quantum model is still lacking.

Throughout this work, we present a comprehensive ED study of the model with
quantum spins. We take Heisenberg AFM coupling strengths obtained from quantum
chemistry calculations [43]:

$(J_{1},\ J_{2},\ J^{\prime}_{1},\ J^{\prime}_{2},\
J^{\prime\prime}_{1})=(7.4,\ 8.3,\ 0.7,\ 0.1,\ 0.3),$ (1)

in the unit of meV. Remarkably, as shown in Fig. 1 (b) and Fig. 2 (a)–(d),
from the order parameters we identify a nonmagnetic ground state for
Sr2CuTe1-xWxO6 in a wide range of $x\in[0.08,0.55]$, matching well with the
NMR observations [44, 45, 46] that there is no magnetic order over $19$ mK for
$x\in[0.1,0.6]$. On top of this, the universal scaling exponent of the QFI
density [32, 36, 37] extracted from the static spin structure factor vanishes
near the phase boundaries, confirming the accuracy of our phase diagram.

Figure 2:  Finite-size scaling of the order parameters $m^{2}_{\text{N}}$,
$m^{2}_{\text{str}}$, $q_{\text{EA}}$, and $q_{\text{VBG}}$ used to
characterize the phases. In the left column, for each value of $x$, we measure
the order parameters with system sizes $N=8,12,16,20,24,28,32$ and extrapolate
to the thermodynamic limit $N\to\infty$. The error bars are standard errors
over random instances and of the extrapolated intercepts. In the right column,
we do the extrapolation for all values of $x$ and plot the results as the
black curves. The shaded regions are standard errors of the intercepts. The
vertical dashed lines indicate the nonmagnetic phase $x\in[0.08,0.55]$. In
(b), when $m_{\text{N}}^{2}$ crosses $0$, $x=0.081\pm 0.012$. In (d), when
$m_{\text{str}}^{2}$ crosses $0$, $x=0.546\pm 0.020$. In (f), over
$x\in[0.08,0.55]$, the average value of $q_{\text{EA}}$ is $0.017\pm 0.004$.

The nature of the intermediate phase is distinct from normal candidates
previously found in quantum magnets subject to quenched randomness, including
the random-singlet state [53, 54, 55, 56, 57, 58], the valence-bond glass [46,
59, 60], and the spin glass [53, 36, 37]. Instead, its main features are more
akin to short-range disordered SLL states in two dimensions, which have an
extremely small spin-freezing parameter and a broad tail in the dynamical spin
structure factor. Intriguingly, deep inside the nonmagnetic phase around
$x=0.3$ where these characteristics become the most predominant, the system
displays zero averaged spin-spin correlation beyond nearest neighbors, and
almost uniformly distributed static spin structure factor, leading to a
uniform QFI density at $T=0$.

The paper is structured as follows. In Section II, based on ED computations,
we perform finite-size scaling of the order parameters and determine the phase
diagram for the random-plaquette $J_{1}$-$J_{2}$ model. We characterize the
emergent short-range disordered SLL state by spin-spin correlations, spin
relaxations, as well as spin structure factors. In Section III, we focus on
the QFI and reveal its universal scaling behaviors at quantum phase
transitions. In Section IV, we advance the understanding of the phase
formation in Sr2CuTe1-xWxO6 qualitatively through the randomness-induced
frustration index.

##  II Properties of the model

Before proceeding to study physical observables, we first briefly introduce
our reweighting method to compute the observables, which enables us to obtain
the disorder-averaged ground-state expectation values with an arbitrary mixing
ratio $x\in[0,1]$. For each system size $N=8,12,16,20,24,28,32$, we construct
a lattice whose geometry allows both Néel and stripe orders under periodic
boundary conditions (PBC). Then for each $k=0,1,\ldots,N$, we randomly
generate $M_{k}$ Hamiltonian instances, each with exactly $k$
$J_{2}$-dominated plaquettes. As $x$ is the probability for a plaquette to be
$J_{2}$-dominated in the disorder average, any physical observable $O$ is
averaged under a binomial distribution, denoted by $[\,\cdots]_{J}$:

$\displaystyle\left[\\!\expectationvalue*{O}\right]_{J}(N,x)=\sum_{k}B(k,N,x)\bar{O}_{k}(N),$
(2)
$\displaystyle\bar{O}_{k}(N)=\frac{1}{M_{k}}\sum_{i}\expectationvalue*{O}_{i},$
(3)

where the binomial distribution takes the form
$B(k,N,x)=\binom{N}{k}x^{k}(1-x)^{N-k}$, and the binomial coefficient reads
$\binom{N}{k}=N!/(k!\,(N-k)!)$. In $\bar{O}_{k}(N)$, the expectation
$\expectationvalue*{O}_{i}$ acts on the ground state of the random instance
$i$. Without reweighting, $\bar{O}_{k}(N)$ alone is an estimator of the
uniform average over all random instances with a fixed $k$, under which
circumstance $x$ can only take discrete values $k/N$ with a finite system size
$N$. Our reweighting method makes it possible to vary $x$ continuously and
more precisely determine the phase boundaries, as shown in the right column of
Fig. 2. More details on the ED computation, including choices of lattice
geometries, reweighting, and estimation of uncertainties in the average and
the extrapolation, can be found in Appendix A.

We comment on the symmetries of the ground state. It inherits the SU$(2)$
symmetry from Heisenberg exchanges, leading to zero magnetization
$\expectationvalue*{S^{z}_{j}}=0$ for each site $j$. The ground state also
respects spin inversion symmetry and selects the even parity sector for all
lattices we use, regardless of the value of $x$ and the configuration of
random plaquettes. Translational symmetry and $D_{4}$ lattice symmetry, by
contrast, are only present at $x=0$ and $1$.

###  II.1 Order parameters and phase diagram

Now, we determine the phase diagram of the random-plaquette $J_{1}$-$J_{2}$
model (see Fig. 1) from four order parameters shown in Fig. 2. While the Néel
and the stripe order parameters set the boundaries of two magnetically ordered
phases when $x$ is close to $0$ and $1$ respectively, the spin-freezing
parameter and the valence-bond-glass parameter suggest that the nonmagnetic
region $x\in[0.08,0.55]$ cannot stabilize a random-singlet state nor a
valence-bond glass.

First, let us take a closer look at two magnetically ordered phases. We choose
the following definitions [53] for the Néel and the stripe order parameters:

$\displaystyle m_{\text{N}}^{2}$
$\displaystyle=\frac{8}{N(N+4)}\left[\sum_{\alpha}\sum_{i,j\in\alpha}\\!\expectationvalue*{\bm{S}_{i}\cdot\bm{S}_{j}}\right]_{J},$
(4) $\displaystyle m_{\text{str}}^{2}$
$\displaystyle=\frac{4}{N(N+4)}\left[\sum_{\nu}\,\sum_{\alpha_{\nu}}\sum_{i,j\in\alpha_{\nu}}\\!\\!\expectationvalue*{\bm{S}_{i}\cdot\bm{S}_{j}}\right]_{J}.$
(5)

In $m_{\text{N}}^{2}$, $\alpha$ denotes each of the two sublattices formed by
the Néel order, as shown in different colors in Fig. 1 (b), and $i,j$ run over
all sites in the sublattice $\alpha$. In $m_{\text{str}}^{2}$, $\nu$ denotes
each of the two directions of stripes, giving rise to two patterns of
sublattices $\alpha_{\nu}$. Figure 1 (b) depicts the horizontal stripes, and
vertical otherwise. The disorder average $[\,\cdots]_{J}$ is evaluated using
our reweighting method in Eq. 2. The normalization factors in Eqs. 4 and 5 are
chosen according to relevant magnetic orderings of the classical Heisenberg
antiferromagnet. For a perfect Néel order, the sum inside $[\,\cdots]_{J}$ in
$m_{\text{N}}^{2}$ containing both $i=j$ and $i\neq j$ amounts to
$2\left(\frac{N}{2}\cdot\frac{3}{4}+\frac{N}{2}(\frac{N}{2}-1)\cdot\frac{1}{4}\right)=N(N+4)/8$.
The factor in $m_{\text{str}}^{2}$ is further multiplied by $2$ taking into
account two types of stripes $\nu$. We then perform a finite-size scaling with
$8\leq N\leq 32$, as shown in Fig. 2 (a)–(d). It turns out that both order
parameters exhibit $N^{-1/2}$ scaling behaviors, consistent with the modified
spin-wave analysis for Heisenberg antiferromagnet on finite square lattices
[61, 62]:

$\expectationvalue*{\bm{S}_{i}\cdot\bm{S}_{j}}=(-1)^{\alpha(i,j)}\left(c_{0}+\frac{c_{1}}{|\bm{r}_{i}-\bm{r}_{j}|}\right)^{2},\quad
i\neq j,$ (6)

where both $c_{0}$ and $c_{1}$ are constant given a finite $N$. The sign
$(-1)^{\alpha(i,j)}=1$ if the sites $i$ and $j$ belong to the same sublattice
$\alpha$, and $-1$ otherwise. Equation 6 leads directly to the scaling
relation:

$m_{\text{N/str}}^{2}=\frac{a}{\sqrt{N}}+b,\quad b\geq 0,$ (7)

where higher order corrections in $O(1/N)$ are neglected. In the thermodynamic
limit $N\to\infty$, if the long-range magnetic order exists,
$m_{\text{N}}^{2}$ and $m_{\text{str}}^{2}$ become positive constants. From
the non-negative extrapolation values of $m_{\text{N}}^{2}$ and
$m_{\text{str}}^{2}$ in Fig. 2 (b) and (d), we identify that in the ground
state, the Néel phase survives from the randomness when $x\in[0,0.08]$ and the
stripe phase survives when $x\in[0.55,1]$. The two quantum phase transition
points crossing zero are resolved as $x_{c,1}=0.081\pm 0.012$,
$x_{c,2}=0.546\pm 0.020$. The extremely narrow region of the Néel order
compared with the stripe order can be explained by reviewing the microscopic
structure of the mixed compound in Fig. 1 (a). When two different kinds of
plaquette encounter, one strong $J_{1}$ bond is replaced by a much weaker
$J^{\prime\prime}_{1}$ bond on the shared edge, while the strong $J_{2}$ bonds
in the $J_{2}$-dominated plaquette remain unaffected.

The consistency of our nonmagnetic phase in absence of long-range order for
$x\in[0.08,0.55]$ with the experimental results $x\in[0.1,0.6]$ [44, 45, 46]
is crucial. It validates the methodology of the current work, including the
Hamiltonian we start with in Eq. 1, our ED computations, as well as our
reweighting method in Eq. 2 and extrapolation method. Next, we proceed to
explore the nature of the nonmagnetic phase by testing two potential
candidates: the random-singlet state and the valence-bond glass.

The random-singlet state can be characterized by the the spin-freezing
parameter [53, 54], also known as the Edwards–Anderson (EA) order parameter
[63, 64], defined by

$q_{\text{EA}}^{2}=\frac{1}{N^{2}}\left[\sum_{i,j}\expectationvalue*{\bm{S}_{i}\cdot\bm{S}_{j}}^{2}\right]_{J},$
(8)

where $\sum_{i}$ sums over a total number of $N$ sites. $q_{\text{EA}}$
detects any type of static spin order, including long-range magnetic orders
and the spin-glass order. In Fig. 2 (e)–(f), $q_{\text{EA}}$ scales with
$N^{-1/2}$ and extrapolates to positive values as $N\to\infty$ for all $x$,
with the nonmagentic region $x\in[0.08,0.55]$ locating a plateau around its
minimum. The inset in (f) shows that the error bars of the extrapolation in
the nonmagnetic phase do not extend to zero, and the average value of the
plateau stays at $0.017\pm 0.004$. Our finding agrees with the linear spin-
wave analysis [52] of the same model in Eq. 1 that the degree of frustration
arising from weak Heisenberg interactions
$J^{\prime}_{1},J^{\prime}_{2},J^{\prime\prime}_{1}$ is not strong enough to
stabilize an all-to-all disordered random-singlet state in Sr2CuTe1-xWxO6.

The valence-bond glass (VBG), on the other hand, can be characterized by an
order parameter that replaces the spin operator in $q_{\text{EA}}$ with the
dimer operator [65]:

$\displaystyle
q_{\text{VBG}}^{2}=\frac{1}{4N^{2}}\left[\sum_{i,j}\sum_{p,q}\expectationvalue*{B_{i}^{p}B_{j}^{q}}^{2}\right]_{J},$
(9) $\displaystyle
B_{i}^{p}=\bm{S}_{i}\cdot\bm{S}_{i+p}-\expectationvalue*{\bm{S}_{i}\cdot\bm{S}_{i+p}},$
(10)

where $B_{i}^{p}$ encodes the dimer-dimer correlation [66] on a nearest-
neighbor bond along the unit vector $p$, and $p$ can be either in $x$ or $y$
direction, i.e., $\sum_{i}\sum_{p}$ sums over a total number of $2N$ nearest-
neighbor bonds. As shown in Fig. 2 (g)–(h), $q_{\text{VBG}}$ also exhibits
$N^{-1/2}$ scaling behavior, but extrapolates to negative values as
$N\to\infty$ for all $x$. We have also tested another dimer operator that sums
over both nearest and next-nearest bonds and obtained the same behavior.
Therefore, the putative valence-bond-glass state from the specific heat
measurement at low temperatures [46], is not observed by our ED results of the
random-plaquette model at $T=0$.

###  II.2 Short-range spin-spin correlations and non-glassy dynamics in the
nonmagnetic phase

Figure 3:  Averaged spin-spin correlation $C_{r}$ of the system with $N=32$
spins. The length of the nearest-neighbor bond is set to $1$. The error bars
include standard errors over random instances and over all spin pairs $(i,j)$
sharing the same $r$ in an instance. At $x=0$ and $1$, the zero error bar is
caused by the absence of randomness and the restored $D_{4}$ lattice symmetry.
Deep inside the nonmagnetic phase around $x=0.3$, the spin-spin correlation
vanishes beyond nearest neighbors.

A small but finite $q_{\text{EA}}$ when $N\to\infty$ could indicate a weak
spin-glass order. However, as we will show in this section through spin-spin
correlations over space and time domains, see Figs. 3 and 4, the spin glass is
also not a candidate for the nonmagnetic phase.

First, let us examine the spin-spin correlation in real space,

$C_{r}=\left[\expectationvalue*{\bm{S}_{i}\cdot\bm{S}_{j}}\right]_{J,\bm{r}},\quad
r=|\bm{r}_{i}-\bm{r}_{j}|.$ (11)

Besides the disorder average $[\,\cdots]_{J}$ over all random instances in Eq.
2, we also take the average over all spin pairs $(i,j)$ sharing the same
distance $r$ in each random instance [67]. Figure 3 shows that while the Néel
(stripe) order is indicated by the dominant negative correlations at nearest-
neighbor $r=1$ (next-nearest-neighbor $r=\sqrt{2}$), the long-range spin-spin
correlations decrease quickly when $x$ enters the nonmagnetic phase. In
particular, deep inside this phase around $x=0.3$, we observe zero spin-spin
correlation beyond nearest neighbors. The remaining AFM spin-spin correlations
on nearest neighbors contribute to the non-zero extrapolation of spin-freezing
parameter $q_{\text{EA}}$ in Fig. 2 (e)–(f), and it is incompatible with the
random domain wall picture of the spin-glass state in the random-bond
$J_{1}$-$J_{2}$ model. According to Ref. [53], the latter would host a
considerable amount of positive ferromagnetic (FM) correlations on nearest-
neighbor spins across the domain walls between two distinct stripe-ordered
states. In our model, the randomness-induced frustration plays a more
important role, of which we leave the discussion to Section IV.

Figure 4:  Temporal auto-correlation $C_{t}$ of a spin in the system with
$N=32$ spins. The shaded regions are standard errors over random instances.
The dynamics is ergodic as the auto-correlation relaxes to zero quickly in the
nonmagnetic phase $x\in[0.08,0.55]$.

Next, we study the temporal auto-correlation of one spin in the ground state,
which can probe the ergodicity breaking in finite-size glassy systems [68]:

$C_{t}=\left[\expectationvalue*{\\{S^{z\dagger}(t),S^{z}(0)\\}}\right]_{J},$
(12)

where $S^{z\dagger}(t)=e^{iHt}S^{z}(0)e^{-iHt}$, and $C_{t}\in\mathbb{R}$
since $S^{z\dagger}(t)=S^{z}(t)$. It is noted that the choice of the spin
$S^{z}$ in the system is arbitrary after averaging over random instances, and
the $z$ orientation is also arbitrary as the model is $\text{SU}(2)$
symmetric. In Fig. 4, the two magnetically ordered phases at $x=0$ and $1$
exhibit oscillating patterns that do not decay at large $t$. By contrast, in
the nonmagnetic phase $x\in[0.08,0.55]$, the auto-correlation relaxes to zero
quickly. We do not observe the non-ergodic behavior in $C_{t}$ that is
characteristic of a spin glass, which would entail a long-lived plateau at
large $t$.

###  II.3 Static and dynamical spin structure factors

Figure 5:  Static spin structure factor $S_{\text{sta}}$ of the system with
$N=32$ spins. From Eq. 19, it also displays the QFI density through the
relation $f_{\text{Q}}(\bm{k})=\frac{4}{3}S_{\text{sta}}(\bm{k})$.

After invalidating three candidates, namely random-singlet state, valence-bond
glass, and spin glass, we proceed to unveil nontrivial SLL signatures of the
nonmagnetic phase through spin structure factors.

The static spin structure factor defined by

$S_{\text{sta}}(\bm{k})=\frac{1}{N}\left[\sum_{i,j}\expectationvalue*{\bm{S}_{i}\cdot\bm{S}_{j}}\cos\left(\bm{k}\cdot(\bm{r}_{i}-\bm{r}_{j})\right)\right]_{J},$
(13)

is resolved in the ground state via ED in Fig. 5, in close resemblance to the
experimental data probed by diffuse polarized neutron diffraction at $T=1.5$ K
[52]. As $x$ increases from $0$ (Néel order) to $1$ (stripe order), the peak
of $S_{\text{sta}}$ moves from $(\pi,\pi)$ to $(\pi,0)$. Meanwhile, Fig. 5
shows a more detailed feature that the peak becomes less sharp once $x$ leaves
the ordered phases and approaches the phase boundaries at $x=0.1$ and $0.5$,
which also agrees with the experimental data. Interestingly, around $x=0.3$,
$S_{\text{sta}}$ becomes almost uniform in a large portion of the Brillouin
zone (BZ) without any visible peak. It indicates a strongly disordered state,
absent of any long-range magnetic order.

Figure 6:  Dynamical spin structure factor $S_{\text{dyn}}$ of the system with
$N=32$ spins, exhibiting broad peaks at $\bm{k}=(\pi,\pi)$ and $(\pi,0)$ in
the nonmagnetic phase $x\in[0.08,0.55]$. The shaded regions are standard
errors over random instances.

Furthermore, as shown from the dynamical spin structure factor in Fig. 6, the
coexistence of a peak centered at a small frequency $\omega$ and a broad tail
at large $\omega$ also supports a SLL state [53, 65] in the nonmagnetic region
around $x=0.3$. In our calculation, we adopt the following definition for the
dynamical spin structure factor:

$\displaystyle S_{\text{dyn}}(\bm{k},\omega)$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}\left[\expectationvalue*{S^{z\dagger}_{\bm{k}}(t)S^{z}_{\bm{k}}(0)}\right]_{J}\\!e^{i\omega
t}\differential t$ (14) $\displaystyle=-\lim_{\eta\to
0}\left[\frac{1}{\pi}\imaginary\expectationvalue*{S^{z\dagger}_{\bm{k}}\frac{1}{\omega+E_{0}-H+i\eta}S^{z}_{\bm{k}}}\right]_{J},$

with $E_{0}$ the ground state energy. The time-domain integral is computed
using continued fraction expansion [69] and smoothed with a sufficiently small
damping factor $\eta=0.1$ meV [53, 65]. The Fourier transform of a spin
operator is taken to be

$S^{z}_{\bm{k}}=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}S_{i}^{z}e^{-i\bm{k}\cdot\bm{r}_{i}}.$
(15)

Figure 6 shows that when the mixing ratio $x$ increases in the nonmagnetic
region, the peak of $S_{\text{dyn}}$ moves in the BZ from the center
$(\pi,\pi)$ to the edge $(\pi,0)$. In the frequency domain, the peak moves
towards smaller $\omega$ as the randomness $x$ increases at
$\bm{k}=(\pi,\pi)$, while at $(\pi,0)$ it stays rather fixed near
$\omega\approx 2$ meV. Most importantly, taking into account the flatness of
the static spin structure factor in the nonmagnetic phase in Fig. 5, the broad
peak of $S_{\text{dyn}}$ extends naturally to other wave vectors in the BZ.
This behavior can be used to distinguish from magnon-like excitations which
introduce peaks only at high-symmetry points of the BZ [65]. We have also
tested finite-size effects and verified that for a fixed $x\in[0.08,0.55]$ at
both wave vectors, the peak reaches lower frequency as the system size
increases.

##  III Quantum Fisher information in 2D random quantum magnets

In this section, we introduce the general relation between the QFI and the
spin structure factors, and demonstrate how to use the QFI density as a
witness of multipartite entanglement in the random-plaquette $J_{1}$-$J_{2}$
model.

###  III.1 Multipartite entanglement witness

Let us start from a generic equilibrium mixed state at the temperature $T$.
Its density matrix has the form
$\rho=\sum_{\lambda}p_{\lambda}\ket{\lambda}\bra{\lambda}$ with the
probability $p_{\lambda}=\frac{1}{Z}e^{-E_{\lambda}/T}$ and the partition
function $Z=\sum_{\lambda}e^{-E_{\lambda}/T}$, where we set $\hbar=k_{B}=1$.
Given a generator $O=\sum_{i=1}^{N}O_{i}$ consisting of local operators
$O_{i}$ that are normalized with unit spectrum width, the QFI is defined by

$f_{\text{Q}}(T)=2\sum_{\lambda,\lambda^{\prime}}\frac{(p_{\lambda}-p_{\lambda^{\prime}})^{2}}{p_{\lambda}+p_{\lambda^{\prime}}}\left|\\!\matrixelement{\lambda}{O}{\lambda^{\prime}}\right|^{2},$
(16)

with $p_{\lambda}+p_{\lambda^{\prime}}>0$. Physically, through the quantum
Cramér–Rao bound [70], $f_{\text{Q}}$ quantifies the maximal precision of a
parameter $\theta$ one can reach after $M$ independent measurements:
$\operatorname{Var}\theta\geq 1/(Mf_{\text{Q}})$. It is a lower bound for the
variance of parameter $\theta$ associated with the generator $O$. Under each
measurement, the system undergoes a unitary transformation:
$\rho^{\prime}=e^{-i\theta O}\rho e^{i\theta O}$. The QFI has been
successfully measured in a single spin system, or a solid-state qubit [71].

Going to multispin systems of higher complexity, an exact relation can be
established [32, 36, 37] between the QFI density $f_{\text{Q}}=f_{\text{Q}}/N$
and the dynamical spin structure factor $S_{\text{dyn}}$:

$f_{\text{Q}}(\bm{k},T)=4\int_{0}^{\infty}\\!\\!\tanh(\frac{\omega}{2T})\left(1-e^{-\frac{\omega}{T}}\right)S_{\text{dyn}}(\bm{k},\omega,T)\,\differential\omega,$
(17)

with the generator chosen to be
$O(\bm{k})=\sqrt{N}S^{z\dagger}_{\bm{k}}=\sum_{i=1}^{N}S_{i}^{z}e^{i\bm{k}\cdot\bm{r}_{i}}$.
Here, $S_{\text{dyn}}$ is obtained from Eq. 14 by applying
$\expectationvalue*{S^{z\dagger}_{\bm{k}}(t)S^{z}_{\bm{k}}(0)}=\Tr[\rho\,S^{z\dagger}_{\bm{k}}(t)S^{z}_{\bm{k}}(0)]$,
which acts on the thermal state when $T>0$ and on the ground state at $T=0$.
Equation 17 expresses the entanglement in an explicit form of quantum
fluctuations. A detailed derivation is provided in Appendix B.

Notably, a quantum state is $(m+1)$-partite entangled when the QFI density
satisfies the inequality [32, 37, 72, 73]:
$f_{\text{Q}}>m(h_{\text{max}}-h_{\text{min}})^{2}$, where $m$ is an integer
and $h_{\text{max}}$ ($h_{\text{min}}$) is the largest (smallest) eigenvalue
of the local generator $O_{i}$. In our case,
$\expectationvalue*{O^{\dagger}_{i}O_{i}}=\expectationvalue*{(S_{i}^{z})^{2}}=1/4$
and $h_{\text{max}}=-h_{\text{min}}=1/2$. Therefore, $(m+1)$-partite
entanglement is witnessed in a multispin system if

$f_{\text{Q}}>m.$ (18)

At the same time, Eq. 18 sets the criterion according to which the quantum
Cramér–Rao bound falls below the classical limits [74].

Figure 7:  QFI density $f_{\text{Q}}$ and multipartite entanglement of the
system at $T=0$ with $N=32$ spins. The vertical dashed lines indicate the
nonmagnetic phase $x\in[0.08,0.55]$. The shaded regions are standard errors
over random instances.

Next, we focus on the features of the QFI in the ground state of the random-
plaquette model. By taking the limit $T\to 0$ in Eq. 17, we can verify that

$f_{\text{Q}}(\bm{k})=4\int_{0}^{\infty}S_{\text{dyn}}(\bm{k},\omega)\,\differential\omega=\frac{4}{3}S_{\text{sta}}(\bm{k}).$
(19)

The extra factor $1/3$ comes from one of the three $\text{SU}(2)$-symmetric
components
$\bm{S}_{i}\cdot\bm{S}_{j}=\sum_{\alpha=x,y,z}S_{i}^{\alpha}S_{j}^{\alpha}$.
Equivalently, we can also evaluate the QFI of a pure ground state
$\ket{\psi}$. The density matrix becomes $\rho=\ket{\psi}\bra{\psi}$, and the
case $p_{\psi}=1$, $p_{\psi^{\prime}}=0$ contributes twice to $f_{\text{Q}}$
in Eq. 16. Applying the identity
$\ket{\psi^{\prime}}\bra{\psi^{\prime}}=\mathds{1}-\ket{\psi}\bra{\psi}$, we
obtain the same relation

$f_{\text{Q}}(\bm{k})=4\left[\expectationvalue*{S^{z}_{\bm{k}}S^{z\dagger}_{\bm{k}}}-|\expectationvalue*{S^{z}_{\bm{k}}}|^{2}\right]_{J}=\frac{4}{3}S_{\text{sta}}(\bm{k}).$
(20)

It is noted that $\expectationvalue*{S^{z}_{\bm{k}}}=0$ respecting the
$\text{SU}(2)$ symmetry.

In Fig. 7, we quantify multipartite entanglement witnessed by $f_{\text{Q}}$
at the two ordering wave vectors using Eqs. 18 and 19. The nonmagnetic phase
is short-range entangled, featuring $2$\- to $4$-partite entanglement on a
finite-size lattice of $N=32$ spins. By contrast, the Néel and the stripe
phases possess much higher entanglement at their respective ordering wave
vectors. This long-range entanglement is also reflected in the sharp peaks of
$S_{\text{sta}}(\bm{k})$ in Fig. 5. When the mixing ratio $x$ increases from
$0$ to $1$, the peak of $f_{\text{Q}}$ moves from $\bm{k}=(\pi,\pi)$ to
$(\pi,0)$ following the same result as $S_{\text{sta}}$. Around $x=0.3$ deep
inside the nonmagnetic phase, $f_{\text{Q}}$ also becomes uniform and the
system displays bipartite entanglement over the BZ.

###  III.2 Universal scaling at quantum critical points

Figure 8:  Scaling of QFI density $f_{\text{Q}}$ at $T=0$. (a) (c) Linear
regression of $\ln f_{\text{Q}}=\nu_{\text{Q}}\ln N+c_{0}$ at specific values
of $\bm{k}$ and $x$, where $\nu_{\text{Q}}$ is the scaling exponent and
$c_{0}$ is a constant. (b) (d) The slope $\nu_{\text{Q}}$ as a function of
$x$. The shaded regions are standard errors of $\nu_{\text{Q}}$. The
horizontal dashed line indicates $\nu_{\text{Q}}=0$, and the vertical dashed
lines indicate the nonmagnetic phase $x\in[0.08,0.55]$. In (a) (b), only the
results of $x<0.6$ are shown, where $f_{\text{Q}}$ is large enough to
numerically estimate the slope. In (b), when $x_{\text{c}}=0.55$,
$\nu_{\text{Q}}=0.02\pm 0.05$. In (d), when $x_{\text{c}}=0.08$,
$\nu_{\text{Q}}=-0.02\pm 0.06$.

Another interesting feature observed from Fig. 7 is that in the finite-size
system with $N=32$, if we choose an ordering wave vector and tunes $x$ towards
the other ordered phase, there is always no multipartite entanglement to be
witnessed at the phase boundary, as $f_{\text{Q}}<1$ at $\bm{k}=(\pi,\pi)$,
$x=0.55$ and $\bm{k}=(\pi,0)$, $x=0.08$. In the following, we will demonstrate
through the scaling analysis of $f_{\text{Q}}$ that in the thermodynamic
limit, the critical points for $f_{\text{Q}}=1$ before entering a separable
(unentangled) region coincide with the phase boundaries of the random-
plaquette model determined by the order parameters.

At $T=0$, we start from a general scaling relation between the QFI density and
the system size:

$\ln f_{\text{Q}}(\bm{k})=\nu_{\text{Q}}(\bm{k})\ln N+c_{0},$ (21)

where $\nu_{\text{Q}}$ denotes the scaling exponent and $c_{0}$ is a constant.
If $\nu_{\text{Q}}=0$, then $f_{\text{Q}}$ converges to a finite value when we
take the thermodynamic limit $N\to\infty$. The finite-size ED computations on
our model lead to Fig. 8. At $\bm{k}=(\pi,\pi)$, $\nu_{\text{Q}}$ is positive
when $x<0.55$ and vanishes as $\nu_{\text{Q}}=0.02\pm 0.05$ around the
critical point $x_{\text{c}}=0.55$. Thus when $N\to\infty$, $f_{\text{Q}}$
converges to a finite value $f_{\text{Q}}\approx 1$ only when $x$ is around
the critical point. The same behaviour of $f_{\text{Q}}$ is observed at
$\bm{k}=(\pi,0)$, which exhibits a vanishing exponent $\nu_{\text{Q}}=-0.02\pm
0.06$ around the critical point $x_{\text{c}}=0.08$.

To conclude, we find that the QFI density is capable to detect the quantum
phase transitions at $T=0$ in the random-plaquette model. Given an ordering
wave vector $\bm{k}$, the scaling exponent of $f_{\text{Q}}$ with respect to
the system size vanishes across the transition point $x_{\text{c}}$:

$\left.\nu_{\text{Q}}(\bm{k})\right|_{x=x_{\text{c}}}\approx 0,$ (22)

which signifies the other magnetically ordered phase. It is noted that the
scaling exponent of $f_{\text{Q}}$ also vanishes at topological phase
transitions, e.g. in the Kitaev superconducting chain [32]. Through Eq. 22, we
extend the universality of this phenomenon to the order-disorder phase
transitions in 2D random quantum magnets.

##  IV Statistics of randomness-induced frustration

Figure 9:  (a) An example of a graph containing only the large couplings
$J_{1}$ and $J_{2}$ with $N=16$ spins and one $J_{2}$-dominated plaquette. Two
examples of frustrated cycles of spin couplings are highlighted in green and
orange respectively. This graph becomes non-frustrated if at least two edges
are removed, therefore it has frustration index $N_{\text{F}}=2$. (b) The
normalized frustration index $n_{\text{F}}=N_{\text{F}}/N$ as a function of
$x$, averaged over random instances with $N=32$ spins. The green region
indicates the nonmagnetic phase $x\in[0.08,0.55]$. The maximum of
$n_{\text{F}}$ is around $x=0.3$, as indicated by the vertical dashed line.
(c) Distribution of nearest-neighbor spin correlations
$\expectationvalue*{\bm{S}_{i}\cdot\bm{S}_{j}}$ with $N=32$ spins. The results
of $x=0.1,0.3,0.5$ are highlighted with thick curves. When $x=0$ and $1$, the
distribution is fully concentrated at
$\expectationvalue*{\bm{S}_{i}\cdot\bm{S}_{j}}=-0.340$ and $0$ respectively,
as indicated by the vertical dashed lines. The probability density function
(PDF) of the distribution is estimated from random instances using Silverman’s
method [75].

Extra insights can be gained on the role of the mixing ratio $x$ in
Sr2CuTe1-xWxO6 by introducing the frustration index of a random instance.
Illustrated in Fig. 9 (a), for each random instance of our model in Eq. 1, we
consider a graph whose edges only contain the large couplings $J_{1}$ and
$J_{2}$, and omit the small ones $J^{\prime}_{1}$, $J^{\prime}_{2}$, and
$J^{\prime\prime}_{1}$. The frustration index $N_{\text{F}}$ is defined to be
the minimal number of edges one needs to remove to make the graph
geometrically non-frustrated (containing only even-length cycles of AFM
couplings). In this simple example, we have $N_{F}=2$. In generic cases, the
method to compute it is equivalent to finding the ground state of a classical
AFM Ising model on the graph, and to solving other well-studied combinatorial
optimization problems such as max-cut [76].

In the random-plaquette model, it is more convenient to consider the
normalized frustration index $n_{\text{F}}=N_{\text{F}}/N$, which converges in
the thermodynamic limit $N\to\infty$. Shown in Fig. 9 (b), when $n_{\text{F}}$
is small ($<0.1$), the system exhibits the Néel or the stripe order. Each
frustrated edge raises the ground state energy approximately by
$2J|S|(|S|+1)$, where $J=J_{1}$ or $J_{2}$ depending on the position of the
frustrated edge, and $|S|=1/2$. Remarkably, when $x$ enters the nonmagnetic
region $[0.08,0.55]$, a high plateau of $n_{\text{F}}$ appears with a peak
around $x=0.3$, where it has so far been identified with most significant
disordered SLL signatures. We can thus gain a phenomenological understanding
on the origin of the nonmagnetic phase in Sr2CuTe1-xWxO6: The mixing ratio $x$
induces additional frustration, which plays a more important role in driving
the system into a short-range disordered SLL phase than the originally present
weak couplings $J^{\prime}_{1}$, $J^{\prime}_{2}$, and $J^{\prime\prime}_{1}$.

Moreover, in Fig. 9 (c) we show the distribution of nearest-neighbour spin
correlations of the random-plaquette model at varying $x$ for a system with
$N=32$ spins. The lower and the upper bounds
$\expectationvalue*{\bm{S}_{i}\cdot\bm{S}_{j}}\in[-3/4,1/4]$ are given by pure
spin singlet and triplet pairs. The two sharp peaks at
$\expectationvalue*{\bm{S}_{i}\cdot\bm{S}_{j}}=-0.340$ and $0$ appear in the
Néel and the stripe phases respectively. There is a smooth transfer of
probability density between them as $x$ changes. Interestingly, in the
nonmagnetic region $x\in[0.08,0.55]$, the distribution spans much broader, and
at $x=0.3$ there appear two peaks of almost equal height. Compared with the
random-bond $J_{1}$-$J_{2}$ model [53], we do not observe any major peak of
$\expectationvalue*{\bm{S}_{i}\cdot\bm{S}_{j}}\in[-0.2,-0.1]$ as in random-
singlet states, and there is also no visible peak at
$\expectationvalue*{\bm{S}_{i}\cdot\bm{S}_{j}}>0$ as in spin-glass states.
This distinguishes our nonmagnetic state from these two candidates.

##  V Discussion

To summarize, in this work we perform extensive ED studies of the random-
plaquette $J_{1}$-$J_{2}$ model with experimentally relevant parameters for
the material Sr2CuTe1-xWxO6. We resolve the phase boundaries of a nonmagnetic
state $x\in[0.08,0.55]$ by the finite-size scaling of various order
parameters. The intermediate phase is absent of long-range Néel and stripe
orders, valence-bond glass order, as well as glassy dynamics. Yet, it still
exhibits short-range spin-spin correlations and a broad tail in
$S_{\text{dyn}}$, which are characteristic of disordered SLL states. Deep
inside the nonmagnetic phase around $x=0.3$, the randomness-induced
frustration index reaches the maximum, bringing a series of non-trivial
features including almost uniform $S_{\text{sta}}$ (or equivalently,
$f_{\text{Q}}$) over the BZ. The broad peaks in $S_{\text{dyn}}$ are also
stretched away from the two ordering wave vectors $\bm{k}=(\pi,\pi)$ and
$(\pi,0)$, marking the difference from magnon-like excitations.

Although the bond frustration from the small couplings
$J^{\prime}_{1},J^{\prime}_{2},J^{\prime\prime}_{1}$ is not strong enough to
support an all-to-all disordered random-singlet state as $q_{\text{EA}}>0$,
the frustration index in the model increases with structured randomness, thus
favoring a disordered SLL state. Finally, from the analysis of the QFI, the
nonmagnetic phase is multipatite entangled, and the vanishing universal
scaling exponent of the QFI density with respect to $N$ across phase
transitions shows the robustness of our phase boundaries in the thermodynamic
limit.

In future, one may generalize the calculation of the QFI to thermal
equilibrium mixed states for $T>0$, using Eq. 17 and the Hams–de Raedt method
[77]. A perfect scaling collapse at low temperatures [32] is expected from
vanishing critical exponent at $T=0$, and will extend the critical-$T$ regime
of the phase diagram in Fig. 1 (b). It is also possible to study the QFI far
from equilibrium in our model using small-size ED clusters, e.g. $N=10$ spins
[78].

Going beyond ED, one may also apply variational Monte Carlo (VMC) to explore
2D random quantum magnets. The recently developed neural quantum states [79,
80], especially the transformer architectures [81, 82], can be expressive
variational ansatzes to capture the rich complexity of the physical system.
There have been successful attempts at applying neural quantum states to
approximate the QSL state in the regular $J_{1}$-$J_{2}$ Heisenberg model on
the square lattice [66, 83, 84, 85], despite the nature of Dirac spin liquid
is still under debate [86]. Given the flexibility of neural networks, the
access to larger system sizes is envisioned, which enables one to further
resolve the dispersion of the dynamical spin structure factor and the spin
excitation gaps, casting more light on the nature of the nonmagnetic state.

## Acknowledgements

We thank Federico Becca, Emil J. Bergholtz, Ellen Fogh, Yuchi He, Riccardo
Rossi, and Luciano L. Viteritti for discussions. We especially thank Henrik M.
Rønnow for pointing our attention to the disordered model studied in this
work. Support from the Swiss National Science Foundation is acknowledged under
Grant No. 200336.

##  Appendix A Finite-size scaling on frustrated random magnets

###  A.1 Choice of finite-size lattices

We choose a series of lattices with $N\leq 32$ to obtain enough data points
from ED and extrapolate to $N\to\infty$, while maintaining the desired
symmetries. Here we formalize the construction of such lattices.

Figure S1:  (a) Shape of a parallelogram unit cell. (b)–(e) The lattices used
with $N=8,12,20,28$. The sites are numbered for the ease of programmatic use.
The dashed circles with the same number are taken as the same site in adjacent
unit cells under PBC.

The infinite square lattice can be covered by repeating parallelogram unit
cells, whose shape is determined by the parameters $a$, $b$, $c$, $d$, and the
number of lattice sites in each unit cell is $N=ac+bd$, as shown in Fig. S1
(a). We denote such a parallelogram unit cell with PBC by the _parallelogram
lattice_.

We want the parallelogram lattices to be as close to square lattices as
possible. The quantities $\delta_{1}=|a^{2}+b^{2}-c^{2}-d^{2}|$ and
$\delta_{2}=|ad-bc|$ are invariant under translations and reflections, where
$\delta_{1}=0$ only if the four edges of the parallelogram have the same
length, and $\delta_{2}=0$ only if the four angles are right. We heuristically
define the _non-squareness_ $\delta=\frac{1}{4}\delta_{1}+\delta_{2}$.

Not all the lattices have the topology to form the Néel and the stripe
orderings. We denote a lattice is _Néel-compatible_ if it is bipartite, and
_stripe-compatible_ if it has exactly two components (connected subgraphs that
are mutually disconnected) with the same size and each component is bipartite.
In corollary, a parallelogram lattice is Néel-compatible only if $a+b$ and
$c+d$ are both even, and it is stripe-compatible only if $a$, $b$, $c$, $d$
are all even.

For each system size $N$, we choose the lattice that is both Néel- and stripe-
compatible, and has minimal $\delta$. The chosen lattices are listed in Table
S1, and the non-trivial shapes are shown in Fig. S1.

$N$ | $(a,c,b,d)$ | $\delta$
---|---|---
8 | (2, 2, 2, 2) | 0
12 | (4, 2, 2, 2) | 7
16 | (4, 4, 0, 0) | 0
20 | (4, 4, 2, 2) | 0
24 | (6, 4, 0, 0) | 5
28 | (6, 4, 2, 2) | 9
32 | (4, 4, 4, 4) | 0
Table S1: The parallelogram lattices used in this paper.

###  A.2 Details of ED computation

The ED is implemented using the software package lattice-symmetries [87],
which is one of the most performant among available software of the kind to
our knowledge. We apply zero total magnetization ($\text{SU}(2)$ symmetry) and
spin inversion symmetry to reduce the Hilbert space, and we have verified that
all random instances have ground states in the even parity sector under spin
inversion. Note that no translational or reflection symmetry can be applied
because of the random plaquettes.

###  A.3 Averaged observables over random instances

We use stratified sampling to generate random instances of the Hamiltonian
with various numbers of $J_{2}$-dominated plaquettes, then reweight them
according to the binomial distribution when computing the averaged observables
given the probability $x$. Therefore, we can reuse the observable values of
the instances for all values of $x$, without recomputing the ED and the
expectations.

For each $k=0,1,\ldots,N$, we generate
$M_{k}=\min(M_{k\text{max}},\binom{N}{k})$ instances randomly without
replacement, each with exactly $k$ $J_{2}$-dominated plaquettes. Here
$M_{k\text{max}}$ is the maximal number of instances for each $k$ limited by
the computation budget, and the actual number $M_{k}$ is also limited by the
number of combinations $\binom{N}{k}$. We take $M_{k\text{max}}=100$ when
$N\leq 20$, $M_{k\text{max}}=30$ when $N=24,28$, and $M_{k\text{max}}=10$ when
$N=32$.

When computing an averaged observable
$[O]_{J}=[\\!\expectationvalue*{O}]_{J}$, we first take a uniform average over
instances with the same $k$, then average over the binomial distribution of
$k$:

$\displaystyle[O]_{J}(N,x)=\sum_{k}B(k,N,x)\bar{O}_{k}(N),$ (S1)
$\displaystyle B(k,N,x)=\binom{N}{k}x^{k}(1-x)^{N-k},$ (S2)
$\displaystyle\bar{O}_{k}(N)=\frac{1}{M_{k}}\sum_{i}O_{i},$ (S3)

where $i$ runs over all $M_{k}$ instances whose number of $J_{2}$-dominated
plaquettes equals $k$, and $O_{i}=\expectationvalue*{O}_{i}$ is the observable
value of this instance. As $O_{i}$ is independent of $x$, we can reuse it for
all values of $x$. For each $x$, we compute $[O]_{J}(N,x)$ with finite values
of $N$, then extrapolate to $[O]_{J}(N\to\infty,x)$, as in Fig. 2.

The standard error of this weighted average is given by

$\displaystyle
O_{\text{err}}^{2}(N,x)=\frac{1}{M_{\text{eff}}(N,x)}\sum_{k}B(k,N,x)O_{k\text{err}}^{2}(N),$
(S4) $\displaystyle M_{\text{eff}}(N,x)=\frac{1}{\sum_{k}B(k,N,x)^{2}},$ (S5)
$\displaystyle
O_{k\text{err}}^{2}(N)=\left(1-\frac{M_{k}}{\binom{N}{k}}\right)\frac{1}{M_{k}^{2}}\sum_{i}(O_{i}-\bar{O}_{k})^{2},$
(S6)

where $M_{\text{eff}}$ is the effective sample size over the binomial
distribution, and $M_{k}/\binom{N}{k}$ is the finite population correction.
The error bars $O_{\text{err}}$ are shown in Fig. 2 when $N<\infty$ and
$0<x<1$.

###  A.4 Uncertainty of extrapolation

We use linear regression to extrapolate observables to $N\to\infty$, and
propagate uncertainties of both the averages over random instances and the
regression.

Consider a linear regression with $P$ data points
$\bm{x}=(x_{1},\ldots,x_{P})^{\top}$, $\bm{y}=(y_{1},\ldots,y_{P})^{\top}$.
The standard errors in $\bm{y}$ are
$\bm{y}_{\text{err}}=(y_{\text{err}1},\ldots,y_{\text{err}P})^{\top}$, and
there is no error in $\bm{x}$. We have

$\bm{y}=\bm{X}\bm{\beta}+\bm{\epsilon},$ (S7)

where $\bm{X}=(\bm{1},\bm{x})$ is a $P\times 2$ matrix including the constant
term, $\bm{\beta}=(\beta_{0},\beta_{1})^{\top}$ contains the intercept and the
slope, and $\bm{\epsilon}$ is the errors of regression. Note that
$\bm{\epsilon}$ is heteroscedastic and we cannot assume that all components of
$\bm{\epsilon}$ are from distributions with the same variance, because the
deviation of $\bm{y}$ from the linear relation is caused by the finite size
effect, whose strength changes with the system size.

The parameters are estimated by least squares:

$\bm{\beta}=\bar{\bm{X}}\bm{y},\quad\bar{\bm{X}}=(\bm{X}^{\top}\bm{X})^{-1}\bm{X}^{\top},$
(S8)

where $\bar{\bm{X}}$ is the pseudoinverse of $\bm{X}$. The variance matrix of
parameters from the errors in $\bm{y}$ is

$\operatorname{Var}^{(1)}\bm{\beta}=\bar{\bm{X}}(\operatorname{diag}\bm{y}_{\text{err}})^{2}\bar{\bm{X}}^{\top},$
(S9)

and another contribution from the errors of regression is

$\displaystyle\operatorname{Var}^{(2)}\bm{\beta}=\frac{P}{P-2}\bar{\bm{X}}(\operatorname{diag}\bm{\epsilon})^{2}\bar{\bm{X}}^{\top},$
(S10) $\displaystyle\bm{\epsilon}=\bm{y}-\bm{X}\bm{\beta}.$ (S11)

The total variance
$\operatorname{Var}\bm{\beta}=\operatorname{Var}^{(1)}\bm{\beta}+\operatorname{Var}^{(2)}\bm{\beta}$,
and the standard error of the intercept $\beta_{0\text{err}}$ is the square
root of its first component. The error bars $\beta_{0\text{err}}$ are shown in
Fig. 2 when $N\to\infty$.

##  Appendix B QFI density from dynamical spin structure factor

Here, we briefly review the approach [32, 36, 37] of obtaining the QFI density
$f_{\text{Q}}$ through the dynamical spin structure in Eq. 17. It assumes
thermal ensembles at the equilibrium and adapts the fluctuation-dissipation
theorem for non-symmetrized $S_{\text{dyn}}(\bm{k},\omega,T)$ [88].

Let us take a generator
$O=\sqrt{N}S^{z\dagger}_{\bm{k}}=\sum_{i=1}^{N}S_{i}^{z}e^{i\bm{k}\cdot\bm{r}_{i}}$.
Considering the definition of QFI density $f_{\text{Q}}=f_{\text{Q}}/N$ and
the form of QFI in Eq. 16, we have

$f_{\text{Q}}(\bm{k},T)=2\sum_{\lambda,\lambda^{\prime}}\frac{(p_{\lambda}-p_{\lambda^{\prime}})^{2}}{p_{\lambda}+p_{\lambda^{\prime}}}\left|\\!\matrixelement{\lambda}{S^{z\dagger}_{\bm{k}}}{\lambda^{\prime}}\right|^{2}.$
(S12)

It is convenient to introduce the dynamical susceptibility

$\chi(\bm{k},\omega,T)=i\int_{0}^{\infty}\Tr\left[\rho[S^{z\dagger}_{\bm{k}}(t),S^{z}_{\bm{k}}(0)]\right]e^{i\omega
t}\differential t,$ (S13)

which, after taking the integral, shares the form

$\displaystyle\chi(\bm{k},\omega,T)=i\pi\sum_{\lambda,\lambda^{\prime}}$
$\displaystyle(p_{\lambda}-p_{\lambda^{\prime}})\left|\\!\matrixelement{\lambda}{S^{z\dagger}_{\bm{k}}(0)}{\lambda^{\prime}}\right|^{2}$
$\displaystyle\times\delta(\omega+E_{\lambda}-E_{\lambda^{\prime}}).$ (S14)

In reaching the form above, we have adopted the density matrix
$\rho=\sum_{\lambda}p_{\lambda}\ket{\lambda}\bra{\lambda}$ with
$p_{\lambda}=\frac{1}{Z}e^{-E_{\lambda}/T}$, and inserted the identity
$\mathds{1}=\sum_{\lambda^{\prime}}\ket{\lambda^{\prime}}\bra{\lambda^{\prime}}$.
Comparing Eq. S14 with Eq. S12 and taking into account

$\frac{p_{\lambda}-p_{\lambda^{\prime}}}{p_{\lambda}+p_{\lambda^{\prime}}}=\tanh(\frac{E_{\lambda^{\prime}}-E_{\lambda}}{2T}),$
(S15)

we arrive at

$f_{\text{Q}}(\bm{k},T)=\frac{4}{\pi}\int_{0}^{\infty}\tanh(\frac{\omega}{2T})\chi^{\prime\prime}(\omega,\bm{k},T)\,\differential\omega,$
(S16)

where $\chi^{\prime\prime}$ denotes the imaginary (dissipative) part of the
dynamical susceptibility.

On the other hand, one can also relate $\chi^{\prime\prime}$ to the non-
symmetrized dynamical spin structure in Eq. 14 through the fluctuation-
dissipation theorem. Starting from the integral definition of $\chi$ in Eq.
S13, it is straightforward to identify that

$\chi^{\prime\prime}(\bm{k},\omega,T)=\frac{1}{2}\left(g_{1}(\bm{k},\omega,T)-g_{2}(\bm{k},\omega,T)\right),$
(S17)

with

$\displaystyle g_{1}(\bm{k},\omega,T)$
$\displaystyle=\int_{-\infty}^{\infty}\\!\expectationvalue*{S^{z\dagger}_{\bm{k}}(t)S^{z}_{\bm{k}}(0)}e^{i\omega
t}\,\differential t,$ $\displaystyle g_{2}(\bm{k},\omega,T)$
$\displaystyle=\int_{-\infty}^{\infty}\\!\expectationvalue*{S^{z}_{\bm{k}}(0)S^{z\dagger}_{\bm{k}}(t)}e^{i\omega
t}\,\differential t.$ (S18)

The notation $\expectationvalue*{\cdots}=\Tr[\rho\cdots]$ is adopted for a
thermal state. Since $g_{1}=2\pi S_{\text{dyn}}$, if we can further associate
$g_{2}$ with $g_{1}$, our goal will be accomplished. It is again more
convenient to compare these two quantities after integration. It turns out
that

$\displaystyle g_{2}(\bm{k},\omega,T)$
$\displaystyle=\sum_{\lambda,\lambda^{\prime}}\frac{e^{-E_{\lambda^{\prime}}/T}}{Z}\left|\\!\matrixelement{\lambda}{S^{z\dagger}_{\bm{k}}(0)}{\lambda^{\prime}}\right|^{2}$
$\displaystyle\phantom{{}=\sum_{\lambda,\lambda^{\prime}}}\times
2\pi\delta(\omega+E_{\lambda}-E_{\lambda^{\prime}})$
$\displaystyle=e^{-\frac{\omega}{T}}g_{1}(\bm{k},\omega,T).$ (S19)

Thus we obtain the adapted fluctuation-dissipation theorem,

$\chi^{\prime\prime}(\bm{k},\omega,T)=\pi\left(1-e^{-\frac{\omega}{T}}\right)S_{\text{dyn}}(\bm{k},\omega,T).$
(S20)

Combining Eqs. S16 and S20, the entanglement witness $f_{\text{Q}}$ is
reconstructed from the quantum fluctuations contained in $S_{\text{dyn}}$, and
we reach Eq. 17:

$f_{\text{Q}}(\bm{k},T)=4\int_{0}^{\infty}\\!\\!\tanh(\frac{\omega}{2T})\left(1-e^{-\frac{\omega}{T}}\right)S_{\text{dyn}}(\bm{k},\omega,T)\,\differential\omega.$
(S21)

## References

  * Anderson [1987] P. W. Anderson, The resonating valence bond state in La2CuO4 and superconductivity, Science 235, 1196 (1987).
  * Kivelson  _et al._ [1987] S. A. Kivelson, D. S. Rokhsar, and J. P. Sethna, Topology of the resonating valence-bond state: Solitons and high-$T_{c}$ superconductivity, Phys. Rev. B 35, 8865 (1987).
  * Kalmeyer and Laughlin [1987] V. Kalmeyer and R. B. Laughlin, Equivalence of the resonating-valence-bond and fractional quantum Hall states, Phys. Rev. Lett. 59, 2095 (1987).
  * Savary and Balents [2016] L. Savary and L. Balents, Quantum spin liquids: A review, Rep. Prog. Phys. 80, 016502 (2016).
  * Zhou  _et al._ [2017] Y. Zhou, K. Kanoda, and T.-K. Ng, Quantum spin liquid states, Rev. Mod. Phys. 89, 025003 (2017).
  * Wen [2017] X.-G. Wen, Colloquium: Zoo of quantum-topological phases of matter, Rev. Mod. Phys. 89, 041004 (2017).
  * Knolle and Moessner [2019] J. Knolle and R. Moessner, A field guide to spin liquids, Annu. Rev. Condens. Matter Phys. 10, 451 (2019).
  * Wen  _et al._ [2019] J. Wen, S.-L. Yu, S. Li, W. Yu, and J.-X. Li, Experimental identification of quantum spin liquids, npj Quantum Mater. 4, 12 (2019).
  * Broholm  _et al._ [2020] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, Quantum spin liquids, Science 367, eaay0668 (2020).
  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Ann. Phys. 321, 2 (2006).
  * Feng  _et al._ [2007] X.-Y. Feng, G.-M. Zhang, and T. Xiang, Topological characterization of quantum phase transitions in a spin-$1/2$ model, Phys. Rev. Lett. 98, 087204 (2007).
  * Le Hur  _et al._ [2017] K. Le Hur, A. Soret, and F. Yang, Majorana spin liquids, topology, and superconductivity in ladders, Phys. Rev. B 96, 205109 (2017).
  * Yao and Kivelson [2007] H. Yao and S. A. Kivelson, Exact chiral spin liquid with non-abelian anyons, Phys. Rev. Lett. 99, 247203 (2007).
  * O’Brien  _et al._ [2016] K. O’Brien, M. Hermanns, and S. Trebst, Classification of gapless $\mathbb{Z}_{2}$ spin liquids in three-dimensional Kitaev models, Phys. Rev. B 93, 085101 (2016).
  * Ma [2023] H. Ma, $\mathbb{Z}_{2}$ spin liquids in the higher spin-$S$ Kitaev honeycomb model: An exact deconfined $\mathbb{Z}_{2}$ gauge structure in a nonintegrable model, Phys. Rev. Lett. 130, 156701 (2023).
  * Wilczek [1982] F. Wilczek, Quantum mechanics of fractional-spin particles, Phys. Rev. Lett. 49, 957 (1982).
  * Majorana [1937] E. Majorana, Teoria simmetrica dell’elettrone e del positrone, Nuovo Cim. 14, 171 (1937).
  * Alicea [2012] J. Alicea, New directions in the pursuit of Majorana fermions in solid state systems, Rep. Prog. Phys. 75, 076501 (2012).
  * Xu and Kee [2024] H. Xu and H.-Y. Kee, Manipulating Majorana zero modes in the spin-1/2 Kitaev ladder, (2024), arXiv:2404.16099 .
  * Kitaev [2003] A. Y. Kitaev, Fault-tolerant quantum computation by anyons, Ann. Phys. 303, 2 (2003).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. D. Sarma, Non-abelian anyons and topological quantum computation, Rev. Mod. Phys. 80, 1083 (2008).
  * Duan  _et al._ [2003] L.-M. Duan, E. Demler, and M. D. Lukin, Controlling spin exchange interactions of ultracold atoms in optical lattices, Phys. Rev. Lett. 91, 090402 (2003).
  * Yang  _et al._ [2018] F. Yang, L. Henriet, A. Soret, and K. Le Hur, Engineering quantum spin liquids and many-body Majorana states with a driven superconducting box circuit, Phys. Rev. B 98, 035431 (2018).
  * Sagi  _et al._ [2019] E. Sagi, H. Ebisu, Y. Tanaka, A. Stern, and Y. Oreg, Spin liquids from Majorana zero modes in a cooper-pair box, Phys. Rev. B 99, 075107 (2019).
  * Semeghini  _et al._ [2021] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, _et al._ , Probing topological spin liquids on a programmable quantum simulator, Science 374, 1242 (2021).
  * Sun  _et al._ [2023] B.-Y. Sun, N. Goldman, M. Aidelsburger, and M. Bukov, Engineering and probing non-abelian chiral spin liquids using periodically driven ultracold atoms, PRX Quantum 4, 020329 (2023).
  * Cookmeyer and Sarma [2024] T. Cookmeyer and S. D. Sarma, Engineering the Kitaev spin liquid in a quantum dot system, Phys. Rev. Lett. 132, 186501 (2024).
  * Balents [2010] L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).
  * Norman [2016] M. R. Norman, Colloquium: Herbertsmithite and the search for the quantum spin liquid, Rev. Mod. Phys. 88, 041002 (2016).
  * Trebst and Hickey [2022] S. Trebst and C. Hickey, Kitaev materials, Phys. Rep. 950, 1 (2022).
  * Song  _et al._ [2012] H. F. Song, S. Rachel, C. Flindt, I. Klich, N. Laflorencie, and K. Le Hur, Bipartite fluctuations as a probe of many-body entanglement, Phys. Rev. B 85, 035409 (2012).
  * Hauke  _et al._ [2016] P. Hauke, M. Heyl, L. Tagliacozzo, and P. Zoller, Measuring multipartite entanglement through dynamic susceptibilities, Nat. Phys. 12, 778 (2016).
  * Laflorencie [2016] N. Laflorencie, Quantum entanglement in condensed matter systems, Phys. Rep. 646, 1 (2016).
  * Irfan  _et al._ [2021] A. A. M. Irfan, P. Blackstone, R. Pynn, and G. Ortiz, Quantum entangled-probe scattering theory, New J. Phys. 23, 083022 (2021).
  * Yang  _et al._ [2020] F. Yang, K. Plekhanov, and K. Le Hur, Valence bond fluctuations in the Kitaev spin model, Phys. Rev. Res. 2, 013005 (2020).
  * Scheie  _et al._ [2021] A. Scheie, P. Laurell, A. M. Samarakoon, B. Lake, S. E. Nagler, G. E. Granroth, S. Okamoto, G. Alvarez, and D. A. Tennant, Witnessing entanglement in quantum magnets using neutron scattering, Phys. Rev. B 103, 224434 (2021).
  * Laurell  _et al._ [2021] P. Laurell, A. Scheie, C. J. Mukherjee, M. M. Koza, M. Enderle, Z. Tylczynski, S. Okamoto, R. Coldea, D. A. Tennant, and G. Alvarez, Quantifying and controlling entanglement in the quantum magnet Cs2CoCl4, Phys. Rev. Lett. 127, 037201 (2021).
  * Fang  _et al._ [2024] Y. Fang, M. Mahankali, Y. Wang, L. Chen, H. Hu, S. Paschen, and Q. Si, Amplified entanglement witnessed in a quantum critical metal, (2024), arXiv:2402.18552 .
  * Ren  _et al._ [2024] T. Ren, Y. Shen, S. F. R. TenHuisen, J. Sears, W. He, M. H. Upton, D. Casa, P. Becker, M. Mitrano, M. P. M. Dean, _et al._ , Witnessing quantum entanglement using resonant inelastic X-ray scattering, (2024), arXiv:2404.05850 .
  * Leiner  _et al._ [2024] J. Leiner, S. Kuhn, S. McKay, J. Jochum, F. Li, A. A. M. Irfan, F. Funama, D. Mettus, L. Beddrich, C. Franz, _et al._ , Spin-energy entanglement of a time-focused neutron, (2024), arXiv:2404.07967 .
  * Preskill [2000] J. Preskill, Quantum information and physics: Some future directions, J. Mod. Opt. 47, 127 (2000).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, New York, 2010).
  * Katukuri  _et al._ [2020] V. M. Katukuri, P. Babkevich, O. Mustonen, H. C. Walker, B. Fåk, S. Vasala, M. Karppinen, H. M. Rønnow, and O. V. Yazyev, Exchange interactions mediated by nonmagnetic cations in double perovskites, Phys. Rev. Lett. 124, 077202 (2020).
  * Mustonen  _et al._ [2018a] O. Mustonen, S. Vasala, E. Sadrollahi, K. P. Schmidt, C. Baines, H. C. Walker, I. Terasaki, F. J. Litterst, E. Baggio-Saitovitch, and M. Karppinen, Spin-liquid-like state in a spin-$1/2$ square-lattice antiferromagnet perovskite induced by $d^{10}$-$d^{0}$ cation mixing, Nat. Commun. 9, 1085 (2018a).
  * Mustonen  _et al._ [2018b] O. Mustonen, S. Vasala, K. P. Schmidt, E. Sadrollahi, H. C. Walker, I. Terasaki, F. J. Litterst, E. Baggio-Saitovitch, and M. Karppinen, Tuning the $S=1/2$ square-lattice antiferromagnet Sr2Cu(Te1-xWx)O6 from Néel order to quantum disorder to columnar order, Phys. Rev. B 98, 064411 (2018b).
  * Watanabe  _et al._ [2018] M. Watanabe, N. Kurita, H. Tanaka, W. Ueno, K. Matsui, and T. Goto, Valence-bond-glass state with a singlet gap in the spin-$\frac{1}{2}$ square-lattice random $J_{1}$-$J_{2}$ Heisenberg antiferromagnet Sr2CuTe1-xWxO6, Phys. Rev. B 98, 054422 (2018).
  * Vasala  _et al._ [2014] S. Vasala, M. Avdeev, S. Danilkin, O. Chmaissem, and M. Karppinen, Magnetic structure of Sr2CuWO6, J. Condens. Matter Phys. 26, 496001 (2014).
  * Babkevich  _et al._ [2016] P. Babkevich, V. M. Katukuri, B. Fåk, S. Rols, T. Fennell, D. Pajić, H. Tanaka, T. Pardini, R. R. P. Singh, A. Mitrushchenkov, _et al._ , Magnetic excitations and electronic interactions in Sr2CuTeO6: A spin-$1/2$ square lattice Heisenberg antiferromagnet, Phys. Rev. Lett. 117, 237203 (2016).
  * Walker  _et al._ [2016] H. C. Walker, O. Mustonen, S. Vasala, D. J. Voneshen, M. D. Le, D. T. Adroja, and M. Karppinen, Spin wave excitations in the tetragonal double perovskite Sr2CuWO6, Phys. Rev. B 94, 064411 (2016).
  * Koga  _et al._ [2016] T. Koga, N. Kurita, M. Avdeev, S. Danilkin, T. J. Sato, and H. Tanaka, Magnetic structure of the $S=1/2$ quasi-two-dimensional square-lattice Heisenberg antiferromagnet Sr2CuTeO6, Phys. Rev. B 93, 054426 (2016).
  * Hong  _et al._ [2021] W. Hong, L. Liu, C. Liu, X. Ma, A. Koda, X. Li, J. Song, W. Yang, J. Yang, P. Cheng, H. Zhang, W. Bao, X. Ma, D. Chen, K. Sun, W. Guo, H. Luo, A. W. Sandvik, and S. Li, Extreme suppression of antiferromagnetic order and critical scaling in a two-dimensional random quantum magnet, Phys. Rev. Lett. 126, 037201 (2021).
  * Fogh  _et al._ [2022] E. Fogh, O. Mustonen, P. Babkevich, V. M. Katukuri, H. C. Walker, L. Mangin-Thro, M. Karppinen, S. Ward, B. Normand, and H. M. Rønnow, Randomness and frustration in a $S=\frac{1}{2}$ square-lattice Heisenberg antiferromagnet, Phys. Rev. B 105, 184410 (2022).
  * Uematsu and Kawamura [2018] K. Uematsu and H. Kawamura, Randomness-induced quantum spin liquid behavior in the $s=\frac{1}{2}$ random $J_{1}$-$J_{2}$ Heisenberg antiferromagnet on the square lattice, Phys. Rev. B 98, 134427 (2018).
  * Ren  _et al._ [2023] H.-D. Ren, T.-Y. Xiong, H.-Q. Wu, D. N. Sheng, and S.-S. Gong, Characterizing random-singlet state in two-dimensional frustrated quantum magnets and implications for the double perovskite Sr2CuTe1-xWxO6, Phys. Rev. B 107, L020407 (2023).
  * Liu  _et al._ [2018] L. Liu, H. Shao, Y.-C. Lin, W. Guo, and A. W. Sandvik, Random-singlet phase in disordered two-dimensional quantum magnets, Phys. Rev. X 8, 041040 (2018).
  * Bhatt and Lee [1982] R. N. Bhatt and P. A. Lee, Scaling studies of highly disordered spin-$1/2$ antiferromagnetic systems, Phys. Rev. Lett. 48, 344 (1982).
  * Fisher [1994] D. S. Fisher, Random antiferromagnetic quantum spin chains, Phys. Rev. B 50, 3799 (1994).
  * Baek  _et al._ [2020] S.-H. Baek, H. W. Yeo, S.-H. Do, K.-Y. Choi, L. Janssen, M. Vojta, and B. Büchner, Observation of a random singlet state in a diluted Kitaev honeycomb material, Phys. Rev. B 102, 094407 (2020).
  * Tarzia and Biroli [2008] M. Tarzia and G. Biroli, The valence bond glass phase, Europhys. Lett. 82, 67008 (2008).
  * Singh [2010] R. R. P. Singh, Valence bond glass phase in dilute kagome antiferromagnets, Phys. Rev. Lett. 104, 177203 (2010).
  * Takahashi [1987] M. Takahashi, Classical Heisenberg ferromagnet in two dimensions, Phys. Rev. B 36, 3791 (1987).
  * Song  _et al._ [2011] H. F. Song, N. Laflorencie, S. Rachel, and K. Le Hur, Entanglement entropy of the two-dimensional Heisenberg antiferromagnet, Phys. Rev. B 83, 224410 (2011).
  * Edwards and Anderson [1975] S. F. Edwards and P. W. Anderson, Theory of spin glasses, J. Phys. F: Met. Phys. 5, 965 (1975).
  * Binder and Young [1986] K. Binder and A. P. Young, Spin glasses: Experimental facts, theoretical concepts, and open questions, Rev. Mod. Phys. 58, 801 (1986).
  * Wu  _et al._ [2019] H.-Q. Wu, S.-S. Gong, and D. N. Sheng, Randomness-induced spin-liquid-like phase in the spin-$1/2$ $J_{1}$-$J_{2}$ triangular Heisenberg model, Phys. Rev. B 99, 085141 (2019).
  * Nomura and Imada [2021] Y. Nomura and M. Imada, Dirac-type nodal spin liquid revealed by refined quantum many-body solver using neural-network wave function, correlation ratio, and level spectroscopy, Phys. Rev. X 11, 031034 (2021).
  * Saha  _et al._ [2023] B. Saha, A. K. Bera, S. M. Yusuf, and A. Hoser, Two-dimensional short-range spin-spin correlations in the layered spin-$\frac{3}{2}$ maple leaf lattice antiferromagnet Na2Mn3O7 with crystal stacking disorder, Phys. Rev. B 107, 064419 (2023).
  * Carleo  _et al._ [2012] G. Carleo, F. Becca, M. Schiró, and M. Fabrizio, Localization and glassy dynamics of many-body quantum systems, Sci. Rep. 2, 243 (2012).
  * Gagliano and Balseiro [1987] E. R. Gagliano and C. A. Balseiro, Dynamical properties of quantum many-body systems at zero temperature, Phys. Rev. Lett. 59, 2999 (1987).
  * Braunstein and Caves [1994] S. L. Braunstein and C. M. Caves, Statistical distance and the geometry of quantum states, Phys. Rev. Lett. 72, 3439 (1994).
  * Yu  _et al._ [2022] M. Yu, Y. Liu, P. Yang, M. Gong, Q. Cao, S. Zhang, H. Liu, M. Heyl, T. Ozawa, N. Goldman, _et al._ , Quantum Fisher information measurement and verification of the quantum Cramér–Rao bound in a solid-state qubit, npj Quantum Inf. 8, 56 (2022).
  * Hyllus  _et al._ [2012] P. Hyllus, W. Laskowski, R. Krischek, C. Schwemmer, W. Wieczorek, H. Weinfurter, L. Pezzé, and A. Smerzi, Fisher information and multiparticle entanglement, Phys. Rev. A 85, 022321 (2012).
  * Tóth [2012] G. Tóth, Multipartite entanglement and high-precision metrology, Phys. Rev. A 85, 022322 (2012).
  * Pezzé and Smerzi [2014] L. Pezzé and A. Smerzi, Quantum theory of phase estimation, in  _Atom Interferometry, Proceedings of the International School of Physics ‘Enrico Fermi’, Course 188_ (IOS Press, Varenna, 2014) pp. 691–741.
  * Silverman [1986] B. W. Silverman, Density estimation for statistics and data analysis, in  _Monographs on Statistics and Applied Probability_ , Vol. 26 (Chapman and Hall, 1986).
  * Aref  _et al._ [2020] S. Aref, A. J. Mason, and M. C. Wilson, A modeling and computational study of the frustration index in signed networks, Networks 75, 95 (2020).
  * Hams and de Raedt [2000] A. Hams and H. de Raedt, Fast algorithm for finding the eigenvalue distribution of very large matrices, Phys. Rev. E 62, 4365 (2000).
  * Baykusheva  _et al._ [2023] D. R. Baykusheva, M. H. Kalthoff, D. Hofmann, M. Claassen, D. M. Kennes, M. A. Sentef, and M. Mitrano, Witnessing nonequilibrium entanglement dynamics in a strongly correlated fermionic chain, Phys. Rev. Lett. 130, 106902 (2023).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Carleo  _et al._ [2019] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Machine learning and the physical sciences, Rev. Mod. Phys. 91, 045002 (2019).
  * Viteritti  _et al._ [2023a] L. L. Viteritti, R. Rende, and F. Becca, Transformer variational wave functions for frustrated quantum spin systems, Phys. Rev. Lett. 130, 236401 (2023a).
  * Viteritti  _et al._ [2023b] L. L. Viteritti, R. Rende, A. Parola, S. Goldt, and F. Becca, Transformer wave function for the Shastry-Sutherland model: Emergence of a spin-liquid phase, (2023b), arXiv:2311.16889 .
  * Roth  _et al._ [2023] C. Roth, A. Szabó, and A. H. MacDonald, High-accuracy variational Monte Carlo for frustrated magnets with deep neural networks, Phys. Rev. B 108, 054410 (2023).
  * Chen and Heyl [2023] A. Chen and M. Heyl, Efficient optimization of deep neural quantum states toward machine precision, (2023), arXiv:2302.01941 .
  * Rende  _et al._ [2023] R. Rende, L. L. Viteritti, L. Bardone, F. Becca, and S. Goldt, A simple linear algebra identity to optimize large-scale neural network quantum states, (2023), arXiv:2310.05715 .
  * Zhang  _et al._ [2024] Y. Zhang, X.-Y. Song, and T. Senthil, Dirac spin liquid as an “unnecessary” quantum critical point on square lattice antiferromagnets, (2024), arXiv:2404.11654 .
  * Westerhout [2021] T. Westerhout, lattice-symmetries: A package for working with quantum many-body bases, J. Open Source Software 16, 3537 (2021).
  * Lovesey [1984] S. W. Lovesey,  _Theory of Neutron Scattering From Condensed Matter_ (Clarendon Press, Oxford, UK, 1984).


* [77] J. Kubo, A. Mondragon, M. Mondragon and E. Rodriguez-Jauregui, _The Flavor symmetry_ , _Prog. Theor. Phys._ 109 (2003) 795–807, [hep-ph/0302196].
  * [78] T. Teshima, _Flavor mass and mixing and S(3) symmetry: An S(3) invariant model reasonable to all_ , _Phys. Rev. D_ 73 (2006) 045019, [hep-ph/0509094].
  * [79] A. Mondragon, M. Mondragon and E. Peinado, _S(3)-flavour symmetry as realized in lepton flavour violating processes_ , _J. Phys. A_ 41 (2008) 304035, [0712.1799].
  * [80] A. Mondragon, M. Mondragon and E. Peinado, _Lepton masses, mixings and FCNC in a minimal S(3)-invariant extension of the Standard Model_ , _Phys. Rev. D_ 76 (2007) 076003, [0706.0354].
  * [81] F. Gonzalez Canales, A. Mondragon and M. Mondragon, _The $S_{3}$ Flavour Symmetry: Neutrino Masses and Mixings_, _Fortsch. Phys._ 61 (2013) 546–570, [1205.4755].
  * [82] E. Ma and B. Melic, _Updated $S_{3}$ model of quarks_, _Phys. Lett. B_ 725 (2013) 402–406, [1303.6928].
  * [83] F. González Canales, A. Mondragón, M. Mondragón, U. J. Saldaña Salazar and L. Velasco-Sevilla, _Quark sector of S3 models: classification and comparison with experimental data_ , _Phys. Rev. D_ 88 (2013) 096004, [1304.6644].
  * [84] D. Das, U. K. Dey and P. B. Pal, _$S_{3}$ symmetry and the quark mixing matrix_, _Phys. Lett. B_ 753 (2016) 315–318, [1507.06509].
  * [85] A. A. Cruz and M. Mondragón, _Neutrino masses, mixing, and leptogenesis in an S3 model_ , 1701.07929.
  * [86] Particle Data Group collaboration, R. L. Workman et al., _Review of Particle Physics_ , _PTEP_ 2022 (2022) 083C01.
  * [87] C. Jarlskog, _Commutator of the Quark Mass Matrices in the Standard Electroweak Model and a Measure of Maximal $CP$ Nonconservation_, _Phys. Rev. Lett._ 55 (1985) 1039.
  * [88] J. Bernabeu, G. C. Branco and M. Gronau, _CP Restrictions on Quark Mass Matrices_ , _Phys. Lett. B_ 169 (1986) 243–247.
  * [89] CMS collaboration, A. M. Sirunyan et al., _Inclusive search for highly boosted Higgs bosons decaying to bottom quark-antiquark pairs in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _JHEP_ 12 (2020) 085, [2006.13251].
  * [90] ATLAS collaboration, G. Aad et al., _Measurements of $WH$ and $ZH$ production in the $H\rightarrow b\bar{b}$ decay channel in $pp$ collisions at 13 TeV with the ATLAS detector_, _Eur. Phys. J. C_ 81 (2021) 178, [2007.02873].
  * [91] ATLAS collaboration, G. Aad et al., _Measurement of the associated production of a Higgs boson decaying into $b$-quarks with a vector boson at high transverse momentum in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, _Phys. Lett. B_ 816 (2021) 136204, [2008.02508].
  * [92] CMS collaboration, A. Tumasyan et al., _Search for flavor-changing neutral current interactions of the top quark and the Higgs boson decaying to a bottom quark-antiquark pair at $\sqrt{s}$ = 13 TeV_, _JHEP_ 02 (2022) 169, [2112.09734].
  * [93] ATLAS collaboration, _Search for flavour-changing neutral current interactions of the top quark and the Higgs boson in events with a pair of $\tau$-leptons in pp collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, 2208.11415.
  * [94] ATLAS collaboration, G. Aad et al., _$CP$ Properties of Higgs Boson Interactions with Top Quarks in the $t\bar{t}H$ and $tH$ Processes Using $H\rightarrow\gamma\gamma$ with the ATLAS Detector_, _Phys. Rev. Lett._ 125 (2020) 061802, [2004.04545].
  * [95] CMS collaboration, _Search for $CP$ violation in ttH and tH production in multilepton channels in proton-proton collisions at $\sqrt{s}$ = 13 TeV_, 2208.02686.
  * [96] CMS collaboration, A. M. Sirunyan et al., _Search for a light charged Higgs boson in the H ± $\to$ cs channel in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _Phys. Rev. D_ 102 (2020) 072001, [2005.08900].
  * [97] ATLAS collaboration, _Search for a light charged Higgs boson in t $\rightarrow$H+b decays, with H+$\rightarrow$cb, in the lepton+jets final state in proton-proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, .
  * [98] A. Pierce and J. Thaler, _Natural Dark Matter from an Unnatural Higgs Boson and New Colored Particles at the TeV Scale_ , _JHEP_ 08 (2007) 026, [hep-ph/0703056].
  * [99] A. Arbey, F. Mahmoudi, O. Stål and T. Stefaniak, _Status of the Charged Higgs Boson in Two Higgs Doublet Models_ , _Eur. Phys. J. C_ 78 (2018) 182, [1706.07414].
  * [100] ALEPH, DELPHI, L3, OPAL, LEP Electroweak collaboration, S. Schael et al., _Electroweak Measurements in Electron-Positron Collisions at W-Boson-Pair Energies at LEP_ , _Phys. Rept._ 532 (2013) 119–244, [1302.3415].
  * [101] M. E. Peskin and T. Takeuchi, _A New constraint on a strongly interacting Higgs sector_ , _Phys. Rev. Lett._ 65 (1990) 964–967.
  * [102] M. E. Peskin and T. Takeuchi, _Estimation of oblique electroweak corrections_ , _Phys. Rev. D_ 46 (1992) 381–409.
  * [103] W. Grimus, L. Lavoura, O. M. Ogreid and P. Osland, _A Precision constraint on multi-Higgs-doublet models_ , _J. Phys._ G35 (2008) 075001, [0711.4022].
  * [104] W. Grimus, L. Lavoura, O. M. Ogreid and P. Osland, _The Oblique parameters in multi-Higgs-doublet models_ , _Nucl. Phys._ B801 (2008) 81–96, [0802.4353].
  * [105] P. J. Fox and N. Weiner, _Light Signals from a Lighter Higgs_ , _JHEP_ 08 (2018) 025, [1710.07649].
  * [106] M. Cepeda, S. Gori, V. M. Outschoorn and J. Shelton, _Exotic Higgs Decays_ , 2111.12751.
  * [107] T. Biekötter, S. Heinemeyer and G. Weiglein, _Mounting evidence for a 95 GeV Higgs boson_ , _JHEP_ 08 (2022) 201, [2203.13180].
  * [108] T. Robens, _A Short Overview on Low Mass Scalars at Future Lepton Colliders_ , _Universe_ 8 (2022) 286, [2205.09687].
  * [109] R. Plantey, O. M. Ogreid, P. Osland, M. N. Rebelo and M. A. Solberg, _Weinberg’s 3HDM potential with spontaneous CP violation_ , 2208.13594.
  * [110] A. Kunčinas, O. M. Ogreid, P. Osland and M. N. Rebelo, _Revisiting two dark matter candidates in $S_{3}$-symmetric three-Higgs-doublet models_, in _8th Symposium on Prospects in the Physics of Discrete Symmetries_ , 1, 2023.  2301.12194.
  * [111] G. C. Branco and I. P. Ivanov, _Group-theoretic restrictions on generation of CP-violation in multi-Higgs-doublet models_ , _JHEP_ 01 (2016) 116, [1511.02764].
  * [112] B. W. Lee, C. Quigg and H. Thacker, _Weak Interactions at Very High-Energies: The Role of the Higgs Boson Mass_ , _Phys. Rev. D_ 16 (1977) 1519.
  * [113] M. Luscher and P. Weisz, _Is There a Strong Interaction Sector in the Standard Lattice Higgs Model?_ , _Phys. Lett._ B212 (1988) 472–478.
  * [114] W. J. Marciano, G. Valencia and S. Willenbrock, _Renormalization Group Improved Unitarity Bounds on the Higgs Boson and Top Quark Masses_ , _Phys. Rev._ D40 (1989) 1725.
  * [115] J. M. Cornwall, D. N. Levin and G. Tiktopoulos, _Derivation of Gauge Invariance from High-Energy Unitarity Bounds on the s Matrix_ , _Phys. Rev._ D10 (1974) 1145.
  * [116] C. E. Vayonakis, _Born Helicity Amplitudes and Cross-Sections in Nonabelian Gauge Theories_ , _Lett. Nuovo Cim._ 17 (1976) 383.