Density-Matrix Mean-Field Theory

Junyi Zhang1, 2 $\star$ and Zhengqian Cheng3

1 William H. Miller III Department of Physics and Astronomy, Johns Hopkins
University, Baltimore, Maryland 21218, USA

2 Institute for Quantum Matter, Johns Hopkins University, Baltimore, Maryland,
21218, USA

3 Department of Applied Physics and Applied Mathematics, Columbia University,
New York, New York 10027, USA

⋆<EMAIL_ADDRESS>
## Abstract

Mean-field theories (MFTs) have proven to be efficient tools for exploring
various phases of matter, complementing alternative methods that are more
precise but also more computationally demanding. Conventional mean-field
theories (MFTs) often fall short in capturing quantum fluctuations, which
restricts their applicability to systems characterized by strong quantum
fluctuations. In this article, we propose a novel mean-field theory, density-
matrix mean-field theory (DMMFT). DMMFT constructs effective Hamiltonians,
incorporating quantum environments shaped by entanglements quantified by the
reduced density matrices. Therefore, it offers a systematic and unbiased
approach to account for effects of fluctuations and entanglements in quantum
ordered phases. As demonstrative examples, we show that DMMFT can not only
quantitatively evaluate the renormalization of order parameters induced by
quantum fluctuations but can even detect the topological order of quantum
phases. Additionally, we discuss the extensions of DMMFT for systems at finite
temperatures and those with disorders. Our work provides a novel and efficient
approach to explore phases exhibiting unconventional quantum orders, which can
be particularly beneficial for investigating frustrated spin systems in high
spatial dimensions.

###### Contents

  1. 1 Introduction
  2. 2 Density-Matrix Mean-Field Theory
    1. 2.1 Conventional Mean-Field Approximation
    2. 2.2 Density-Matrix Mean-Field Approximation
    3. 2.3 Comparison of the Mean-Field Approximations
  3. 3 Applications
    1. 3.1 Affleck-Kennedy-Lieb-Tasaki Model
    2. 3.2 Antiferromagnetic Heisenberg Model on Triangular Lattices
  4. 4 Discussions
    1. 4.1 Comparisons with Alternative Methods
    2. 4.2 Systems at Finite Temperatures
    3. 4.3 Systems with Disorders
  5. 5 Conclusion
  6. A Iterative Algorithm for Mean-Field Equations

##  1 Introduction

Frustrated Hubbard and Heisenberg models [1, 2, 3] have continued to capture
research attention over the last half-century due to their potential to host
various intriguing quantum phases [4, 5, 6, 7, 8], as well as their relevance
to high $T_{c}$ superconductors [9] and their applications in quantum
computations [10]. Determining the ground states of frustrated Hubbard and
Heisenberg models is often a challenging task.

Exact approaches frequently encounter limitations posed by the exponential
wall. In exact diagonalizations (ED), the dimension of the Hilbert space
increases exponentially with the system size. In the density-matrix
renormalization group (DMRG) [11], the area law provides relief from the
exponential wall issue for gapped systems in one spatial dimension ($d=1$), it
still faces exponential scaling challenges for gapless systems or in higher
dimensions ($d\geq 2$)[12, 13, 14]. On the other hand, quantum Monte Carlo
(QMC) methods are less constrained by the system size. Nevertheless, the
notorious sign problem often plagues fermionic systems and frustrated magnetic
systems [15, 16].

Approximation methods, serving as complements to exact approaches, prove to be
useful and efficient tools for exploring various phases of the systems. The
usual mean-field theories (MFTs) already provide insights into non-trivial
effects arising resulting from the interactions, e.g., the formation of the
local moments in metals [17], and the BCS theory for the superconductivity
[18]. Conventional MFTs achieve the simplification by neglecting the
fluctuations, Consequently, they tend to exhibits bias towards ordered states
and overlook the nuanced effects stemming from the fluctuations.

Beyond conventional MFTs, various approximations have been proposed for
fermionic systems [19, 20, 21, 22, 23, 24, 25]. In dynamical mean-field theory
(DMFT), a lattice model is mapped to a local impurity model, and the effective
action is constructed using Green’s functions as dynamical mean fields [19],
which allows DMFT to capture the quantum features of the metal-insulator
transitions [20]. While DMFT precisely describes systems in infinite
dimensions, the computational demands of solving the local impurity problems
with continuous baths necessitate ongoing efforts to further simplify the
DMFT. Recently, a quantum embedding method called density matrix embedding
theory (DMET) has been introduced to enhance the efficiency of the DMFT by
taking the advantage of the frequency-independent local density matrix [21].
More recently, another simplification of DMFT known as variational discrete
action theory (VDAT) has also been proposed utilizing sequential product
density matrix to variationally determine ground states [22, 23, 24, 25].

Despite of the successes of DMFT and its simplifications for fermionic
systems, a gap remains in methods beyond semiclassical mean-field
approximations for spin systems. Although employing the Holstein-Primakoff
transformation [26] allows spins to be mapped to bosons, to which the DMFT may
be adapted in principle, the transformation is nonlinear, and the
semicalssical large-$S$ expansion becomes less controllable for $S=1/2$ in the
quantum limit, posing challenges for semiclassical methods applied to quantum
spins. Similar challenges may be encountered with alternative methods. For
instance, in the Schwinger boson representation [27], one needs to set
$\mathcal{N}=2$ for the quantum limit after a saddle-point mean-field
approximation with a large $\mathcal{N}$ .

In this article, we propose a new mean-field method beyond conventional MFTs,
which we call density-matrix mean-field theory (DMMFT). DMMFT constructs
effective Hamiltonians, incorporating quantum environments shaped by
entanglements quantified by the reduced density matrices without presumed
semiclassical orders. Therefore, it offers an unbiased approach to account for
effects of fluctuations and entanglements in quantum ordered phases. In
contrast to QMC and DMFT, DMMFT is generically applicable to systems of
fermions, bosons as well as spins, regardless of frustrations. More
importantly, by gauging the quantum fluctuations with the reduced DM, DMMFT
can detect not only the symmetry-breaking phases in Landau’s paradigm but also
the topological phases with the help of entanglement spectra. Our work
provides a novel and efficient approach to explore phases exhibiting
unconventional quantum orders, Particularly, it fills the gap left by the MFTs
in studying quantum ordered phases in frustrated spin systems, where
semiclassical methods become less controllable in the quantum limit, and QMC
methods fail due to the sign problem.

The rest of the article is organized as follows. In Sec. 2, we present the
generic formulation of the DMMFT. The mean-field equations of DMMFT are
derived in parallel to those of conventional MFTs. Moreover, we demonstrate
that DMMFT becomes equivalent to conventional MFTs when the hyperparameter,
which gauge the quantum fluctuations, is minimized. In Sec. 3, we apply the
DMMFT to two demonstrative examples, 1) the Affleck-Khomoto-Lieb-Tasaki (AKLT)
model [28, 29, 30], and 2) the antiferromagnetic Heisenberg model on
triangular lattices (AFHTL) [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42,
43, 44, 45, 46, 47]. In the AKLT model, DMFT identifies the topological ground
states through their entanglement spectra. In the AFHTL, DMMFT reveals that
quantum fluctuations not only renormalize the order parameters but also shift
the phase boundaries. In Sec. 4, we compare the DMMFT with the DMRG and the
DMFT. Additionally, we discuss the extensions of DMMFT for systems at finite
temperatures and those with disorders. Finally, in Sec. 5, we draw conclusions
and discuss potential avenues for future research

##  2 Density-Matrix Mean-Field Theory

In this section, we formulate the DMMFT for a generic Hamiltonian. Let
$\mathcal{I}=\\{i\\}$ be a set of all sites. On each site, there is a
collection of local operators $\mathcal{O}_{i}=\\{O_{i}^{\alpha}\\}$. A
generic Hamiltonian of local operators can be organized as follows

$\begin{split}H[\mathcal{O}_{\mathcal{I}}]=\sum_{c\in\mathcal{C}}H_{c}[\mathcal{O}_{c}]+\sum_{(c,c^{\prime})}H_{c,c^{\prime}}[\mathcal{O}_{c},\mathcal{O}_{c^{\prime}}]+\sum_{(c,c^{\prime},c^{\prime\prime})}H_{c,c^{\prime},c^{\prime\prime}}[\mathcal{O}_{c},\mathcal{O}_{c^{\prime}},\mathcal{O}_{c^{\prime\prime}}]+\dots,\end{split}$
(1)

where $\mathcal{C}=\\{c\\}$ is a partition of $\mathcal{I}$, and
$\mathcal{O}_{\mathcal{S}}=\cup_{i\in\mathcal{S}}\mathcal{O}_{i}$ is the
collection of local operators over the sites within a set $\mathcal{S}$. In
Eq.(1), $H_{c}$ depends only on the local operators within a cluster $c$,
while $H_{c,c^{\prime}},H_{c,c^{\prime},c^{\prime\prime}},\dots$ describe the
inter-cluster couplings. For systems with finite-range interactions, there are
proper partitions such that interactions only involve finitely many clusters.

Now, we formulate the DMMFT for the Hamiltonian in Eq.(1) Since no further
assumptions are made specifying the local operators or the microscopic details
of the coupling terms, DMMFT is generically applicable to fermions, bosons as
well as spins, irrespective of the presence of frustration. In conventional
MFTs, a local cluster is separated from the system and coupled to an effective
environment, where the environment is assumed to be classical, and the
correlated fluctuations between the cluster and the environment are neglected.
DMMFT improves the conventional mean-field approximation by including the
essential quantum fluctuations in the environment, where the reduced DM is
used to gauge the quantum fluctuations and select a Hilbert subspaces
approximating the effective environment.

In the following subsections, we will begin by reviewing the conventional
mean-field approximation in Sec. 2.1. Then, we will develop the mean-field
equations and self-consistency conditions of DMMFT in parallel to those in
conventional MFTs in Sec. 2.2. Following this, in Sec. 2.3, we will conduct a
comprehensive comparison between the DMMFT and the conventional MFTs. Within
this comparison, we will identify a hyperparameter, $n_{c}$, that interpolates
the DMMFT and the conventional MFTs. Particularly, when $n_{c}$ attains its
minimal value of $1$, DMMFT becomes equivalent to the conventional MFTs.

###  2.1 Conventional Mean-Field Approximation

We review the approximations employed in conventional MFTs before delving into
the development of DMMFT. In conventional MFTs, the mean-field decoupling
localizes the operator products to individual Hilbert subspaces by neglecting
the correlated fluctuations. More precisely, consider an operator product
$O_{i}^{\alpha}O_{j}^{\beta}$ acting on the Hilbert space
$\mathscr{H}_{i}\otimes\mathscr{H}_{j}$. The conventional mean-field
approximation decouples the operator product as follwos

$\begin{split}O_{i}^{\alpha}O_{j}^{\beta}\approx O_{i}^{\alpha}\langle
O_{j}^{\beta}\rangle+\langle O_{i}^{\alpha}\rangle O_{j}^{\beta}-\langle
O_{i}^{\alpha}\rangle\langle O_{j}^{\beta}\rangle,\end{split}$ (2)

where the terms on the right-hand side act on the subspace $\mathscr{H}_{i}$
or $\mathscr{H}_{j}$, or act trivially as an additive $c$-number.
Consequently, the correlated fluctuations $\langle\delta O_{i}^{\alpha}\delta
O_{j}^{\beta}\rangle=\langle O_{i}^{\alpha}O_{j}^{\beta}\rangle-\langle
O_{i}^{\alpha}\rangle\langle O_{j}^{\beta}\rangle$ vanish in this
approximation. Alternatively, from the perspective of the quantum states, as
the operator product factorizes in product states,
$\langle\phi_{i}|\otimes\langle\phi_{j}|O_{i}^{\alpha}O_{j}^{\beta}|\phi_{i}\rangle\otimes|\phi_{j}\rangle=\langle\phi_{i}|O_{i}^{\alpha}|\phi_{i}\rangle\langle\phi_{j}|O_{j}^{\beta}|\phi_{j}\rangle$,
the conventional MFTs implicitly assume the product structure of the states
and neglect quantum entanglements.

Keeping this consideration in mind, we formulate the conventional mean-field
approximation for the generic Hamiltonian in Eq.(1) as follows. Given a
cluster $c$, let $\mathcal{C}^{\prime}_{c}=\cup c^{\prime}$ be the collection
of the cluster connecting to $c$, referred to as the environment surrounding
$c$. The associated Hilbert spaces are $\mathscr{H}_{c}$ and
$\mathscr{H}_{\mathcal{C}^{\prime}_{c}}$ for the cluster and the environment,
respectively, where
$\mathscr{H}_{\mathcal{S}}=\otimes_{i\in\mathcal{S}}\mathscr{H}_{i}$ for a set
of sites $\mathcal{S}$. Retaining the terms within the extended cluster
$\bar{c}=c\cup\mathcal{C}^{\prime}_{c}$, the local Hamiltonian is

$\begin{split}H[\mathcal{O}_{\bar{c}}]=&H_{c}[\mathcal{O}_{c}]+H_{c,\mathcal{C}^{\prime}_{c}}[\mathcal{O}_{c},\mathcal{O}_{\mathcal{C}^{\prime}_{c}}]+H_{\mathcal{C}^{\prime}_{c}}[\mathcal{O}_{\mathcal{C}^{\prime}_{c}}],\\\
H_{c,\mathcal{C}^{\prime}_{c}}[\mathcal{O}_{c},\mathcal{O}_{\mathcal{C}^{\prime}_{c}}]=&\sum_{(c,c^{\prime}):c^{\prime}\in\mathcal{C}^{\prime}_{c}}H_{c,c^{\prime}}[\mathcal{O}_{c},\mathcal{O}_{c^{\prime}}]+\sum_{(c,c^{\prime},c^{\prime\prime}):c^{\prime},c^{\prime\prime}\in\mathcal{C}^{\prime}_{c}}H_{c,c^{\prime},c^{\prime\prime}}[\mathcal{O}_{c},\mathcal{O}_{c^{\prime}},\mathcal{O}_{c^{\prime\prime}}]+\dots,\\\
H_{\mathcal{C}^{\prime}_{c}}[\mathcal{O}_{\mathcal{C}^{\prime}_{c}}]=&\sum_{c^{\prime}\in\mathcal{C}^{\prime}_{c}}H_{c^{\prime}}[\mathcal{O}_{c^{\prime}}]+\sum_{(c^{\prime},c^{\prime\prime}):c^{\prime},c^{\prime\prime}\in\mathcal{C}^{\prime}_{c}}H_{c^{\prime},c^{\prime}}[\mathcal{O}_{c^{\prime}},\mathcal{O}_{c^{\prime\prime}}]+\dots,\end{split}$
(3)

where $H_{c}$ is the Hamiltonian of the focused cluster $c$,
$H_{\mathcal{C}^{\prime}_{c}}$ is the Hamiltonian of the environment
$\mathcal{C}^{\prime}_{c}$, and $H_{c,\mathcal{C}^{\prime}_{c}}$ represents
the couplings between the cluster $c$ and its environment
$\mathcal{C}^{\prime}_{c}$. If the target state can be approximated by a
product state locally, i.e.,

$\begin{split}|\phi_{\bar{c}}\rangle\approx|\phi_{\bar{c}}^{\text{MF}}\rangle=|\phi_{c}\rangle\otimes|\phi_{\mathcal{C}^{\prime}_{c}}\rangle,\end{split}$
(4)

a local effective Hamiltonian for the cluster $c$ can be obtained by
substituting the operators over $\mathscr{H}_{\mathcal{C}^{\prime}_{c}}$ with
their expectation values in $|\phi_{\mathcal{C}^{\prime}_{c}}\rangle$, i.e.,

$\begin{split}H_{\text{MF}}^{(c)}[\mathcal{O}_{c}]=&H_{c}[\mathcal{O}_{c}]+\langle\phi_{\mathcal{C}^{\prime}_{c}}|H_{c,\mathcal{C}^{\prime}_{c}}[\mathcal{O}_{c},\mathcal{O}_{\mathcal{C}^{\prime}_{c}}]|\phi_{\mathcal{C}^{\prime}_{c}}\rangle+\langle\phi_{\mathcal{C}^{\prime}_{c}}|H_{\mathcal{C}^{\prime}_{c}}[\mathcal{O}_{\mathcal{C}^{\prime}_{c}}]|\phi_{\mathcal{C}^{\prime}_{c}}\rangle,\end{split}$
(5)

which acts only on the Hilbert subspace $\mathscr{H}_{c}$. One solves
$|\phi_{c}\rangle$ as an eigenstate with respect to the effective Hamiltonian
$H_{\text{MF}}^{(c)}$ for each cluster $c$ locally. If one further assumes
that the target state of the entire system can also be approximated with a
product state, i.e.,

$\begin{split}|\Phi_{\text{MF}}\rangle=\otimes_{c}|\phi_{c}\rangle,\end{split}$
(6)

then, for each extended cluster,
$|\phi_{\mathcal{C}^{\prime}_{c}}\rangle=\otimes_{c^{\prime}\in\mathcal{C}^{\prime}_{c}}|\phi_{c^{\prime}}\rangle$.
Therefore, Eq.(4), Eq.(5), and Eq.(6) form a closed set of coupled mean-field
equations for the conventional MFT.

The mean-field equations can often be further simplified when systems have
additional symmetries. If the partition $\mathcal{C}=\\{c\\}$ respects some
symmetry of the system, there exist symmetry transformations relating the
clusters $T_{c^{\prime},c}:c\mapsto c^{\prime}$. A symmetric mean-field Ansatz
state can be constructed simply as
$|\Phi_{\text{MF}}\rangle=\otimes_{c}T_{c,c_{0}}|\phi_{c_{0}}\rangle$.
Particularly, for a translationally invariant state, all $|\phi_{c}\rangle$
can be chosen equal to the same state $|\phi_{c_{0}}\rangle$ in the local
Hilbert space $\mathscr{H}_{c_{0}}$, then the coupled mean-field equations for
the clusters reduce to a single set of self-consistent mean-field equations
for $|\phi_{c_{0}}\rangle$.

The major assumption in the conventional MFTs described above is that the
target states can be approximated with product states [Eqs.(4) and (6)], which
is, nevertheless, not justified a priori. For many ordered states, the
correlated fluctuations arising from the quantum entanglements are not
negligible, particularly at short range. Prototypically, in frustrated
magnetic systems and symmetry-protected topologically ordered systems, such
quantum fluctuations and quantum entanglements play essential roles. Improved
treatments of quantum fluctuations and quantum entanglements are needed to
better understand the emerging quantum ordered phases.

###  2.2 Density-Matrix Mean-Field Approximation

The simplifications achieved in conventional MFTs stem from the _separability_
, i.e., reducing the challenging task of studying the total Hamiltonian $H$
over an exponentially large Hilbert space $\mathscr{H}_{\mathcal{C}}$ to the
more manageable tasks of studying the local Hamiltonians $H^{(c)}_{\text{MF}}$
over local $\mathscr{H}_{c}$. However, it is not _necessary_ for a system to
be in a product state for two local subsystems to be separable. The
assumptions of the product states [Eqs.(4) and (6)] can thus be relaxed. More
precisely, for gapped systems, the correlated fluctuations $\langle\delta
O_{i}^{\alpha}\delta O_{j}^{\beta}\rangle\rightarrow 0$ as
$|i-j|\rightarrow\infty$. The absence of long-range entanglements ensures
separability, allowing the local physics to be approximated with effective
local systems. However, the short-range entanglements encode the quantum
fluctuations, demanding a more faithful treatment.

Instead of Eq.(2) for mean-field decoupling, it is instructive to recognize
that the entanglement of a local cluster with an (infinite or finitely large)
environment can always be _faithfully_ reproduced within a finite extension of
the cluster. More precisely, consider a generic state
$|\Psi\rangle\in\mathscr{H}_{\mathcal{C}}$. The entanglement of the state
$|\Psi\rangle$ over the cluster $c$ and the rest parts of the system
$\mathcal{C}\backslash c$ can be characterized by the reduced DM

$\begin{split}\rho_{c}=&\text{Tr}_{\mathcal{C}\backslash
c}|\Psi\rangle\langle\Psi|.\end{split}$ (7)

The entanglement is controlled by $\mathscr{H}_{c}$ than
$\mathscr{H}_{\mathcal{C}\backslash c}$ provided
$D_{\mathscr{H}_{c}}<D_{\mathscr{H}_{\mathcal{C}\backslash c}}$, where
$D_{\mathscr{H}}$ denotes the dimension of a Hilbert space $\mathscr{H}$.
According to the purification theorem [48], there exists a state
$|\tilde{\Psi}\rangle$ in
$\tilde{\mathscr{H}}=\mathscr{H}_{c}\otimes\tilde{\mathscr{H}}_{\tilde{c}}$
such that
$\text{Tr}_{\tilde{c}}|\tilde{\Psi}\rangle\langle\tilde{\Psi}|=\tilde{\rho}_{c}$
is equivalent to $\rho_{c}$, and the dimension of the Hilbert subspace
$\tilde{\mathscr{H}}_{\tilde{c}}$ is bounded by
$\tilde{D}^{E}_{c}=\lceil\exp({SE}_{c})\rceil$, where

$\begin{split}{SE}_{c}=-\text{Tr}_{c}\left[\rho_{c}\ln(\rho_{c})\right],\end{split}$
(8)

is the entanglement entropy and $\lceil q\rceil$ represents the smallest
integer larger than or equal to $q$. In this context, DMMFT seeks for an
effective Hamiltonian $\tilde{H}_{\text{MF}}^{(c)}$ over some
$\tilde{\mathscr{H}}$ (to be specified below), such that the reduced DM of the
state obtained from the the Hamiltonian,

$\begin{split}\rho_{c}^{\text{MF}}=\text{Tr}_{\tilde{c}}|\tilde{\phi}\rangle\langle\tilde{\phi}|,\end{split}$
(9)

well approximates the reduced DM of the target state
$\rho_{c}=\text{Tr}_{\mathcal{C}\backslash c}|\Phi\rangle\langle\Phi|$.

Instead of Eq.(4), DMMFT assumes that _for systems without long-range
entanglements, $\tilde{\mathscr{H}}$ can be found as a Hilbert subspace of
$\mathscr{H}_{\bar{c}}$ over the extended cluster $\bar{c}$._ Mathematically,
there exists a homomorphism

$\begin{split}\Pi=\text{Id}_{c}\otimes\Pi_{\mathcal{C}^{\prime}_{c}}:\tilde{\mathscr{H}}\rightarrow\mathscr{H}_{\bar{c}},\end{split}$
(10)

and its pseudo-inverse
$\Pi^{\dagger}:\mathscr{H}_{\bar{c}}\rightarrow\tilde{\mathscr{H}}$ is a
projector.

To avoid confusion, it is essential to emphasize that the
$\tilde{\mathscr{H}}$ constructed in DMMFT differs from that in DMET. The
purification theorem only ensures the existence of $\tilde{\mathscr{H}}$ with
its dimension bounded by $\tilde{D}_{c}^{E}$ from below. In DMET, an optimal
$\tilde{H}$ with the lowest possible dimension is employed, but at the cost of
a less straightforward construction of the embedding Hamiltonian. In contrast,
DMMFT can intuitively construct an effective Hamiltonian by restricting the
Hamiltonian for the extended cluster to $\tilde{\mathscr{H}}$, i.e.,

$\begin{split}\tilde{H}_{\text{MF}}^{(c)}[\mathcal{O}_{c},\tilde{\mathcal{O}}_{\mathcal{C}^{\prime}_{c}}]=\Pi^{\dagger}H[\mathcal{O}_{\bar{c}}]\Pi,\end{split}$
(11)

where the local operators in $\tilde{\mathcal{O}}_{\mathcal{C}^{\prime}_{c}}$
are

$\begin{split}\tilde{O}_{j}^{\beta}=&\Pi_{\mathcal{C}^{\prime}_{c}}^{\dagger}{O}_{j}^{\beta}\Pi_{\mathcal{C}^{\prime}_{c}},\forall
j\in\mathcal{C}^{\prime}_{c}.\end{split}$ (12)

The reduced density matrix of the target state is approximated with
$\rho_{c}^{\text{MF}}$ [Eq.(9)] for an eigenstate $|\tilde{\phi}\rangle$ of
$\tilde{H}_{\text{MF}}^{(c)}$.

Reciprocally, for each $\rho_{c}$, we may construct the local projectors
$\Pi_{c}$ as follows. Let
$\left\\{\lambda_{i}^{(c)},|\lambda_{i}^{(c)}\rangle\right\\}$ be the spectral
decomposition of $\rho_{c}$, and assume the state vectors are arranged in
decreasing order according to their eigenvalues, i.e.,
$\lambda_{i}^{(c)}\geq\lambda_{j}^{(c)},\forall i<j$. Define

$\begin{split}\Pi_{c}=\sum_{i=1}^{n_{c}}|\lambda_{i}^{(c)}\rangle\langle\lambda_{i}^{(c)}|,\end{split}$
(13)

where $n_{c}\leq D_{\mathscr{H}_{c}}$ is a cut-off parameter. For the extended
cluster $\bar{c}=c\cup\mathcal{C}^{\prime}_{c}$, we choose

$\begin{split}\Pi_{\mathcal{C}^{\prime}_{c}}=\otimes_{c^{\prime}\in\mathcal{C}^{\prime}_{c}}\Pi_{c^{\prime}}.\end{split}$
(14)

Parallel to the mean-field equations for the conventional MFTs, Eqs.(9) – (14)
constitute a closed set of coupled mean-field equations for the DMMFT. The
procedures for implementing the DMMFT algorithm are outlined in Appendix A. In
the presence of symmetries, the mean-field equations can be further
simplified. Specifically, for a translationally invariant target state
$|\Phi\rangle$, the reduced DMs of the local clusters are all identical to the
same reduced DM $\rho_{c_{0}}$ (as are the projectors). Consequently, the
coupled mean-field equations for the clusters reduce to a single set of self-
consistent mean-field equations for $\rho_{c_{0}}$.

###  2.3 Comparison of the Mean-Field Approximations

Comparing conventional MFTs with DMMFT, we find the cut-off parameter $n_{c}$
in Eq.(13) can be regarded as a hyperparameter interpolating the DMMFT and the
conventional MFT. To demonstrate this point, we first observe that when
$D_{\mathscr{H}_{\mathcal{C}^{\prime}_{c}}}=1$, DMMFT is equivalent to a
conventional MFT. If the ground state is a product state [Eq.(6)], the local
reduced density matrix of cluster $c$ is
$\rho_{c}=|\phi_{c}\rangle\langle\phi_{c}|$. The expectation values of the
local observables agree

$\begin{split}\langle
O_{i}^{\alpha}\rangle_{c}=\langle\phi_{c}|O_{i}^{\alpha}|\phi_{c}\rangle=\text{Tr}_{c}(\rho_{c}O_{i}^{\alpha}).\end{split}$
(15)

Moreover, since
$\Pi=\text{Id}_{c}\otimes\left(\otimes_{c^{\prime}}|\phi_{c}\rangle\langle\phi_{c}|\right)$,
$\tilde{\mathscr{H}}$ is isomorphic to $\mathscr{H}_{c}$, and the effective
Hamiltonians $\tilde{H}_{\text{MF}}^{(c)}$ are $H_{\text{MF}}^{(c)}$ identical
under this isomorphism.

Furthermore, from Eq.(14), we have
$D_{\mathscr{H}_{\mathcal{C}^{\prime}_{c}}}=\prod_{c^{\prime}\in\mathcal{C}^{\prime}_{c}}n_{c^{\prime}}$.
Therefore, we may take $n_{c}$ as a hyperparameter interpolating the
conventional MFT and the DMMFT. When $n_{c}=1$ (set to its minimal value), the
DMMFT simply reduces to the conventional MFT, which underestimates the quantum
fluctuations. When $n_{c}=D_{\mathscr{H}_{c}}$ (set to its maximal value), the
DMMFT describes a collection of (overlapped) extended cluster $\bar{c}$, which
overestimate the quantum fluctuations compared to the infinite system in the
thermodynamic limit. Therefore, as $n_{c}$ increases from $1$ to
$D_{\mathscr{H}_{c}}$, more and more quantum fluctuations are included.

Since $n_{c}$ gauges the amount of quantum fluctuations included in DMMFT from
underestimation to overestimation, it seems reasonable to expect an optimal
choice of its value, which leads to the following _conjecture_.

######  Conjecture.

The optimal choice of the hyperparameters is given by

$\begin{split}n_{c}^{*}=\tilde{D}^{E}_{c}=\lceil\exp({SE}_{c})\rceil.\end{split}$
(16)

We do not have a mathematical proof of this conjecture, which nevertheless
will not harm the practicality of the DMMFT. Particularly, one may take
$n_{c}=\tilde{D}^{E}_{c}$ as a rule-of-thumb choice and treat $n_{c}$ as a
_variational_ hyperparameter.

It is easy to observe that $\rho_{c}$ generically has non-vanishing
entanglement entropy whenever $D_{\mathscr{H}_{\mathcal{C}^{\prime}_{c}}}>1$.
This is exactly the short-range entanglement captured by the DMMFT beyond the
conventional MFTs. Not surprisingly, one expects that the DMMFT may detect
short-range entangled topological phases with a properly chosen $n_{c}$. In
Sec. 3.1, we demonstrate with the AKLT model that DMMFT can correctly detect
the topological ground states in terms of the entanglement spectra without
being biased to the symmetry-breaking states.

Lastly, we comment on the separability in DMMFT. In contrast to Eq.(6) for the
conventional MFTs, DMMFT does not assume the product structure of the state.
Instead, the weaker version of the separability adopted by DMMFT is

$\begin{split}\rho^{(2)}_{c_{1},c_{2}}=&\text{Tr}_{\mathcal{C}\backslash(c_{1}\cup
c_{2})}|\Phi\rangle\langle\Phi|\approx\rho_{c_{1}}\otimes\rho_{c_{2}},\end{split}$
(17)

for any two clusters $c_{1}$ amd $c_{2}$ that are not connected.

##  3 Applications

In this section, we apply the DMMFT to two demonstrative examples, the AKLT
model and the AFHTL. These two models have been extensively studied in the
literatures (cf. Ref. 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41,
42, 43, 44, 45, 46, 47 and the references therein). The AKLT model is an
exactly solvable model with topological ground states, despiet the vanishing
semiclassical (symmetry breaking) order parameters [29, 30]. On the other
hand, the AFHTL represents a prototypical model for frustrated magnetism. With
only nearest neighbor couplings and easy axis anisotropy, the AFHTL is known
to order in a three-sublattice ordered state homologous to the semiclassical
Néel state. However, the quantum fluctuations play non-trivial roles in
stabilizing the ordered state [45, 46, 47] and significantly renormalize the
order parameter [31]. In Sec. 3.1, we demonstrate that DMMFT can correctly
detects the topological ground states based on the entanglement spectra. In
Sec. 3.2, we compare the phases of the AFHTL in an external magnetic field
determined by DMMFT and conventional MFT.

###  3.1 Affleck-Kennedy-Lieb-Tasaki Model

We apply the DMMFT to the AKLT model. The AKLT model is a one-dimensional
spin-$1$ chain defined by the Hamiltonian

$\begin{split}H_{\text{AKLT}}=\sum_{i}P_{i,i+1}^{(S=2)}=\sum_{i}\frac{1}{2}\left[\mathbf{S}_{i}\cdot\mathbf{S}_{i+1}+\beta(\mathbf{S}_{i}\cdot\mathbf{S}_{i+1})^{2}+\gamma\right],\end{split}$
(18)

where $\beta$ describes the biquadratic coupling and $\gamma$ is an additive
constant, $\mathbf{S}_{i}$ are spin-$1$ operators whose representations in the
$S^{z}$-basis are

$\begin{split}S_{i}^{x}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&1&0\\\ 1&0&1\\\
0&1&0\end{pmatrix},S_{i}^{y}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&-\mathrm{i}&0\\\
\mathrm{i}&0&-\mathrm{i}\\\
0&\mathrm{i}&0\end{pmatrix},S_{i}^{z}=\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&-1\end{pmatrix}.\end{split}$ (19)

The coefficients $\beta$ and $\gamma$ are determined by the condition that
each term in the Hamiltonian is a projector to the Hilbert subspace of total
spin $S=2$. More precisely, when $\beta=1/3$, and $\gamma=2/3$,

$\begin{split}P^{(S=2)}_{i,i+1}(x)=\frac{1}{2}(x+\beta
x^{2}+\gamma)=&\begin{cases}1,&x=1,\\\ 0,&x=-1,-2.\\\ \end{cases}\end{split}$
(20)

where $x=\mathbf{S}_{i}\cdot\mathbf{S}_{i+1}=1,-1,-2$ in total spin sectors
$\mathscr{H}^{(S=2,1,0)}_{\\{i,i+1\\}}$ respectively. By construction,
$H_{\text{AKLT}}$ describes antiferromagnetic nearest neighbor couplings
penalizing only on neighboring spins being in $S=2$ states.

The AKLT model is an exactly solvable system, showcasing non-trivial
topological order. Notably, it features spin-$1/2$ (fractional to spin-$1$)
edge states and exhibits $4$-fold degeneracy of the ground states for an open
chain. This makes the AKLT model a valuable benchmark example.

Figure 1: Spectra of the two-site and single-site reduced DMs for the AKLT
model. Blue circles denote exact values, while red dots with error bars
represent DMMFT results.

Let us consider a two-site cluster $c$ situated in the bulk of an infinite
chain. By translational symmetry, any two-site cluster $c$ within the chain
shares the same reduced DM $\rho_{c}^{(2)}$. The spectrum of the two-site
reduced DM is exactly known and given by
$\text{eig}\rho_{c}^{(2)}=\left\\{\frac{1}{3},\frac{2}{9},\frac{2}{9},\frac{2}{9},0,0,0,0,0\right\\}$.
Furthermore, if we trace out one of the sites within the cluster, we obtain
the single-site reduced DM whose spectrum is exactly known,
$\text{eig}\rho_{c}^{(1)}=\left\\{\frac{1}{3},\frac{1}{3},\frac{1}{3}\right\\}$.

We apply the DMMFT to the AKLT model, focusing the two-site cluster $c$ and
solving the reduced DM $\rho_{c}^{(2)}$ self-consistently according to the
mean-field equations as described in Sec. 2.2. Specifically, the reduced DM is
calculated over an extended cluster encompassing the left and right nearest
neighboring two-site clusters to $c$. The associated Hilbert space is denoted
as
$\tilde{\mathscr{H}}_{\bar{c}}=\tilde{\mathscr{H}}_{l}\otimes\mathscr{H}_{c}\otimes\tilde{\mathscr{H}}_{r}$,
where $\tilde{\mathscr{H}}_{l,r}\sim\Pi_{c}(\mathscr{H}_{c})$ are the Hilbert
subspaces selected by the reduced DM, as per Eq.(13). We choose the cut-off
parameter $n_{c}=4$.

In Fig. 1, the spectra of the two-site and single-site reduced DMs are
presented. Blue circles represent exact values, while red dots depict DMMFT
results. The spectrum of the two-site reduced DM obtained with the DMMFT
reasonably aligns with the exact values [Fig. 1(a)], notably capturing the
correct degeneracies. The entanglement entropy of the two-site cluster
evaluated byDMMFT is $SE_{c,\text{MF}}^{(2)}=1.58$. Although it is slightly
larger than the exact value $SE_{c,\text{exact}}^{(2)}=1.3689$, this
discrepancy is expected due to the cut-off parameter $n_{c}$ chosen being
marginally larger than $\exp(SE_{c,\text{exact}}^{(2)})=3.9310$. Fig. 1(b)
displays the spectrum of the single-site reduced DM, where the DMMFT results
align excellently with the exact values. The single-site entanglement entropy
$SE_{c}^{(1)}=1.0986=\ln(3)$ also matches. Furthermore, the expectation values
of two spins within the cluster $\langle
S_{1,2}^{\alpha}\rangle=\text{Tr}_{c}(S_{1,2}^{\alpha}\rho_{c})$ both vanish.
This implies that the ground states found by DMMFT are not biased towards the
semiclassical Néel states, which agrees with the exact results.

###  3.2 Antiferromagnetic Heisenberg Model on Triangular Lattices

We apply the DMMFT to another illustrative example, the AFHTL. The Hamiltonian
for the AFHTL involving nearest neighbor anisotropic exchange interactions is
given by

$\begin{split}H_{AFH}=&\sum_{i,j}\left[J_{xy}\left(S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y}\right)+J_{z}S_{i}^{z}S_{j}^{z}\right]\\\
=&\sum_{i,j}J\left[\mathcal{A}\left(S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y}\right)+S_{i}^{z}S_{j}^{z}\right],\end{split}$
(21)

where $J_{z}$ and $J_{xy}$ are strengths of longitudinal and transverse
exchange interactions. For antiferromagnetic couplings, it is more convenient
to parameterize the model with $J=J_{z}>0$ and $\mathcal{A}=J_{xy}/J_{z}$. The
phases of the system only depend on the dimensionless parameter $\mathcal{A}$
that specifies the exchange anisotropy. Our focus is on the easy-axis case,
i.e., $0\leq\mathcal{A}\leq 1$. The two limiting points of $\mathcal{A}$ are:
1) $\mathcal{A}=0$, where the transverse exchange terms vanish (often referred
to as the Ising limit), and 2) $\mathcal{A}=1$, where the exchange
interactions have full rotational symmetry (often referred to as the
Heisenberg limit). Both of the limits limits have been extensively studied in
the literatures. Particularly, in the Ising limit, quantum fluctuations are
additionally introduced by the transverse magnetic field, and the sign problem
of the diagonal (longitudinal) frustrations can be mitigated in the QMC
simulations [47, 49]; while in the Heisenberg limit the system orders in the
so-called $120^{\circ}$-Néel state, as confirmed by various numerical methods
[31, 32, 33, 41, 42].

While the phase diagram is believed to be bounded by the Ising and Heisenberg
limits for a generic anisotropy parameter $0<\mathcal{A}<1$, computing the
details of quantum phases becomes challenging, particularly due to the QMC
sign problem when $\mathcal{A}\neq 0$. A study in Ref. 42 utilized large-size
cluster mean-field theory (belonging to the class of conventional MFTs) with a
scaling scheme, which provides a reasonable benchmark for the quantum phase
diagram. With DMMFT, we can more efficiently and quantitatively capture the
effects of quantum fluctuations with a small cluster of three sublattice
sites.

For illustrative purposes, we focus on the phases of the AFHTL with easy-axis
anisotropy subjected to a longitudinal external field, i.e.,

$\begin{split}H_{\text{AFH,Z}}=H_{\text{AFH}}+H_{Z},\end{split}$ (22)

where $H_{Z}=-\sum_{i}g\mu_{B}h_{z}S_{i}^{z}$ is the Zeeman coupling of the
spins to the magnetic field $h_{z}$. In our DMMFT calculations, we set the
gyromagnetic factor $g\mu_{B}=1$ and $h_{z}$ is measured in the corresponding
natural unit.

Four phases emerge in the AFHTL with generic easy-axis anisotropy under a
longitudinal external field. They are the coplanar ’“Y”-shaped phase, the
collinear up-up-down (UUD) phase, the coplanar “V”-shaped phase, and the
collinear polarized phase in the order of increasing $h_{z}$. All these phases
are compatible with a cluster partition of three-sublattice sites and
translational symmetries. A simple order parameter distinguishing these phases
is the magnetization along the easy axis, i.e.,

$\begin{split}M^{z}_{c}=\sum_{\alpha\in
c}S_{\alpha}^{z}=S_{A}^{z}+S_{B}^{z}+S_{C}^{z},\end{split}$ (23)

where the subscript $\alpha=A,B,C$ labels three sublattice sites within the
cluster $c$. The polarized and the UUD phases are characterized by
magnetization plateaux at $M^{z}_{c}=M_{z}^{c,(s)}$ and
$M^{z}_{c}=M_{z}^{c,(s)}/3$, respectively, where $M^{z}_{c,(s)}=\frac{3}{2}$
is the saturation magnetization. The “V”-shaped phase is sandwiched in between
the UUD and the polarized phases, while the “Y”-shaped phase appears at low
fields. A distinguishing order parameter characterizing the non-collinearity
of the “Y”-shaped phase is the vector chirality

$\begin{split}\bm{\kappa}_{V}=\frac{8}{3\sqrt{3}}\left(\mathbf{S}_{A}\times\mathbf{S}_{B}+\mathbf{S}_{B}\times\mathbf{S}_{C}+\mathbf{S}_{C}\times\mathbf{S}_{A}\right),\end{split}$
(24)

where the normalization factor is chosen such that $|\bm{\kappa}_{V}|=1$ for
the $120^{\circ}$-Néel state of $S=\frac{1}{2}$ spins without being
renormalized by quantum fluctuations.

Figure 2: Phase diagram of AFHTL in a longitudinal magnetic field
($\mathcal{A}=0.9$). (a) Vector chirality. (b) Magnetization. (c) Entanglement
Entropy. Red dots represent DMMFT results, while blue curves depict
conventional mean-field theory (CMFT) results. Vertical dashed lines indicate
phase boundaries. $h_{c1}$ distinguishes the “Y”-shaped phase from the UUD
phase, where DMMFT and CMFT are in agreement. $h_{c2}$ distinguishes the UUD-
shaped phase from the “V”-shaped phase. Determined from the magnetization
curves, the value of $h_{c2}$ in blue from CMFT exhibits a discernible
difference compared to $h^{\prime}_{c2}$ in red from DMMFT.

In our DMMFT calculations, we choose an extended cluster $\bar{c}$ comprising
of four three-sublattice clusters with periodic boundary conditions. Utilizing
periodic boundary conditions is to minimize the boundary effects, since the
Hilbert subspaces selected by the reduced DM is sensitive to the boundary
conditions as known in DMRG [50], While all four phases are homologous to
their corresponding semiclassical ordered states, the quantum fluctuations
play crucial roles in the “Y”-shaped and UUD phases. Specifically, quantum
fluctuations lift the accidental degeneracies in the semiclassical ground
state manifold, which is known as quantum order-by-disorder [45, 46], and
significantly renormalize the magnitudes of the ordered spins [31, 32, 33].

In Fig. 2 we present the dependence of the vector chirality, the
magnetization, and the entanglement entropy of the ground states of AFHTL
(with anisotropy parameter $\mathcal{A}=0.9$) on the longitudinal magnetic
field, calculated with the DMMFT (red) and the conventional MFT (blue). The
cut-off parameter $n_{c}$, as discussed in Sec. 2.3 interpolating the
conventional MFT and the DMMFT, is set to $n_{c}=4$ for the DMMFT and
$n_{c}=1$ for the conventional MFT.

In Fig. 2(a), the vector chirality calculated with conventional MFT decreases
monotonically as $h_{z}$ increases in the “Y”-shaped phase and vanishes in the
UUD and the “V”-shaped phases. In sharp contrast, the vector chirality
calculated with DMMFT, does not decrease monotonically in the “Y”-shaped
phase. Notably, the vector chirality at $h_{z}=0$ in DMMFT is only about
$30\%$ of that in conventional MFT. The discrepancy reflects the
renormalization of the ordered spins due to the quantum fluctuations.
Examining the magnitude of the ordered spins more closely , we find
$|\langle\mathbf{S}\rangle|=0.269\pm 0.014$ for small anisotropy
$\mathcal{A}=0.9$, This value is only about $50\%$ of $S=\frac{1}{2}$ but
agrees with the expectation from previous reported values [31, 32, 33] . Since
$\bm{\kappa}_{V}$ scales as $S^{2}$, a reduction factor of about $25\%$ is
expected for the magnitude of the vector chirality at $h_{z}=0$. The curves of
the DMMFT and the conventional MFT converge as $h_{z}$ increases and both
vanish beyond $h_{c1}$ (the critical field separating the “Y”-shaped and the
UUD phases determined according to the magnetization curves). With stronger
$h_{z}$, the spins are aligned closer to the easy axis, and the effect of the
quantum fluctuations weakens, which is clearly evidenced from the
monotonically decreasing entanglement entropy in the “Y”-shaped phase as shown
in Fig. 2(c). Therefore, the suppression of the quantum fluctuations and the
aligning of the spins towards the collinear UUD configuration by the
longitudinal magnetic field jointly lead to the non-monotonic dependence of
the vector chirality $|\bm{\kappa}_{V}|$ on $h_{z}$.

As $h_{z}$ surpasses $h_{c2}$, the system transitions towards the “V”-shaped
phase, deviating from the $\frac{1}{3}$ magnetization plateau of the UUD
phase. Fig. 2(b) demonstrates a discernible difference in the critical field
determined from the magnetization curves in DMMFT ($h^{\prime}_{c2}$) and in
convetional MFT ($h_{c2}$). Notably, the range of the field stabilizing the
UUD phase in DMMFT is larger than that in conventional MFT. Although the UUD
configuration is homologous to the semiclassical correspondence,quantum
fluctuations aid in stabilizing the UUD phase [42]. Especially in the
Heisenberg limit, the UUD phase is only stable at one point of
$h_{z}=h_{c2}=\frac{1}{3}h_{z}^{(s)}=1.5$ (where $h_{z}^{(s)}$ is the
saturation field) in the classical phase diagram, while it extends to a finite
range of the field in the quantum phase diagram. A similar effect of quantum
fluctuations is also expected when $\mathcal{A}<1$. Since $n_{c}$ used in
DMMFT is larger than that in conventional MFT, more quantum fluctuations are
included, a larger field range ($h^{\prime}_{c2}>h_{c2}$) that stabilizes the
UUD phase is expected.

Additionally, a self-consistency check for the choice of $n_{c}=4$ in DMMFT
can be confirmed from the entanglement entropy in the “Y”shaped phase. At
$h_{z}=0$, we have $SE_{c}\approx 1.8\ln(2)$. Therefore,
$\tilde{D}^{E}_{c}=\lceil\exp({SE}_{c})\rceil=4$.

##  4 Discussions

In Sec. 4.1, we undertake a comparative analysis of DMMFT with DMRG and DMFT.
Additionally, in Sec.4.2, we extend the DMMFT for systems at finite
temperatures, and in Sec. 4.3, we extend the DMMFT for systems with disorders.

###  4.1 Comparisons with Alternative Methods

We first compare DMMFT with DMRG. The constructions of the effective
Hamiltonian appear similar, where the reduced DM is employed to select Hilbert
subspace significant with respect to the target state. The resemblance between
the DMMFT and the infinite DMRG is most pronounced when the spatial dimension
$d=1$ as illustrated in the AKLT model in Sec. 3.1. The primary difference
lies in where the Hilbert space
$\tilde{\mathscr{H}}=\mathscr{H}_{c}\otimes\tilde{\mathscr{H}}_{\tilde{c}}$ is
chosen. In DMMFT, the projector $\Pi_{\mathcal{C}^{\prime}_{c}}$ is
iteratively optimized solely over the Hilbert space
$\mathscr{H}_{\mathcal{C}^{\prime}_{c}}$ without changing the cluster size of
the environment. Conversely, in infinite DMRG, the environment is constructed
iteratively over the Hilbert space of the semi-infinite chains. Gauged by the
quantum entanglement, the dimension of the Hilbert subspace selected for
constructing the effective Hamiltonian is bounded by $\tilde{D}^{E}_{c}$
[Eq.(16)] in both cases. As $\tilde{D}^{E}_{c}$ is bounded by the size of the
cluster $c$, DMRG can feasibly select a finite significant Hilbert subspace
without being hindered by the infinite size of the environment asymptotically.
On the contrary, in DMMFT, we assume that the significant Hilbert subspace
resides within $\mathscr{H}_{\mathcal{C}^{\prime}_{c}}$, which is a reasonable
approximation, especially for short-range entangled systems. Ignoring
practical constraints provided we choose the cluster size being larger than
the typical entanglement range. the DMMFT and the DMRG are equivalent in
$d=1$.

The advantages of DMMFT become more significant in higher dimensions, $d\geq
2$. In DMRG, a $d$-dimensional system is compactified to a quasi-one-
dimensional system, and an artificial one-dimensional chain of clusters winds
through the system. Consequently, even in the case of a short-range entangled
system, artificial long-range entanglements emerge for sites physically close
but distant along the chain. In contrast, DMMFT makes use of the separability
[Eq.(17)] enables the study of local problems over extended clusters without
introducing artificial long-range entanglements arising from the transverse
dimensions. From this perspective, DMMFT may be regarded as a mean-field
approximation of DMRG and proves particularly valuable for short-range
entangled systems in higher dimensions.

Now, we turn to compare DMMFT with DMFT. In DMFT, designed for fermionic
systems, the effective action is constructed over a local cluster, employing
the single-particle Green’s function as a dynamical mean-field to capture
fluctuations [19, 20]. This local problem, entailing a dynamical mean-field,
is subsequently mapped to an effective impurity model and solved using ED or
QMC kernels [19]. In contrast to the Anderson impurity model tailored for
fermions, establishing a comparable impurity model for spins is not a
straightforward task. Furthermore, even if a spin impurity model aligned with
the principles of DMFT could be formulated, the challenge persists in
efficiently solving it. Constructing a generic spin impurity (if not
impossible) can introduce significant frustrations, which leads to the failure
of the QMC kernel due to the sign problem. On the other hand, achieving an
accurate representation of the environment requires incorporating a large
number of spins beyond those in the local cluster, which causes the ED kernel
to rapidly encounter the exponential wall. A simplification of DMFT proposed
in DMET for fermionic systems involves utilizing the local density matrix, a
static observable, instead of the dynamical Green’s function [21]. While a
formal generalization of DMET to spin systems is feasible, the challenge lies
in identifying a suitable reference wave-function that facilitates
straightforward Schmidt decomposition and the construction of the embedding
Hamiltonian. On the contrary, in DMMFT, the reduced density matrix is
generally defined for many-body states, making DMMFT naturally applicable to
fermions, bosons, as well as spins.

An alternative approach to investigating the ground state of a quantum many-
body system is to work directly with the wave function. One example can be the
Gutzwiller method, which uses Gutzwiller wave-function as an variational
Ansatz [51, 52, 53, 54]. The Gutzwiller approximation can be viewed as a
quantum embedding approach in a unified perspective with DMFT and DMET [55].
More recently, VDAT has been proposed to solve the shortcomings of Gutzwiller
Approximation and simplify the computational complexity of DMFT. In VDAT, the
variational ansatz is the sequential product density matrix, which is
evaluated via the discrete action theory [22, 23]. The VDAT Ansatz is
controlled by an integral parameter $\mathcal{N}$, and it is asymptotically
exact as $\mathcal{N}$ tends to infinity but also becomes more computationally
demanding. While it has been demonstrated that VDAT can yield accurate ground
state wave functions for fermionic systems, even with $\mathcal{N}$ as small
as $3$ [24, 25], it faces challenges when applied to spin systems.
Constructing the sequential product density matrix for the spin system is
straightforward; however, the absence of Wick’s theorem for spin systems
presents a challenge in efficiently applying discrete action theory to
evaluate the Ansatz.

Integrating wave function approaches with DMMFT raises the question of
determining a global state from reduced density matrices of local clusters.
This challenge is recognized as the quantum marginal problem [56]. Finding a
general solution to this problem is currently beyond the scope of this paper
and remains an open and challenging research question [57].

###  4.2 Systems at Finite Temperatures

The extension of DMMFT for systems at finite temperatures is straightforward.
In Eq.(9), we use the ground state of the effective Hamiltonian over the
extended cluster to compute the reduced DM. At finite temperatures, assuming
local thermal equilibrium, we can average the reduced DMs of the eigenstates
weighted by the Boltzmann distribution. Specifically, let
$\\{E_{\alpha},|\tilde{\phi}_{\alpha}\rangle\\}$ be the eigensystems of
$\tilde{H}_{\text{MF}}^{(c)}$. Then, the thermally averaged reduced DM is
given by

$\begin{split}[\rho_{c}^{\text{MF}}]_{T}=&\sum_{\alpha}w_{\alpha}\rho_{c}^{\alpha},\\\
w_{\alpha}=&\frac{\exp[-E_{\alpha}/(k_{B}T)]}{\sum_{\beta}\exp[-E_{\beta}/(k_{B}T)]},\\\
\rho_{c}^{\alpha}=&\text{Tr}_{\tilde{c}}|\tilde{\phi}_{\alpha}\rangle\langle\tilde{\phi}_{\alpha}|,\end{split}$
(25)

where $k_{B}$ is the Boltzmann constant, and the square bracket with a
subscript $T$ indicates the thermal average. As discussed in Sec. 2.3, when
$D_{\mathscr{H}_{\mathcal{C}^{\prime}_{c}}}=1$, DMMFT reduces to the CMFTs at
zero temperature. It is straightforward to check that, since
$[\rho_{c}^{\text{MF}}]_{T}\sim\exp[-\tilde{H}_{\text{MF}}^{(c)}/(k_{B}T)]$
when $D_{\mathscr{H}_{\mathcal{C}^{\prime}_{c}}}=1$, the extension of DMMFT
according Eq.(25) also aligns with the conventional MFTs at finite
temperatures.

However, it is worth noting that the thermal averaging in Eq.(25) does not
fully capture all thermal fluctuations. Specifically, at low temperatures, it
is often the gapless long wave length modes (Goldstone modes) that dominate
the thermal fluctuations. However, due to the mean-field approximation
inherent in DMMFT, these thermal fluctuations are neglected. Therefore,
similar to all other mean-field methods, DMMFT tends to underestimate the
impact of thermal fluctuations at finite temperatures.

###  4.3 Systems with Disorders

Disorders are an inevitable aspect of experiments. Even minuscule amounts of
disorders can alter the nature of fragile quantum phases, particularly in the
frustrated systems. Consequently, it becomes crucial to incorporate disorders
into the DMMFT framework.

In the presence of the disorders, the Hamiltonian of the extended cluster
undergoes direct modification, denoted as

$\begin{split}H[\mathcal{O}_{\bar{c}}]&\rightarrow
H^{(\delta)}[\mathcal{O}_{\bar{c}}],\end{split}$ (26)

where the superscript $\delta$ labels the types of the disorders. However, the
process of averaging over disorders is more intricate than the thermal average
discussed in Sec. 4.2.

Consider a straightforward scenario involving a single magnetic vacancy in the
cluster $c$, wherein there can be either no vacancies or one vacancy on a site
$i\in c$. In this case, we denote $\delta\in\Delta=\\{0\\}\cup c$, where
$\delta=0$ signifies no vacancies in the cluster $c$. We assume the
probability measure of the disorder configurations to be
$P_{\Delta}(\delta)\mathrm{d}\delta$. Analogous to to Eq.(25), it is intuitive
to to define an average over the disorders as

$\begin{split}[\rho_{c}^{\text{MF}}]_{\Delta,\text{ann}}=&\int_{\Delta}\rho_{c}^{(\delta)}P_{\Delta}(\delta)\mathrm{d}\delta,\\\
\rho_{c}^{(\delta)}=&\text{Tr}_{\tilde{c}}|\tilde{\phi}_{(\delta)}\rangle\langle\tilde{\phi}_{(\delta)}|,\end{split}$
(27)

where $\tilde{\phi}_{(\delta)}$ is the ground state of the effective
Hamiltonian $\tilde{H}_{\text{MF}}^{(c),(\delta)}$, and the square bracket
with a subscript $\Delta$ indicates the average over the disorder
configurations. The physical meaning of
$[\rho_{c}^{\text{MF}}]_{\Delta,\text{ann}}$ so-defined requires further
clarification. When $[\rho_{c}^{\text{MF}}]_{\Delta,\text{ann}}$ is employed
to set the effective environment $\mathcal{C}^{\prime}_{c}=\\{c^{\prime}\\}$
for the cluster $c$, each $c^{\prime}$ is already averaged over various
disorder configurations according to Eq.(27). In other words, the disorders
are _annealed_. For this reason, we also add the subscript “ann” in Eq.(27).

In experiments, magnetic vacancies often exhibit behavior more akin to static
disorders than fluctuating disorders, particularly when the vacancies are not
thermally migrating. Disorders as such are often referred to as _quenched_
disorders, in contrast to the annealed disorders. In this case, we consider
all possible disorder configurations over the extended cluster $\bar{c}$.
Assume $\bar{c}$ consisting of $\nu$ simple clusters, including the focused
cluster $c$ and $\nu-1$ clusters in the environment. If the simple clusters
are identical, all possible disorder configurations of the quenched disorders
form a set
$\bm{\Delta}=\Delta^{\nu}=(\delta_{1},\delta_{2},\dots,\delta_{\nu})$. We
obtain coupled mean-field equations for
$\\{\rho_{c}^{(\bm{\delta})}\\}_{\bm{\Delta}}$ . Once
$\\{\rho_{c}^{(\bm{\delta})}\\}_{\bm{\Delta}}$ is solved, the expectation
values of local observables averaged over the quenched disorders should be
calculated as

$\begin{split}[\langle
O_{i}^{\alpha}\rangle]_{\bm{\Delta},\text{quen}}=\int_{\bm{\Delta}}\text{Tr}_{c}\left[O_{i}^{\alpha}\rho_{c}^{(\bm{\delta})}\right]P_{\bm{\Delta}}(\bm{\delta})(\mathrm{d}^{\nu}\bm{\delta}),\end{split}$
(28)

where $P_{\bm{\Delta}}(\bm{\delta})(\mathrm{d}^{\nu}\bm{\delta})$ is the
probability measure over the configuration space of the quenched disorders
$\bm{\Delta}$. In the special case where each vacancy is independent and
subjected to an identical distribution $P_{\Delta}(\delta)\mathrm{d}\delta$,
then
$P_{\bm{\Delta}}(\bm{\delta})(\mathrm{d}^{\nu}\bm{\delta})=\prod_{n=1}^{\nu}P_{\Delta}(\delta_{n})\mathrm{d}\delta_{n}$.

In real systems, the effects of disorders can be more complicated. For
instance, there can be vacancies clusters due to lower formation free energy,
and the probability measure of the vacancies no longer takes a simple product
form. Moreover, disorders can manifest as inhomogeneity of the samples,
especially when scale of the probe is small, and the effects of disorders are
not self-averaged in the experiments. Despite of the complications of
disorders arising in real systems, DMMFT can, in principle, be adapted to
include the disorder effects accordingly.

##  5 Conclusion

In this study, we introduce a novel mean-field approach, DMMFT, specifically
designed to capture quantum fluctuations beyond conventional MFTs. Combining
the strengths of DMRG and DMFT, DMMFT stands out as an efficient and unbiased
method applicable to fermions, bosons, and spins, even in the presence of
frustrations. A noteworthy feature of DMMFT is its ability to discern
topological phases, a capability further enhanced by its utilization of the
reduced density matrix. Additionally, DMMFT seamlessly extends to systems at
finite temperatures and with disorders. In conclusion, DMMFT provides an
effective tool for exploring phases exhibiting unconventional quantum orders,
particularly beneficial for investigating frustrated spin systems in high
spatial dimensions.

## Acknowledgements

This work is supported as part of the Institute for Quantum Matter, an Energy
Frontier Research Center funded by the U.S. Department of Energy, Office of
Science, Office of Basic Energy Sciences, under Award DESC0019331. This work
is also partially supported by the NSF CAREER Grant DMR-1848349. JYZ is
additionally supported by TIPAC from the Department of Physics and Astronomy,
Johns Hopkins University. ZC is supported by the Columbia Center for
Computational Electrochemistry. JYZ acknowledges motivating discussions with
Collin Broholm, Oleg Tchernyshyov, and Tong Chen. JYZ also thanks Salvia
Felixsen Skogen for the hospitality at Domini Katzen Institute.

##  Appendix A Iterative Algorithm for Mean-Field Equations

In Section 2.2, we established a closed set of coupled mean-field equations
for DMMFT, delineated by Eqs. (9) to (14). Here, we present an iterative
algorithm, as outlined in Table 1, to achieve self-consistent solutions
numerically.

Step number | Step description
---|---
0 | Initialize a set of trial reduced DMs $\\{\rho_{c}\\}_{c\in\mathcal{C}}$.
For each cluster $c$: |
1 | Diagonalize $\rho_{c}$, and construct $\Pi_{c}$ according to Eq.(13).
2 | Use Eq.(14) to construct $\Pi_{\mathcal{C}^{\prime}_{c}}$.
3 | Use Eq.(10) to construct $\Pi$.
4 | Use Eq.(11) to construct $\tilde{H}_{\text{MF}}^{(c)}$.
5 | Diagonalize $\tilde{H}_{\text{MF}}^{(c)}$ and select the target state $|\tilde{\phi}_{c}\rangle$.
6 | Calculate new $\rho_{c}$ according to Eq.(9).
0’ | Use new $\\{\rho_{c}\\}_{c\in\mathcal{C}}$ to update the trial reduced DMs,
| and repeat steps 1 - 6 until convergence.
Table 1: DMMFT implementation steps.

There remain some noteworthy considerations. In Step 0, while a stable self-
consistent solution should be independent of initial values, providing an
initial guess close to the self-consistent solution often aids convergence. An
empirical approach is to start with a trivial ${\rho_{c}}$ that can be
constructed from solutions of conventional MFTs. In Step 0’, a linear
interpolation of the reduced DM

$\begin{split}\rho_{c}^{(n+1)}=(1-\alpha)\rho_{c}^{(n)}+\alpha\rho_{c}^{(n,\text{new})},\end{split}$
(29)

can be employed for updating $\rho_{c}$. Here, $\rho_{c}^{(n)}$ represents the
reduced DM used in Step 0 at the $n$-th iteration, and
$\rho_{c}^{(n,\text{new})}$ is the new reduced DM obtained in Step 0’, The
learning rate $\alpha$ varies in the range $(0,1)$. Larger values of $\alpha$
lead to a more rapid update, while smaller values contribute to stabilizing
the iterative process.

Finally, the hyperparameter $n_{c}$ used in Eq.(13) in Step 1 can be compared
to $\tilde{D}^{E}_{c}=\lceil\exp({SE}_{c})\rceil$, where the entanglement
entropy is calculated according to Eq.(8). It is important to note that
setting $n_{c}=1$ should recover the results of conventional MFTs.

#### Funding information

JYZ is supported by IQM EFRC DOE under Award DESC0019331, NSF CAREER under
Grant DMR-1848349, and TIPAC PHA JHU. ZC is supported by the Columbia Center
for Computational Electrochemistry.

## References

  * [1] L. Balents,  _Spin liquids in frustrated magnets_ ,  Nature 464(7286), 199 (2010),  10.1038/nature08917.
  * [2] M. J. P. Gingras and P. A. McClarty,  _Quantum spin ice: a search for gapless quantum spin liquids in pyrochlore magnets_ ,  Reports on Progress in Physics 77(5), 056501 (2014),  10.1088/0034-4885/77/5/056501.
  * [3] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman and T. Senthil,  _Quantum spin liquids_ ,  Science 367(6475), eaay0668 (2020),  10.1126/science.aay0668,  https://www.science.org/doi/pdf/10.1126/science.aay0668.
  * [4] P. W. Anderson,  _Resonating valence bonds: A new kind of insulator?_ ,  Mater. Res. Bull. 8, 153 (1973),  10.1016/0025-5408(73)90167-0.
  * [5] P. W. A. P. Fazekas,  _On the ground state properties of the anisotropic triangular antiferromagnet_ ,  Philos. Mag. 30, 423 (1974),  10.1080/14786439808206568.
  * [6] V. Kalmeyer and R. B. Laughlin,  _Equivalence of the resonating-valence-bond and fractional quantum hall states_ ,  Phys. Rev. Lett. 59, 2095 (1987),  10.1103/PhysRevLett.59.2095.
  * [7] A. Kitaev,  _Anyons in an exactly solved model and beyond_ ,  Annals of Physics 321(1), 2 (2006),  https://doi.org/10.1016/j.aop.2005.10.005,  January Special Issue.
  * [8] C. Castelnovo, R. Moessner and S. L. Sondhi,  _Magnetic monopoles in spin ice_ ,  Nature 451, 42 (2008),  10.1038/nature06433.
  * [9] P. W. Anderson,  _The resonating valence bond state in $la_{2}cuo_{4}$ and superconductivity_,  Science 235(4793), 1196 (1987),  10.1126/science.235.4793.1196,  https://www.science.org/doi/pdf/10.1126/science.235.4793.1196.
  * [10] A. Kitaev,  _Fault-tolerant quantum computation by anyons_ ,  Annals of Physics 303(1), 2 (2003),  https://doi.org/10.1016/S0003-4916(02)00018-0.
  * [11] S. R. White,  _Density matrix formulation for quantum renormalization groups_ ,  Phys. Rev. Lett. 69, 2863 (1992),  10.1103/PhysRevLett.69.2863.
  * [12] H. C. Jiang, Z. Y. Weng and D. N. Sheng,  _Density matrix renormalization group numerical study of the kagome antiferromagnet_ ,  Phys. Rev. Lett. 101, 117203 (2008),  10.1103/PhysRevLett.101.117203.
  * [13] S. Yan, D. A. Huse and S. R. White,  _Spin-liquid ground state of the $s=1/2$ kagome Heisenberg antiferromagnet_,  Science 332(6034), 1173 (2011),  10.1126/science.1201080,  https://www.science.org/doi/pdf/10.1126/science.1201080.
  * [14] Z. Zhu, P. A. Maksimov, S. R. White and A. L. Chernyshev,  _Topography of spin liquids on a triangular lattice_ ,  Phys. Rev. Lett. 120, 207203 (2018),  10.1103/PhysRevLett.120.207203.
  * [15] G. Senatore and N. H. March,  _Recent progress in the field of electron correlation_ ,  Rev. Mod. Phys. 66, 445 (1994),  10.1103/RevModPhys.66.445.
  * [16] A. W. Sandvik,  _Computational Studies of Quantum Spin Systems_ ,  AIP Conference Proceedings 1297(1), 135 (2010),  10.1063/1.3518900,  https://pubs.aip.org/aip/acp/article-pdf/1297/1/135/11407753/135_1_online.pdf.
  * [17] P. W. Anderson,  _Localized magnetic states in metals_ ,  Phys. Rev. 124, 41 (1961),  10.1103/PhysRev.124.41.
  * [18] J. Bardeen, L. N. Cooper and J. R. Schrieffer,  _Theory of superconductivity_ ,  Phys. Rev. 108, 1175 (1957),  10.1103/PhysRev.108.1175.
  * [19] A. Georges, G. Kotliar, W. Krauth and M. J. Rozenberg,  _Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions_ ,  Rev. Mod. Phys. 68, 13 (1996),  10.1103/RevModPhys.68.13.
  * [20] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet and C. A. Marianetti,  _Electronic structure calculations with dynamical mean-field theory_ ,  Rev. Mod. Phys. 78, 865 (2006),  10.1103/RevModPhys.78.865.
  * [21] G. Knizia and G. K.-L. Chan,  _Density matrix embedding: A simple alternative to dynamical mean-field theory_ ,  Phys. Rev. Lett. 109, 186404 (2012),  10.1103/PhysRevLett.109.186404.
  * [22] Z. Cheng and C. A. Marianetti,  _Variational discrete action theory_ ,  Phys. Rev. Lett. 126, 206402 (2021),  10.1103/PhysRevLett.126.206402.
  * [23] Z. Cheng and C. A. Marianetti,  _Foundations of variational discrete action theory_ ,  Phys. Rev. B 103, 195138 (2021),  10.1103/PhysRevB.103.195138.
  * [24] Z. Cheng and C. A. Marianetti,  _Precise ground state of multiorbital mott systems via the variational discrete action theory_ ,  Phys. Rev. B 106, 205129 (2022),  10.1103/PhysRevB.106.205129.
  * [25] Z. Cheng and C. A. Marianetti,  _Gauge constrained algorithm of variational discrete action theory at $\mathcal{N}=3$ for the multiorbital hubbard model_,  Phys. Rev. B 108, 035127 (2023),  10.1103/PhysRevB.108.035127.
  * [26] T. Holstein and H. Primakoff,  _Field dependence of the intrinsic domain magnetization of a ferromagnet_ ,  Phys. Rev. 58, 1098 (1940),  10.1103/PhysRev.58.1098.
  * [27] D. P. Arovas and A. Auerbach,  _Functional integral theories of low-dimensional quantum Heisenberg models_ ,  Phys. Rev. B 38, 316 (1988),  10.1103/PhysRevB.38.316.
  * [28] F. D. M. Haldane,  _Nonlinear field theory of large-spin Heisenberg antiferromagnets: Semiclassically quantized solitons of the one-dimensional easy-axis néel state_ ,  Phys. Rev. Lett. 50, 1153 (1983),  10.1103/PhysRevLett.50.1153.
  * [29] I. Affleck, T. Kennedy, E. H. Lieb and H. Tasaki,  _Rigorous results on valence-bond ground states in antiferromagnets_ ,  Phys. Rev. Lett. 59, 799 (1987),  10.1103/PhysRevLett.59.799.
  * [30] I. Affleck, T. Kennedy, E. H. Lieb and H. Tasaki,  _Valence bond ground states in isotropic quantum antiferromagnets_ ,  Communications in Mathematical Physics 115(3), 477 (1988).
  * [31] D. A. Huse and V. Elser,  _Simple variational wave functions for two-dimensional Heisenberg spin-½ antiferromagnets_ ,  Phys. Rev. Lett. 60, 2531 (1988),  10.1103/PhysRevLett.60.2531.
  * [32] T. Jolicoeur, E. Dagotto, E. Gagliano and S. Bacci,  _Ground-state properties of the s=1/2 Heisenberg antiferromagnet on a triangular lattice_ ,  Phys. Rev. B 42, 4800 (1990),  10.1103/PhysRevB.42.4800.
  * [33] S. R. White and A. L. Chernyshev,  _Néel order in square and triangular lattice Heisenberg models_ ,  Phys. Rev. Lett. 99, 127004 (2007),  10.1103/PhysRevLett.99.127004.
  * [34] Z. Zhu and S. R. White,  _Spin liquid phase of the $s=\frac{1}{2}\phantom{\rule{4.0pt}{0.0pt}}{J}_{1}-{J}_{2}$ Heisenberg model on the triangular lattice_,  Phys. Rev. B 92, 041105 (2015),  10.1103/PhysRevB.92.041105.
  * [35] W.-J. Hu, S.-S. Gong, W. Zhu and D. N. Sheng,  _Competing spin-liquid states in the spin- $\frac{1}{2}$ Heisenberg model on the triangular lattice_,  Phys. Rev. B 92, 140403 (2015),  10.1103/PhysRevB.92.140403.
  * [36] Y. Iqbal, W.-J. Hu, R. Thomale, D. Poilblanc and F. Becca,  _Spin liquid nature in the Heisenberg ${J}_{1}-{J}_{2}$ triangular antiferromagnet_,  Phys. Rev. B 93, 144411 (2016),  10.1103/PhysRevB.93.144411.
  * [37] A. Wietek and A. M. Läuchli,  _Chiral spin liquid and quantum criticality in extended $s=\frac{1}{2}$ Heisenberg models on the triangular lattice_,  Phys. Rev. B 95, 035141 (2017),  10.1103/PhysRevB.95.035141.
  * [38] S. Hu, W. Zhu, S. Eggert and Y.-C. He,  _Dirac spin liquid on the spin- $1/2$ triangular Heisenberg antiferromagnet_,  Phys. Rev. Lett. 123, 207203 (2019),  10.1103/PhysRevLett.123.207203.
  * [39] S.-S. Gong, W. Zheng, M. Lee, Y.-M. Lu and D. N. Sheng,  _Chiral spin liquid with spinon fermi surfaces in the spin- $\frac{1}{2}$ triangular Heisenberg model_,  Phys. Rev. B 100, 241111 (2019),  10.1103/PhysRevB.100.241111.
  * [40] Y.-F. Jiang and H.-C. Jiang,  _Nature of quantum spin liquids of the $s=\frac{1}{2}$ Heisenberg antiferromagnet on the triangular lattice: A parallel dmrg study_,  Phys. Rev. B 107, L140411 (2023),  10.1103/PhysRevB.107.L140411.
  * [41] L. Seabra, T. Momoi, P. Sindzingre and N. Shannon,  _Phase diagram of the classical Heisenberg antiferromagnet on a triangular lattice in an applied magnetic field_ ,  Phys. Rev. B 84, 214418 (2011),  10.1103/PhysRevB.84.214418.
  * [42] D. Yamamoto, G. Marmorini and I. Danshita,  _Quantum phase diagram of the triangular-lattice $xxz$ model in a magnetic field_,  Phys. Rev. Lett. 112, 127203 (2014),  10.1103/PhysRevLett.112.127203.
  * [43] A. V. Chubukov and D. I. Golosov,  _Quantum theory of an antiferromagnet on a triangular lattice in a magnetic field_ ,  Journal of Physics: Condensed Matter 3(1), 69 (1991),  10.1088/0953-8984/3/1/005.
  * [44] A. V. Chubukov and T. Jolicoeur,  _Order-from-disorder phenomena in Heisenberg antiferromagnets on a triangular lattice_ ,  Phys. Rev. B 46, 11137 (1992),  10.1103/PhysRevB.46.11137.
  * [45] B. Kleine, E. Müller-Hartmann, K. Frahm and P. Fazekas,  _Spin-wave analysis of easy-axis quantum antiferromagnets on the triangular lattice_ ,  Zeitschrift für Physik B Condensed Matter 87, 103 (1992),  10.1007/BF01308264.
  * [46] Q. Sheng and C. L. Henley,  _Ordering due to disorder in a triangular Heisenberg antiferromagnet with exchange anisotropy_ ,  Journal of Physics: Condensed Matter 4(11), 2937 (1992),  10.1088/0953-8984/4/11/020.
  * [47] R. Moessner and S. L. Sondhi,  _Ising models of quantum frustration_ ,  Phys. Rev. B 63, 224401 (2001),  10.1103/PhysRevB.63.224401.
  * [48] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information_ ,  Cambridge University Press,  10.1017/CBO9780511976667 (2010).
  * [49] Y. Da Liao, H. Li, Z. Yan, H.-T. Wei, W. Li, Y. Qi and Z. Y. Meng,  _Phase diagram of the quantum ising model on a triangular lattice under external field_ ,  Phys. Rev. B 103, 104416 (2021),  10.1103/PhysRevB.103.104416.
  * [50] S. R. White and R. M. Noack,  _Real-space quantum renormalization groups_ ,  Phys. Rev. Lett. 68, 3487 (1992),  10.1103/PhysRevLett.68.3487.
  * [51] M. C. Gutzwiller,  _Correlation of electrons in a narrow $s$ band_,  Phys. Rev. 137, A1726 (1965),  10.1103/PhysRev.137.A1726.
  * [52] W. Metzner and D. Vollhardt,  _Analytic calculation of ground-state properties of correlated fermions with the gutzwiller wave function_ ,  Phys. Rev. B 37, 7382 (1988),  10.1103/PhysRevB.37.7382.
  * [53] X. Deng, L. Wang, X. Dai and Z. Fang,  _Local density approximation combined with gutzwiller method for correlated electron systems: Formalism and applications_ ,  Phys. Rev. B 79, 075114 (2009),  10.1103/PhysRevB.79.075114.
  * [54] N. Lanatà, H. U. R. Strand, X. Dai and B. Hellsing,  _Efficient implementation of the gutzwiller variational method_ ,  Phys. Rev. B 85, 035133 (2012),  10.1103/PhysRevB.85.035133.
  * [55] T. Ayral, T.-H. Lee and G. Kotliar,  _Dynamical mean-field theory, density-matrix embedding theory, and rotationally invariant slave bosons: A unified perspective_ ,  Phys. Rev. B 96, 235139 (2017),  10.1103/PhysRevB.96.235139.
  * [56] C. Schilling,  _The quantum marginal problem_ ,  https://arxiv.org/abs/1404.1085.
  * [57] X.-D. Yu, T. Simnacher, N. Wyderka, H. C. Nguyen and O. Gühne,  _A complete hierarchy for the pure state marginal problem in quantum mechanics_ ,  Nature Communications 12(1), 1012 (2021),  10.1038/s41467-020-20799-5.