# Frustrated spin-$\frac{1}{2}$ Heisenberg magnet on an $AA$-stacked honeycomb
bilayer: High-order study of the collinear magnetic phases of the
$J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model

P H Y Li1,2 and R F Bishop1,2<EMAIL_ADDRESS><EMAIL_ADDRESS>1 School of Physics and Astronomy, Schuster
Building, The University of Manchester, Manchester, M13 9PL, UK  2 School of
Physics and Astronomy, University of Minnesota, 116 Church Street SE,
Minneapolis, Minnesota 55455, USA

###### Abstract

The zero-temperature phase diagram of the frustrated spin-$\frac{1}{2}$
$J_{1}$–$J_{2}$–$J_{1}^{\perp}$ Heisenberg magnet on an $AA$-stacked honeycomb
bilayer lattice is studied using the coupled cluster method implemented to
very high orders. On each monolayer the spins interact via nearest-neighbor
(NN) and frustrating next-nearest-neighbor isotropic antiferromagnetic
Heisenberg interactions with respective strength parameters $J_{1}>0$ and
$J_{2}\equiv\kappa J_{1}>0$. The two layers are coupled such that NN
interlayer pairs of spins also interact via a similar isotropic Heisenberg
interaction of strength $J_{1}^{\perp}\equiv\delta J_{1}$, which may be of
either sign. In particular, we locate with high accuracy the complete phase
boundaries in the $\kappa$-$\delta$ half-plane with $\kappa>0$ of the two
quasiclassical collinear antiferromagnetic phases with Néel or Néel-II
magnetic order in each monolayer, and the interlayer NN pairs of spins either
aligned (for $\delta<0$) or anti-aligned (for $\delta>0$) to one another.
Compared to the two-sublattice Néel order, in which all NN intralayer pairs of
spins are antiparallel to one another, the four-sublattice Néel-II order is
characterized by NN intralayer pairs of spins on the honeycomb lattice being
antiparallel to one another along zigzag (or sawtooth) chains in a specified
direction from among the three equivalent honeycomb-lattice directions, and
parallel to one another for the corresponding interchain pairs.

###### keywords:

honeycomb bilayer lattice , coupled cluster method , antiferromagnetism ,
regions of stability , collinear phases

††journal: Journal of magnetism and magnetic materials

##  1 Introduction

Frustrated quantum spin-lattice models in two dimensions, in which a variety
of different phases can emerge from zero-point quantum fluctuations even at
zero temperature ($T=0$), are prototypical models in which to study quantum
phase transitions as some system parameter is varied at $T=0$. The possible
phases obviously include quasiclassical states with definite magnetic long-
range order (LRO), such that SU(2) spin-rotational symmetry is broken and the
average local on-site magnetization (proportional to
$\langle\mathbf{s}_{j}\rangle$, where $\mathbf{s}_{j}$ is the spin operator on
lattice $j$) is nonzero. More interestingly, however, they also include
quantum paramagnetic phases with no classical counterparts, for which
$\langle\mathbf{s}_{j}\rangle=0$, and hence all magnetic LRO has melted.

Such quantum paramagnetic phases that preserve the SU(2) spin-rotational
symmetry may or may not still break other lattice symmetries. The former case
includes a variety of valence-bond crystalline (VBC) phases in which one or
more lattice symmetries are broken by the formation of spin-singlet states
involving static spin complexes arranged in some regular pattern. For example,
a plaquette VBC phase typically breaks only translational symmetry, while a
dimer VBC phase breaks both translational and rotational lattice symmetries.
By contrast, quantum spin liquid (QSL) phases [1] conserve all lattice
symmetries, and any other symmetries of the Hamiltonian, by definition. QSL
phases themselves may be either gapless or gapped, the latter typically with
some topological order (e.g., of the $Z_{2}$ type). They are archetypal
disordered ground states that possess massive many-body entanglement.

A typical generic system in which strong quantum fluctuations can combine with
spin frustration to produce magnetically disordered ground-state (GS) phases
is the $J_{1}$–$J_{2}$ Heisenberg model on a given two dimensional (2D)
lattice, in which both nearest-neighbor (NN) and next-nearest-neighbor (NNN)
Heisenberg exchange interactions, of respective strengths $J_{1}>0$ and
$J_{2}\equiv\kappa J_{1}>0$, compete with one another. By now there have been
very many calculations on a variety of 2D lattices (e.g., square, honeycomb,
triangular, kagome, and other Archimedean lattices) that show that each of
these $J_{1}$–$J_{2}$ spin-lattice systems can exhibit such quantum
paramagnetic GS phases, typically in the vicinity of the value for the
frustration parameter $\kappa$ that denotes maximal frustration in the
classical limit $s\to\infty$, where $s$ is the spin quantum number of each of
the lattice spins.

Generally, for a lattice of given dimensionality (here taken to be equal to
2), quantum fluctuations tend to be larger for lower values of both the spin
quantum number $s$ and the lattice coordination number $z$. For these reasons
the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$ Heisenberg magnet on the honeycomb
lattice occupies a special niche in the present context since it is the
simplest (in the sense of being composed of only one type of polygon, i.e.,
the hexagon) of the four 2D Archimedean lattices (all comprised of
arrangements of regular polygons, with every site equivalent to all others)
that take the lowest coordination, $z=3$. Each Archimedean lattice is uniquely
defined by specifying the ordered sequence of polygons that surround each of
the equivalent vertices. The other three Archimedean lattices with $z=3$ all
comprise lattices with more than one type of polygon. For example, of the
remaining two bipartite lattices with $z=3$, the CaVO ($4\cdot 8^{2}$) lattice
comprises squares and octagons, while the SHD ($4\cdot 6\cdot 12$) lattice
comprises squares, hexagons, and dodecagons. The non-bipartite Archimedean
lattice with $z=3$ is the star ($3\cdot 12^{2}$) lattice, which comprises
triangles and dodecagons. We note that of the four Archimedean lattices with
lowest coordinations, $z=3$, the honeycomb ($6^{3}$) lattice is special in
that it is the only one in which all of the edges are also equivalent. The
remaining three, namely the CaVO, SHD, and star lattices, all contain two
different types of NN bonds.

The honeycomb lattice is also unlike the other well-studied square,
triangular, and kagome lattices for the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$
model in that it is non-Bravais. Unlike these others it comprises two sites
per unit cell, with the structure of two interlacing triangular Bravais
sublattices. An immediate consequence is that the spin-$\frac{1}{2}$
$J_{1}$–$J_{2}$ model on the honeycomb lattice is free from the restrictions,
which apply to the corresponding model on Bravais lattices, that are imposed
by the Lieb-Schultz-Mattis (LSM) theorem [2] and its generalizations due to
Hastings [3] and others [4, 5, 6]. Broadly speaking, the extension by Hastings
of the LSM theorem implies that any such spin-lattice system that has a half-
odd-integer spin per unit cell cannot have both a unique ground state and a
gap in the excitation spectrum. Thus, for spin-$\frac{1}{2}$ quantum magnets
on Bravais lattices, such as the square, triangular, and kagome lattices, a
gapped spectrum strictly implies a degenerate GS phase. This could be, for
example, either a VBC phase caused by a lattice symmetry breaking, or a QSL
phase (e.g., of the $Z_{2}$ variety) caused by topological degeneracy. By
contrast, for spin-$\frac{1}{2}$ models on the honeycomb lattice, with its two
sites per unit cell, in principle one can have a gapped quantum paramagnetic
QSL GS phase that does not break any symmetry.

For all of the above reasons the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$ model on
the honeycomb-lattice monolayer has received an enormous amount of attention
in recent years, with a large number of theoretical techniques being applied
to it [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. Despite this intense activity, there
is still no overall consensus on the structure of its ($T=0$) quantum phase
diagram, particularly in the paramagnetic regime, as a function of the
frustration parameter, $\kappa\equiv J_{2}/J_{1}$, as we discuss in more
detail in Sec. 2. While almost all studies concur that the system retains Néel
magnetic LRO for sufficiently weak frustration, $\kappa<\kappa_{c_{1}}^{>}$,
there remain disagreements both over the nature of the stable GS phases for
$\kappa>\kappa_{c_{1}}^{>}$ and the precise numerical value of the critical
parameter $\kappa_{c_{1}}^{>}$ at which Néel order melts, although in the
latter regard most recent calculations using theoretical techniques of high
potential accuracy do give values for $\kappa_{c_{1}}^{>}$ in the approximate
range $0.19\lesssim\kappa_{c_{1}}^{>}\lesssim 0.23$.

If we restrict further discussion about these uncertainties to potentially
high-accuracy methods performed either at high orders in some well-defined
sequence of approximations or in large-scale numerical implementations, there
is no doubt that they are in large part due to the twin facts that almost all
such methods are biased in favor of some predetermined GS phase and/or are not
performed from the outset in the infinite-lattice limit in which we are
interested. In the latter regard, for example, such potentially high-accuracy
techniques as the exact diagonalization (ED) of finite-sized lattices
comprising $N$ spins and the density-matrix renormalization group (DMRG)
method always require some form of finite-size scaling to extrapolate to the
thermodynamic limit, $N\to\infty$. Such extrapolations have been shown
explicitly (see, e.g., Ref. [35]) to be capable of containing large
uncertainties. While this is particularly true in cases where theoretical
considerations provide little or no information on which rigorously to base
the extrapolation scheme, it can also even be the case when some such guidance
exists.

Within this context the coupled cluster method (CCM) [36, 37, 38, 39, 40, 41,
42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54] has come to occupy a
special role in recent years since it is one of the very few high-accuracy
methods of modern quantum many-body theory that is applied from the very
outset in the thermodynamic limit, $N\to\infty$, and hence for which any need
for finite-size scaling is always obviated. It is in large part for that
reason why it is the method we utilize in the present work. Furthermore, it is
now widely accepted that the CCM offers one of the most flexible, most widely
applicable, and most accurate at a given level of computational resource, of
all ab initio techniques that are available for dealing with a diverse range
of problems in microscopic quantum many-body theory. In particular, the CCM
has already been applied to many different spin-lattice systems in the broad
arena of quantum magnetism (see e.g., Refs. [17, 23, 24, 25, 50, 51, 54, 55,
56, 57, 58, 59, 60] and references contained therein), including both the
spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$ model on the honeycomb lattice [23, 25] of
interest here, as well as its various extensions, e.g., to the spin-1
counterpart [59], the case of anisotropic ($XXZ$) couplings [56], and the
isotropic $XY$ version of the model [55].

Given the still unresolved nature of the $T=0$ quantum phase diagram of the
spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$ model on the honeycomb-lattice monolayer,
it seems worthwhile to examine larger classes of models, to which they reduce
as a special case. One such direction that has already received much attention
(see, e.g., Refs. [19, 24] and references contained therein) is to remain with
the spin-$\frac{1}{2}$ honeycomb-lattice monolayer, but also now to include
next-next-nearest-neighbor Heisenberg interactions of strength $J_{3}$,
resulting in the so-called $J_{1}$–$J_{2}$–$J_{3}$ model. Another, potentially
even more revealing, extension is to the corresponding spin-$\frac{1}{2}$
$J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on a honeycomb bilayer, which is the
model of interest here. Each layer comprises an identical frustrated
$J_{1}$–$J_{2}$ system, and the two layers are now connected by NN Heisenberg
exchange bonds of strength $J_{1}^{\perp}\equiv\delta J_{1}$, with the layers
arranged in $AA$ stacking [i.e., with each site of one (horizontal) monolayer
placed vertically above its equivalent on the other].

Figure 1: The $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the bilayer honeycomb
lattice, showing (a) the two layers $A$ (red) and $B$ (blue), the nearest-
neighbor (NN) bonds ($J_{1}=$ —–; $J_{1}^{\perp}=$ \- - -) and the four sites
($1_{A},2_{A},1_{B},2_{B}$) of the unit cell; (b) the intralayer bonds
($J_{1}=$ —–; $J_{2}=$ \- - -) on each monolayer and the monolayer Néel state;
and (c) the triangular Bravais lattice vectors ${\bf a}$ and ${\bf b}$, and
one of the three equivalent monolayer Néel-II states. Sites ($1_{A},2_{B}$)
and ($2_{A},1_{B}$) on the two triangular lattices of each monolayer are shown
by filled and empty circles respectively, and the spins are represented by the
(green) arrows on the lattice sites. For both states on the bilayer, spins on
NN sites between the two layers are antiparallel (parallel) for $\delta>0$
$(\delta<0)$.

The $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model is particularly interesting since
the $J_{2}$ and $J_{1}^{\perp}$ bonds act in quite different ways to destroy
the Néel magnetic LRO that exists when the (antiferromagnetic) NN
$J_{1}\;(>0)$ bonds prevail. Thus, while the inclusion of (antiferromagnetic)
NNN $J_{2}\;(>0)$ bonds acts to frustrate the NN $J_{1}$ bonds by tending to
favor another form of magnetic LRO, the NN $J_{1}^{\perp}$ bonds of either
sign do not  directly frustrate the Néel order. However, they do  compete with
the $J_{1}$ bonds by promoting interlayer NN dimers to develop, which are
spin-singlet (spin-triplet) pairs when $J_{1}^{\perp}>0\;(J_{1}^{\perp}<0)$.
The inclusion of interlayer $J_{1}^{\perp}$ bonds thus promotes a competition
with the intralayer $J_{1}$ bonds between a phase with Néel magnetic LRO on
each monolayer and a nonclassical paramagnetic phase of the dimerized VBC
variety, wherein the (independent) dimers are now between NN interlayer pairs.
This effect is purely quantum-mechanical in origin, since in the classical
($s\to\infty$) limit the $J_{1}^{\perp}$ bonds have no effect at all on the
intralayer Néel ordering. A priori, the effect is expected to be maximal for
the case $s=\frac{1}{2}$ considered here.

While the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on a
honeycomb bilayer (with $J_{2}\equiv\kappa J_{1};J_{1}^{\perp}\equiv\delta
J_{1};J_{1}>0$) has received some attention in the last few years, almost all
of the effort has focused on calculating the boundary of the Néel phase in the
quadrant of the $\kappa$-$\delta$ plane with $\kappa>0$, $\delta>0$ [61, 62,
63, 64]. Nevertheless, the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$ model on a
honeycomb monolayer has also been shown to have another form of magnetic LRO
in the range $\kappa_{c_{2}}^{<}<\kappa<\kappa_{c_{2}}^{>}$ , where
$\kappa_{c_{2}}^{<}>\kappa_{c_{1}}^{>}$, which is not present in the classical
($s\to\infty$) version of the model except precisely at the point
$\kappa=\frac{1}{2}$ where it is degenerate with a class of spiral states.
This phase, the so-called Néel-II phase, discussed more fully in Sec. 2, and
which corresponds to phase IV in Ref. [9], is actually the $T=0$ GS phase in
the classical $J_{1}$–$J_{2}$–$J_{3}$ honeycomb-lattice model only for values
$J_{3}<0$ (when $J_{1}>0$) [7, 9], where, although it is then actually
degenerate with an infinite manifold of noncoplanar spin configurations [9],
it is then stabilized by quantum fluctuations [9] via the order-by-disorder
mechanism [65, 66].

It has been shown by several authors [19, 24] that for the spin-$\frac{1}{2}$
case of the honeycomb-lattice $J_{1}$–$J_{2}$–$J_{3}$ model the Néel-II phase
is actually stabilized over a finite region of the phase space with $J_{i}>0$,
for $i=1,2,3$ (and see, in particular, Fig. 2 in each of Refs. [19, 24]). For
this reason it is of considerable interest to investigate also the phase
boundary of the Néel-II phase for the spin-$\frac{1}{2}$
$J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the honeycomb bilayer, which to our
knowledge, has only been previously investigated [67] in the half-plane
$\delta>0$, as also for the Néel phase. Again, to the best of our knowledge,
there have been no prior investigations of the stability of both of these
collinear magnetic orderings (i.e., Néel and Néel-II) on each monolayer of the
bilayer model, in the case of ferromagnetic (FM) interlayer coupling
($\delta<0$), despite the fact that this case is particularly interesting in
its own right, for reasons that we now elaborate.

Thus, in the limit $\delta\to-\infty$, the system clearly simulates exactly a
spin-1 $J_{1}$–$J_{2}$ model on a honeycomb-lattice monolayer. While far fewer
theoretical studies exist for this model than for its spin-$\frac{1}{2}$
counterpart, two recent calculations, one using the density-matrix
renormalization group (DMRG) on cylindrical systems [68] and another using the
CCM employed here [59] both agree that each of the quasiclassical Néel and
Néel-II antiferromagnetic (AFM) phases also form stable GS phases in the
spin-1 honeycomb-lattice monolayer version of the $J_{1}$–$J_{2}$ model, just
as in the spin-$\frac{1}{2}$ case discussed above. Thus, we can use the
spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ bilayer model as a tool to
interpolate, as a function of the interlayer coupling $\delta$, between the
spin-$\frac{1}{2}$ ($\delta=0$) and spin-1 ($\delta=-\infty$) $J_{1}$–$J_{2}$
honeycomb-lattice monolayer models, with the possibility in so doing of
shedding more light on the collinear AFM phases of both models.

##  2 The model

The $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on a honeycomb-lattice bilayer,
shown in Fig. 1, has a Hamiltonian given by

$\displaystyle H$ $\displaystyle=J_{1}\sum_{{\langle
i,j\rangle},\alpha}\mathbf{s}_{i,\alpha}\cdot\mathbf{s}_{j,\alpha}+J_{2}\sum_{{\langle\langle
i,k\rangle\rangle},\alpha}\mathbf{s}_{i,\alpha}\cdot\mathbf{s}_{k,\alpha}$ (1)
$\displaystyle\quad+J_{1}^{\perp}\sum_{i}\mathbf{s}_{i,A}\cdot\mathbf{s}_{i,B}$
$\displaystyle\equiv J_{1}h(\kappa,\delta)\,;\quad\kappa\equiv
J_{2}/J_{1}\,,\quad\delta\equiv J_{1}^{\perp}/J_{1}\,,$

where the sums over $\langle i,j\rangle$ and $\langle\langle
i,k\rangle\rangle$ run, respectively, over all NN and NNN pairs of spins in
each (horizontal) monolayer, counting each pair once and once only, and where
the layer index $\alpha$ runs over the two monolayers, labelled $A$ and $B$.
The two monolayers are arranged in $AA$ stacking, such that each spin
$\mathbf{s}_{i,A}$ on layer $A$ lies vertically above its counterpart
$\mathbf{s}_{i,B}$ on layer $B$. The sites on each monolayer are arranged in a
regular honeycomb-lattice pattern, and each site $s_{i,\alpha}$ carries a
spin-$s$ particle described by the usual SU(2) spin operators
$\mathbf{s}_{i,\alpha}\equiv(s^{x}_{i,\alpha},s^{y}_{i,\alpha},s^{z}_{i,\alpha})$,
with $\mathbf{s}_{i,\alpha}^{2}=s(s+1)\mathbbm{1}$.

We shall restrict attention here to the extreme quantum case, $s=\frac{1}{2}$.
We further restrict ourselves to the most interesting case where both
intralayer Heisenberg bonds (i.e., NN bonds and NNN bonds) are AFM in nature
(i.e., $J_{1}>0$ and $J_{2}\equiv\kappa J_{1}>0$), such that they tend to
frustrate one another. The parameter $J_{1}$ then merely acts to set the
overall energy scale, and the Hamiltonian may thus be written as
$H=J_{1}h(\kappa,\delta)$, as in the last line of Eq. (1), thereby explicitly
demonstrating that the relevant model parameters are only $\kappa$ and
$\delta$. Our aim is thus to examine the $T=0$ quantum phase diagram of the
spin-$\frac{1}{2}$ version of the $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model
Hamiltonian of Eq. (1) on a honeycomb bilayer in the $\kappa$-$\delta$ half-
plane with $\kappa>0$. Although we also concentrate here on the hitherto
unexplored regime of FM interlayer coupling ($\delta<0$), we shall also
present some results with $\delta>0$ for completeness.

Our primary interest will be to investigate the regimes of stability in the
$\kappa$-$\delta$ plane of the two collinear AFM phases (i.e., the
quasiclassical Néel and Néel-II phases) on each monolayer as the interlayer
coupling parameter $\delta$ is introduced. The $AA$-stacked bilayer honeycomb-
lattice itself is shown in Fig. 1(a), and the spin arrangements of the Néel
and Néel-II phases on each monolayer are illustrated respectively in Figs.
1(b) and 1(c). Clearly, the Néel state has AFM sawtooth (or zigzag) chains
(i.e., with spins that alternate in direction along the chains) in each of the
three equivalent honeycomb-lattice directions. The Néel-II state, by contrast
has such AFM sawtooth chains along only one of the three directions, as shown
in Fig. 1(c), for example, along the $x$ direction. The NN spins between
adjacent such AFM sawtooth chains in the Néel-II state are aligned parallel to
one another. The pattern of such adjacent parallel spin pairs in the Néel-II
state is identical to that of the dimers in the staggered-dimer VBC (SDVBC)
state, such that the Néel-II and SDVBC state break the same lattice symmetries
(i.e., rotational and translational). For this reason many theoretical
investigations find it particularly difficult to differentiate these two
states in the $T=0$ quantum phase diagram of any such model that supports one
or other of them (e.g., Ref. [19]). The SDVBC itself is also sometimes known
as the lattice nematic state (see, e.g., Ref. [10]).

##  3 Methodology

Since the CCM formalism has been discussed extensively elsewhere in the
literature (see, e.g., Refs. [48, 49, 51] and references contained therein),
we confine ourselves here to a brief discussion of its most important features
for the problem at hand. The exact ket and bra GS wave functions are defined
to be $|\Psi\rangle$ and
$\langle\tilde{\Psi}|\equiv\langle\Psi|/\langle\Psi|\Psi\rangle$. The system
under study has a Hilbert space that we assume may be described in terms of a
normalised model (or reference) state $|\Phi\rangle$,
$\langle\Phi|\Phi\rangle=1$, which acts as a cyclic vector for a corresponding
set of mutually commuting  multiconfigurational creation operators
$\\{C_{I}^{+}\\}$, $[C_{I}^{+},C_{J}^{+}]=0$, such that $|\Phi\rangle$ is a
generalized vacuum state with respect to them,
$C_{I}^{-}|\Phi\rangle=0=\langle\Phi|C_{I}^{+}\,,\forall I\neq 0$, in a
notation where we define $C_{I}^{-}\equiv(C_{I}^{+})^{{\dagger}}$ and
$C_{0}^{+}\equiv\mathbbm{1}$. The GS ket $|\Psi\rangle$ is chosen within the
CCM to have the intermediate normalization, $\langle\Phi|\Psi\rangle=1$.

The index $I$ is a set index and, in general, $C_{I}^{+}$ comprises a product
of single-particle operators, as we describe more fully below in the current
context of spin-lattice models. The set $\\{I\\}$ is complete in the usual
sense that the set of states $\\{C_{I}^{+}|\Phi\rangle\\}$ provides a complete
basis for the ket Hilbert space. Although not vital, it is also convenient to
choose the basis to be orthonormalized, such that
$\langle\Phi|C_{I}^{-}C_{J}^{+}|\Phi\rangle=\delta_{IJ}$, where $\delta_{IJ}$
is an appropriately defined Kronecker symbol.

The GS wave functions $|\Psi\rangle$ and $\langle\tilde{\Psi}|$ are now
formally parametrized in the CCM, in a independent fashion, in terms of the
distinctive exponentiated forms involving the correlation operators,

$|\Psi\rangle=e^{S}|\Phi\rangle\,,\quad\langle\tilde{\Psi}|=\langle\Phi|\tilde{S}e^{-S}\,;$
(2)

$S=\sum_{I\neq 0}{\cal
S}_{I}C_{I}^{+}\,,\quad\tilde{S}=\mathbbm{1}+\sum_{I\neq 0}\tilde{\cal
S}_{I}C_{I}^{-}\,,$ (3)

which are a hallmark of the method. Since $S$ and $\tilde{S}$ are henceforth
treated as independent operators, the Hermiticity relation,
$\langle\tilde{\Psi}|=(|\Psi\rangle)^{{\dagger}}/\langle\Psi|\Psi\rangle$, may
be violated when subsequent approximations are made, typically by truncating
the complete set of configurations $\\{I\\}$ in the sums of Eq. (3), as
described more fully below. This possible shortcoming is far outweighed in
practice by the fact that the CCM parametrizations of Eqs. (2) and (3) always
exactly  satisfy the very important Hellmann-Feynman theorem at all  such
levels of approximation [49]. It is precisely this feature that guarantees the
robustness and accuracy of numerical results obtained within the CCM
framework.

The GS expectation value $\bar{X}$ of an arbitrary operator $X$ may thus be
expressed within the CCM purely in terms of the $c$-number correlation
coefficients $\\{{\cal S}_{I},\tilde{{\cal S}}_{I}\\}$ as

$\bar{X}=\bar{X}({\cal S}_{I},\tilde{{\cal
S}}_{I})\equiv\langle\Phi|\tilde{S}e^{-S}Xe^{S}|\Phi\rangle\,.$ (4)

The coefficients $\\{{\cal S}_{I},\tilde{{\cal S}_{I}}\\}$ are themselves
obtained, by minimizing the GS energy functional, $\bar{H}=\bar{H}({\cal
S}_{I},\tilde{{\cal S}}_{I})$, as obtained from Eq. (4) by replacing the
arbitrary operator $X$ by the Hamiltonian $H$, with respect to all members of
the set with $I\neq 0$ separately. The GS expectation value $\bar{X}$ (and
$\bar{H}$ in particular) is evaluated in practice by employing the nested
commutation expansion for the CCM similarity transform $e^{-S}Xe^{S}$ of the
operator $X$,

$e^{-S}Xe^{S}=\sum_{n=0}^{\infty}\frac{1}{n!}[X,S]_{n}\,,$ (5)

in terms of the $n$-fold nested commutators $[X,S]_{n}$, which are defined
iteratively as follows,

$[X,S]_{n}\equiv[[X,S]_{n-1},S],\quad[X,S]_{0}=X\,.$ (6)

The CCM parametrizations of Eqs. (2) and (3) are specifically chosen so that
the otherwise infinite sum in Eq. (5) will actually terminate exactly  at a
low finite order for any operator $X$ all of whose terms involve only a
product of a finite number of single-particle operators, such as the present
Hamiltonian of Eq. (1). The reason for this termination lies simply in the
twin facts that all components of $S$ in the expansion of Eq. (3) commute with
one another, and the basic single-particle operators form a Lie algebra,
together with the model state $|\Phi\rangle$ being a vacuum state for all
operators $C_{I}^{-}$. This same key feature of the CCM is also responsible
for all terms in the expansion of $\bar{H}$ from Eq. (4) being linked. It is
precisely this fact that implies immediately that the CCM automatically
satisfies exactly the Goldstone linked-cluster theorem at any level of
truncation in the expansions of Eq. (3) for the correlation operators. In turn
this also immediately guarantees the size-extensivity of the method at all
such levels of approximation. Thus, in the solution of the explicit CCM
equations for the basic amplitudes $\\{{\cal S}_{I},\tilde{{\cal S}}_{I}\\}$
that follow from the extremization of $\bar{H}$ (viz.,
$\delta\bar{H}/\delta\tilde{{\cal S}}_{I}=0=\delta\bar{H}/\delta{\cal
S}_{I},\forall I\neq 0$), the only approximation that is ever made in any CCM
calculation is to decide which configurations $\\{I\\}$ are to be retained in
the expansion of Eq. (3), as described for the problem at hand below.

A simple choice of CCM model state $|\Phi\rangle$ for any spin-lattice problem
is any quasiclassical state with perfect magnetic LRO. All such states are
specified uniquely by fixing the spin on each lattice site independently in
terms of its projection onto some given spin quantization axis. For the
present application we will thus use both of the Néel and Néel-II states on
each honeycomb-lattice monolayer as our CCM model states. In order to treat
each lattice site and each spin for any  such quasiclassical model state
$|\Phi\rangle$ as being completely equivalent to one another it is now useful
to choose a local spin quantization axis (or, equivalently, by making a
suitable passive rotation in spin space) on every lattice site independently
so that in such local axes the model state is one in which every spin points
in the same downwards direction (i.e., along the negative $z$ direction),
$|\Phi\rangle=|\\!\downarrow\downarrow\downarrow\cdots\downarrow\rangle$. Such
a description has the effect that all such spin-lattice calculations may, from
this point on, be treated on an equal footing and via a universal
computational code. All that then distinguishes one case from another is the
resulting spin Hamiltonian, which hence needs to be re-expressed in terms of
the local spin axes now specified separately for each model state.

The magnetic order parameter is now taken to be the sublattice magnetization
or average local on-site magnetization, $M$. In the local  rotated spin axes
defined above, this simply takes the universal form

$M=-\frac{1}{N}\sum_{i=1}^{N}\langle\Phi|\tilde{S}e^{-S}s_{l_{i}}^{z}e^{S}|\Phi\rangle\,,$
(7)

where the index $l_{i}\equiv(k_{i},\alpha)$ labels the sites on both bilayers
$\alpha=A,B$ of the system, and $N(\rightarrow\infty)$ is the total number of
spins on the bilayer.

The above choice of local spin axes now ensures that the multiconfigurational
creation operators $C_{I}^{+}$ may be chosen to be products of single spin-
raising operators $s_{l_{i}}^{+}\equiv s_{l_{i}}^{x}+is_{l_{i}}^{y}$. The set
indices $I$ hence simply become sets of bilayer lattice site indices,
$\\{{I}\\}\rightarrow\\{{l_{1},l_{2},\cdots,l_{n};n=1,2,\cdots,2sN}\\}$, in
which any given site index $l_{i}$ may be repeated, but so that it appears at
most $2s$ times in the most general case where each site carries a spin-$s$
particle. Thus, we have $C_{I}^{+}\rightarrow s_{l_{1}}^{+}s_{l_{2}}^{+}\cdots
s_{l_{n}}^{+}$, with $n=1,2,\cdots,2sN$. For the present paper we restrict
ourselves to the case $s=\frac{1}{2}$, for which quantum effects are expected
to be the largest, so that no lattice site in any multiconfigurational index
$I$ may appear more than once.

As we have already intimated above, the only  approximation that we now ever
make is to restrict the configurations $\\{{I}\\}$ that we retain in the
expansions of Eq. (3) for the CCM correlation operators, $S$ and $\tilde{S}$.
For the present calculations we shall utilize the well-tested and widely used
localized (lattice-animal-based subsystem) LSUB$n$ hierarchy of
approximations. At a specified $n$th order in the LSUB$n$ scheme one retains
all multispin-flip correlations corresponding to all possible configurations
of spins on the lattice that are confined to clusters of at most $n$
contiguous sites. A set $I$ of sites is defined to be contiguous in this sense
if every site of the set is a NN to at least one other site of the set, in
some specified geometry.

In order to reduce the number of independent  multispin-flip configurations at
a given LSUB$n$ order to the minumum, $N_{f}=N_{f}(n)$, at a given LSUB$n$
order, we utilize all of the space- and point-group symmetries of both the
system Hamiltonian and the particular CCM model state being used. In the same
vein we also employ any relevant conservation laws. For example, both the
quasiclassical Néel and Néel-II AFM states are eigenstates of the operator
$s_{T}^{z}\equiv\sum_{k,\alpha}s_{k,\alpha}^{z}$, using global  spin axes,
which represents the total $z$ component of spin for the system as a whole.
Both states have corresponding eigenvalue $s_{T}^{z}=0$, and accordingly all
of the multispin-flip configurations $I$ retained in the expansions of the CCM
correlation operators $S$ and $\tilde{S}$ given by Eq. (3) are chosen to
follow this conservation law.

Even after all such symmetries and conservation laws have been fully
incorporated the minimum number of independent LSUB$n$ configurations,
$N_{f}(n)$, grows rapidly with the truncation index $n$, typically
exponentially as $n$ becomes large. The available computer power then
typically determines the maximum order $n$ that can be calculated in practice
for a given model. In the present case, for the spin-$\frac{1}{2}$
$J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the $AA$-stacked honeycomb bilayer
lattice, by employing both massive parallelization and large-scale
supercomputing resources, as well as a purpose-built and customized computer
algebra package [69] to both derive and solve [50] the large sets of $N_{f}$
coupled CCM equations for the GS ket- and bra-state coefficients $\\{{{\cal
S}_{I},\tilde{{\cal S}}_{I}}\\}$, we can perform LSUB$n$ calculations up to
the very high order $n=10$. Thus, for the model under consideration, when the
two monolayers are coupled ferromagnetically (i.e., when $\delta<0$), we have
$N_{f}(10)=64\,780\;(183\,939)$ for the GS properties when the CCM model state
is chosen so that each monolayer has Néel (Néel-II) AFM order. The
corresponding numbers for the case when the two monolayers are coupled
antiferromagnetically (i.e., when $\delta>0$) are
$N_{f}(10)=70\,118\;(197\,756)$.

Figure 2: CCM results for the GS energy per spin $E/N$ (in units of $J_{1}$)
versus the intralayer frustration parameter, $\kappa\equiv J_{2}/J_{1}$, for
the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the bilayer
honeycomb lattice (with $J_{1}>0$), for a selected value ($\delta=-1.2$) of
the scaled interlayer exchange coupling constant, $\delta\equiv
J_{1}^{\perp}/J_{1}$. Results based on (a) the Néel state and (b) the Néel-II
state on each monolayer, and the two layers coupled so that NN spins between
them are parallel to one another, as CCM model states are shown in LSUB$n$
approximations with $n=2,4,6,8,10$, together with the corresponding
LSUB$\infty$ extrapolated result based on Eq. (8) and the LSUB$n$ data sets
$n=\\{2,6,10\\}$. Subfigure (c) shows an exploded view of the intermediate
region.

It is clear that our LSUB$n$ approximation becomes asymptotically exact as the
truncation index grows, $n\to\infty$. Since no approximations are ever made in
calculating any given LSUB$n$ approximant for any GS observable, it then
follows that the only source of errors in our CCM calculations is in the very
last step where the sequences of LSUB$n$ approximants are themselves
extrapolated to the exact, $n\to\infty$, limit. A great deal of empirical
evidence now exists on how to perform such extrapolations, based on a large
corpus of applications to diverse spin-lattice systems. For example, a well
tested and highly accurate extrapolation scheme for the GS energy per spin,
$E/N$, and see, e.g., Refs. [70, 71, 72, 51, 73, 74, 75, 76, 77, 78, 79, 17,
80, 23, 81, 24, 82, 83, 25, 84, 57, 58, 85, 59]) is given by

$\frac{E(n)}{N}=e_{0}+e_{1}n^{-2}+e_{2}n^{-4}\,.$ (8)

By fitting our LSUB$n$ approximants $E(n)/N$ to Eq. (8) we may thereby obtain
the extrapolated (exact LSUB$\infty$) value $e_{0}$.

The GS expectation values of other physical observables typically converge
slower than the energy, as is wholly to be expected. Thus, for example, the
LSUB$n$ approximants $M(n)$ to the magnetic order parameter of strongly
frustrated systems, particularly for those with a quantum phase transition
between states with and without magnetic LRO, such as our present model, have
been found (and see, e.g., Refs. [86, 87, 88, 89, 90, 53, 74, 75, 67]) to be
accurately fitted by the well tested extrapolation scheme

$M(n)=\mu_{0}+\mu_{1}n^{-1/2}+\mu_{2}n^{-3/2}\,.$ (9)

Once again, the extrapolated (exact LSUB$\infty$) value $\mu_{0}$ for $M$ may
then be obtained by fitting a sequence of LSUB$n$ approximants $M(n)$ to Eq.
(9).

When making use of extrapolation schemes such as those in Eqs. (8) and (9) one
needs to be aware of any “staggering effects” that may be present in the
approximant sequences. For example, it is well known in perturbation theory
that there usually exists an odd/even or $(2m-1)/2m$ (where
$m\in\mathbb{Z}^{+}$ is a positive integer) staggering effect in the sequence
of $n$th order approximants for various physical observables. In such cases,
where we often know exact  extrapolation schemes, both the $n=(2m-1)$ and
$n=2m$ subsequences obey an extrapolation scheme of the same sort (i.e., with
identical exponents in the leading and sub-leading terms), but where the
respective coefficients of the terms other than that corresponding to the
extrapolated $(n\to\infty)$ value itself may differ. In such cases one should
clearly not mix odd- and even-order terms together in a single extrapolation
sequence unless the staggering is separately incorporated properly. In general
it is difficult to include the staggering explicitly in a robust fashion, and
one then usually separately extrapolates the odd terms and the even terms. Our
CCM SUB$n$ sequences of approximants for all physical observables also display
such an odd/even staggering to a greater or lesser degree, dependant on both
the model and the observable. Since the Hamiltonian of Eq. (1) is itself
bilinear in the single-spin operators, it is then natural to confine ourselves
to extrapolating the even-order $(n=2m)$ approximants for any physical
observable in any CCM spin-lattice calculation for the model.

While the above $(2m-1)/2m$ staggering in LSUB$n$ sequences of approximants
for any physical quantity is common to all spin-lattice models on all
lattices, honeycomb lattices tend to exhibit an additional  $(4m-2)/4m$
staggering in the even subsequences, as has been noted elsewhere [23, 25, 59,
67, 91], and which now seems to originate in the non-Bravais nature of the
honeycomb lattice [67], which itself comprises two interlacing triangular
Bravais sublattices. Each of these displays the above-mentioned $(2m-1)/2m$
staggering, and the composite honeycomb lattice then magnifies the effect
twofold into the observed $(4m-2)/4m$ staggering of the even $(n=2m)$
subsequence and, presumably, also a $(4m-3)/(4m-1)$ staggering of the odd
$(n=2m-1)$ subsequence. Since, as we have noted above, we are restricted for
the present model to LSUB$n$ approximants with $n\leq 10$, in order to take
all of these staggering effects into account, we restrict all extrapolations
such as those based on Eqs. (8) and (9) to LSUB$n$ data sets with
$n=\\{2,6,10\\}$ for the results that we present in Sec. 4.

##  4 Results

Figure 3: CCM results for the GS magnetic order parameter $M$ versus the
intralayer frustration parameter, $\kappa\equiv J_{2}/J_{1}$, for the
spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the bilayer
honeycomb lattice (with $J_{1}>0$), for a selected value ($\delta=-1.2$) of
the scaled interlayer exchange coupling constant, $\delta\equiv
J_{1}^{\perp}/J_{1}$. Results based on (a) the Néel state and (b) the Néel-II
state on each monolayer, and the two layers coupled so that NN spins between
them are parallel to one another, as CCM model states are shown in LSUB$n$
approximations with $n=2,4,6,8,10$, together with the corresponding
LSUB$\infty$ extrapolated result based on Eq. (9) and the LSUB$n$ data sets
$n=\\{2,6,10\\}$.

We first present results for the GS energy per spin (in units of $J_{1}$,
$E/(NJ_{1})$, for our bilayer system, based separately on each of the Néel and
Néel-II states on each monolayer as our CCM model state, in the regime of FM
interlayer coupling, so that in the overall model state the two layers are
coupled such that NN interlayer spin pairs are parallel to one another. We
show in Fig. 2 results at a typical FM value, $\delta=-1.2$, of the interlayer
coupling parameter as a function of the intralayer frustration parameter
$\kappa$. One clearly observes the very rapid convergence of the CCM LSUB$n$
sequences in each case. We also show the corresponding LSUB$\infty$
extrapolated results, $e_{0}/J_{1}$, which are obtained from Eq. (8) and the
corresponding LSUB$n$ data sets with $n=\\{2,6,10\\}$ as input to circumvent
the $(4m-2)/4m$ staggering effect discussed above and which we will discuss
below in more detail in connection with the corresponding results for the
magnetic order parameter $M$, where we will also consider the natural
transition points of the LSUB$n$ results that can be observed, for example, in
Fig. 2(c). We note from Fig. 2 the extremely rapid convergence of the set of
LSUB$n$ approximants to the GS energy, as the truncation index $n$ is
increased, for both phases shown. The results presented in Fig. 2(c) also
clearly foreshadow the existence of an intermediate GS (non-magnetic) phase
between the Néel and Néel-II bilayer phases for this model.

We turn next to our results for the magnetic order parameter $M$. Thus, we
show in Figs. 3(a) and 3(b) results at the same value $\delta=-1.2$ of the
interlayer coupling as those shown in Figs. 2(a) and 2(b), respectively, for
the GS energy. We note immediately that, as is true in all practical
implementations of the CCM, the LSUB$n$ approximants $M(n)$ extend beyond the
actual (LSUB$\infty$) transition point(s) for the phase in question, for all
finite values of the truncation index $n$, out to some corresponding
termination point(s), beyond which no real solution exists for the respective
LSUB$n$ equations. Furthermore, these LSUB$n$ termination points clearly
appear to converge uniformly to the LSUB$\infty$ estimates of the
corresponding quantum critical points (QCPs) as $n$ is increased, as has been
observed many times before (and see, e.g., Refs. [23, 25, 51, 80]). The
$(4m-2)/4m$ staggering effect is clearly observed in the LSUB$n$ results for
$M(n)$, and is particularly marked for the Néel-II results shown in Fig. 3(b).
Thus, the corresponding LSUB$\infty$ results $\mu_{0}$ for $M$ shown in Fig. 3
are obtained from fitting Eq. (9) to the respective LSUB$n$ data sets with
$n=\\{2,6,10\\}$ only.

We note from Fig. 3 that at all orders in our LSUB$n$ expansions of the CCM
cluster operators the magnetic order parameter $M$ is overestimated, but as
the truncation index $n$ is increased $M$ decreases. This is fully as
expected, since as one moves to higher values of the truncation index $n$ one
is including more correlations. These correlations simply involve clusters of
“wrong” spins with respect to the corresponding perfectly ordered (Néel or
Néel-II) classical state. Clearly, the effect of systematically incorporating
more and more such clusters of “wrong” spins will thus lead naturally to a
systematic reduction in the order parameter, precisely as is observed in Fig.
3 for both ordered antiferromagnetic phases.

Results such as these shown in Fig. 3(a) clearly indicate the existence of an
upper QCP, $\kappa^{>}_{c_{1}}=\kappa^{>}_{c_{1}}(\delta)$, beyond which the
Néel phase ceases to be a stable GS phase. Thus, for example, we see from Fig.
3(a) that $\kappa_{c_{1}}^{>}(\delta=-1.2)\approx 0.261$, which may be
compared with the corresponding CCM result for the $J_{1}$–$J_{2}$ honeycomb-
lattice monolayer, $\kappa_{c_{1}}^{>}(\delta=0)\approx 0.183$, using the same
LSUB$n$ data set with $n=\\{2,6,10\\}$ as input to Eq. (9) [63]. Similarly,
results such as those shown in Fig. 3(b) also indicate both a lower QCP,
$\kappa_{c_{2}}^{<}=\kappa_{c_{2}}^{<}(\delta)$, and an upper QCP,
$\kappa^{>}_{c_{2}}=\kappa_{c_{2}}^{>}(\delta)$, such that the Néel-II phase
is a stable GS phase only in the regime
$\kappa_{c_{2}}^{<}(\delta)<\kappa<\kappa_{c_{2}}^{>}(\delta)$ for a given
value of the scaled interlayer coupling parameter $\delta$. Thus, from Fig.
3(b) we see that $\kappa^{<}_{c_{2}}(\delta=-1.2)\approx 0.302$ and
$\kappa_{c_{2}}^{>}(\delta=-1.2)\approx 3.563$, which values may be compared
with the corresponding CCM results for the $J_{1}$–$J_{2}$ honeycomb-lattice
monolayer, $\kappa_{c_{2}}^{<}(\delta=0)\approx 0.449$ and
$\kappa_{c_{2}}^{>}\approx 1.487$, again using the same LSUB$n$ data sets with
$n=\\{2,6,10\\}$ as input to Eq. (9) [67].

Figure 4: CCM results for the GS magnetic order parameter $M$ versus the
intralayer frustration parameter, $\kappa\equiv J_{2}/J_{1}$, for the
spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the bilayer
honeycomb lattice (with $J_{1}>0$), for a variety of ferromagnetic values of
the scaled interlayer exchange coupling constant, $\delta\equiv
J_{1}^{\perp}/J_{1}$. In each case we show extrapolated results, based on the
Néel state as CCM model state on each monolayer, and the two layers coupled so
that NN spins between them are parallel to one another, obtained from using
Eq. (9) with the corresponding LSUB$n$ data sets $n=\\{2,6,10\\}$. See the
text for an explanation of the curve labelled $\delta=-\infty$ $(s=1)$.

In Fig. 4 we show a set of LSUB$\infty$ extrapolated curves for the magnetic
order parameter $M=M(\kappa)$ in the Néel GS phase, analogous to that shown in
Fig. 3(a) for the specific value $\delta=-1.2$, but now for a variety of
values of $\delta$ in the regime where the two layers are coupled
ferromagnetically (so that NN interlayer pairs align parallel to one another),
i.e., for $\delta<0$. In each case we plot the value $\mu_{0}$ obtained from
fitting Eq. (9) to the corresponding LSUB$n$ data for $M(n)$ with
$n=\\{2,6,10\\}$. We also show in Fig. 4 the corresponding CCM result for the
Néel GS phase of the spin-1 $J_{1}$–$J_{2}$ model on a honeycomb-lattice
monolayer [59]. We note that in this case, instead of the LSUB$n$
approximation scheme, the alternative so-called SUB$n$-$n$ scheme has been
used, as we now explain.

The CCM SUB$n$-$m$ approximation retains all multispin configurations $I$ in
Eq. (3) for the CCM correlation operators $S$ and $\tilde{S}$ that involve $n$
or fewer spin flips (with respect to the model state $|\Phi\rangle$) spanning
a range of no more than $m$ contiguous sites. In this context a single spin
flip requires the action of the spin-raising operator $s^{+}_{k}$ acting once.
Clearly, for the spin-$\frac{1}{2}$ systems, the two truncation schemes
LSUB$n$ and SUB$n$-$n$ are identical, since no more than one spin flip per
site is possible. However, for general spin-$s$ systems, up to $2s$ spin flips
per site are possible, so that LSUB$n$ $\equiv$ SUB$2sn$-$n$. Thus, for the
spin-1 system, the SUB$n$-$n$ truncation contains fewer configurations than
its LSUB$n$ counterpart for general values of the truncation index $n$.

Since the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on a
honeycomb-lattice bilayer in the limit $\delta\rightarrow-\infty$ clearly
corresponds exactly to a spin-1 $J_{1}$–$J_{2}$ model on a honeycomb-lattice
monolayer, the spin-1 monolayer result shown in Fig. 4 should provide a good
limiting-case check of the present model. Thus, we see that
$\kappa_{c_{1}}^{>}(\delta\rightarrow-\infty)\approx 0.248$ from the spin-1
model using the SUB$n$-$n$ data set with $n=\\{2,6,10\\}$, which may be
compared, for example, with the value $\kappa_{c_{1}}^{>}(\delta=-2.5)\approx
0.271$, as shown in Fig. 4, and which is also obtained from the LSUB$n$ data
set with $n=\\{2,6,10\\}$. The closeness of the two values clearly
demonstrates both the inherent accuracy of the CCM and the compatibility of
the two different approximation schemes, LSUB$n$ and SUB$n$-$n$.

From curves such as those shown in Fig. 4, from which we determined the values
$\kappa_{c_{1}}^{>}(\delta)$ for various values of the scaled interlayer
coupling parameter $\delta$, we can now plot the region in the
$\kappa$-$\delta$ plane where the Néel state on each monolayer forms the
stable GS phase, as shown in Fig. 5.

Figure 5: $T=0$ Néel phase diagram of the spin-$\frac{1}{2}$
$J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the bilayer honeycomb lattice with
$J_{1}>0$, $\delta\equiv J_{1}^{\perp}/J_{1}$, and $\kappa\equiv J_{2}/J_{1}$.
The darker (blue) region is the quasiclassical GS phase with AFM Néel order in
each monolayer (with the two layers coupled so that NN spins across them are
aligned for $\delta<0$ and anti-aligned for $\delta>0$), while in the lighter
(grey) region Néel order is absent. The filled and empty square symbols are
points at which the extrapolated GS magnetic order parameter $M$ for the Néel
phase vanishes, for specified values of $\delta$ and $\kappa$, respectively.
The Néel state on each monolayer is used as CCM model state, with the two
layers coupled so that NN spins between them are parallel (antiparallel) to
one another when $\delta<0\;(\delta>0)$, and Eq. (9) is used for the
extrapolations with the corresponding LSUB$n$ data sets $n=\\{2,6,10\\}$. The
vertical line at $\kappa=\kappa_{c_{1}}^{>}(\delta\rightarrow-\infty)\approx
0.248$ is taken from the $s=1$ honeycomb lattice monolayer result for
$\kappa_{c_{1}}^{>}$, as explained in the text.

To the present results obtained for the case $\delta<0$, where the CCM model
state has interlayer NN pairs aligned parallel to one another, we also display
corresponding results obtained for the case $\delta>0$, for which the CCM
model state has the same pairs aligned antiparallel to one another [63].
Perhaps the most noteworthy feature of Fig. 5 is the extremely sharp cusp in
the phase diagram for the Néel phase, centered at $\delta=0$. This feature
clearly explains the great sensitivity that is observed in practice in most
theoretical calculations of $\kappa_{c_{1}}^{>}(\delta=0)$, corresponding to
the QCP for the vanishing of Néel order in the spin-$\frac{1}{2}$
$J_{1}$–$J_{2}$ model on a honeycomb-lattice monolayer.

Figure 6: CCM results for the GS magnetic order parameter $M$ versus the
intralayer frustration parameter, $\kappa\equiv J_{2}/J_{1}$, for the
spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the bilayer
honeycomb lattice (with $J_{1}>0$), for a variety of ferromagnetic values of
the scaled interlayer exchange coupling constant, $\delta\equiv
J_{1}^{\perp}/J_{1}$. In each case we show extrapolated results, based on the
Néel-II state as CCM model state on each monolayer, and the two layers coupled
so that NN spins between them are parallel to one another, obtained from using
Eq. (9) with the corresponding LSUB$n$ data sets $n=\\{2,6,10\\}$. See the
text for an explanation of the curve labeled $\delta=-\infty$ $(s=1$). Figure
7: $T=0$ Néel-II phase diagram of the spin-$\frac{1}{2}$
$J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the bilayer honeycomb lattice with
$J_{1}>0$, $\delta\equiv J_{1}^{\perp}/J_{1}$, and $\kappa\equiv J_{2}/J_{1}$.
The darker (pink) region is the quasiclassical GS phase with AFM Néel-II order
in each monolayer (with the two layers coupled so that NN spins across them
are aligned for $\delta<0$ and anti-aligned for $\delta>0$), while in the
lighter (grey) region Néel-II order is absent. The filled and empty square
symbols are points at which the extrapolated GS magnetic order parameter $M$
for the Néel phase vanishes, for specified values of $\delta$ and $\kappa$,
respectively. The Néel-II state on each monolayer is used as CCM model state,
with the two layers coupled so that NN spins between them are parallel
(antiparallel) to one another when $\delta<0$ ($\delta>0$), and Eq. (9) is
used for the extrapolations with the corresponding LSUB$n$ data sets
$n=\\{2,6,10\\}$. The vertical lines at
$\kappa=\kappa_{c_{2}}^{<}(\delta\rightarrow-\infty)\approx 0.343$ and
$\kappa=\kappa_{c_{2}}^{>}(\delta=-\infty)\approx 1.274$ are taken from the
$s=1$ honeycomb-lattice monolayer results for $\kappa_{c_{2}}^{<}$ and
$\kappa_{c_{2}}^{>}$, respectively, as explained in the text.

Similar to Fig. 4 for the Néel phase, we show in Fig. 6 a set of LSUB$\infty$
extrapolated curves for the magnetic order parameter $M=M(\kappa)$ in the
Néel-II GS phase, analogous to the one shown in Fig. 3(b) for the particular
value $\delta=-1.2$. We again display results for various values of
$\delta(<0)$ in the regime where the two layers are coupled ferromagnetically
(i.e., with NN interlayer pairs aligned parallel to each other. Once again we
plot the value $\mu_{0}$ obtained from fitting Eq. (9) to the corresponding
data set for $M(n)$ with $n=\\{2,6,10\\}$. As in Fig. 4 for the Néel phase, we
also shown in Fig. 6 the corresponding CCM result for the Néel-II GS phase of
the spin-1 $J_{1}$–$J_{2}$ model on a honeycomb-lattice monolayer, where we
have again used Eq. (9) to extract the extrapolated value $\mu_{0}$ from fits
to the SUB$n$-$n$ data set with $n=\\{2,6,10\\}$. In this case we have
extended the results given in Ref. [59], which were presented only for the
range $0\leq\kappa\leq 1$ there (and see Fig. 5 of Ref. [59]), so as to also
extract the value $\kappa_{c_{2}}^{>}$ for the model. We find that, for this
model, $\kappa_{c_{2}}^{<}\approx 0.343$ and $\kappa_{c_{2}}^{>}\approx
1.274$, which values should now also provide good estimates for the limiting
values, $\kappa_{c_{2}}^{<}(\delta\rightarrow-\infty$) and
$\kappa_{c_{2}}^{>}(\delta\rightarrow-\infty)$, of the current
spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ honeycomb-lattice bilayer
model. For the purposes of comparison we quote the corresponding results,
$\kappa_{c_{2}}^{<}(\delta=-30)\approx 0.335$ and
$\kappa_{c_{2}}^{>}(\delta=-30)\approx 1.763$ for $\delta=-30$, and
$\kappa_{c_{2}}^{<}(\delta=-100)\approx 0.337$ and
$\kappa_{c_{2}}^{>}(\delta=-100)\approx 1.643$ for $\delta=-100$. Clearly, the
totally independent spin-1 monolayer results are in very good accord with the
current spin-$\frac{1}{2}$ bilayer results.

Using curves such as those shown in Fig. 6 to extract the values
$\kappa_{c_{2}}^{<}(\delta)$ and $\kappa_{c_{2}}^{>}(\delta)$ for a variety of
values of $\delta$, we can now also plot the region in the $\kappa$-$\delta$
plane where the Néel-II state on each monolayer forms the stable GS phase. The
result is shown in Fig. 7 where, in addition to the present results obtained
for the case $\delta<0$, where the CCM model has interlayer NN spins aligned
parallel to each other, we also display corresponding results obtained for the
case $\delta>0$, for which the CCM model state has the same interlayer NN
pairs aligned antiparallel to each other [67]. Just as in Fig. 4 for the Néel
phase, so the phase diagram in Fig. 7 for the Néel-II phase exhibits two sharp
cusps centered at $\delta=0$. These explain the sensitivity in providing
accurate estimates for $\kappa_{c_{2}}^{<}(\delta=0)$ and
$\kappa_{c_{2}}^{>}(\delta=0)$, which correspond to the QCPs for the
disappearance of Néel-II AFM order in the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$
model on a honeycomb-lattice monolayer.

Figure 8: $T=0$ phase diagram of the spin-$\frac{1}{2}$
$J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the bilayer honeycomb lattice with
$J_{1}>0$, $\delta\equiv J_{1}^{\perp}/J_{1}$, and $\kappa\equiv J_{2}/J_{1}$.
The darker (blue and pink) shaded regions are the quasiclassical GS phases
with AFM Néel and Néel-II orders in each monolayer, respectively (with the two
layers coupled so that NN spins across them are aligned for $\delta<0$ and
anti-aligned for $\delta>0$), while in the lighter (grey) shaded region
quasiclassical collinear order is absent. The filled and empty squares are
points at which the extrapolated GS magnetic order parameter for the two
quasicallical AFM phases vanishes, for specified values of $\delta$ and
$\kappa$, respectively. In each case the Néel or Néel-II state on each layer
is used as CCM model state, with the two layers coupled so that NN spins
between them are parallel (antiparallel) to one another when $\delta<0$
($\delta>0$), and Eq. (9) is used for the extrapolations with the
corresponding LSUB$n$ data sets $n=\\{2,6,10\\}$. The vertical lines at
$\kappa=\kappa_{c_{1}}^{>}(\delta=-\infty)\approx 0.248$,
$\kappa=\kappa_{c_{2}}^{<}(\delta=-\infty)\approx 0.343$, and
$\kappa=\kappa_{c_{2}}^{>}(\delta=-\infty)\approx 1.274$ are taken from the
$s=1$ honeycomb-lattice monolayer results for $\kappa_{c_{1}}^{>}$,
$\kappa_{c_{2}}^{<}$, and $\kappa_{c_{2}}^{>}$, respectively, as explained in
the text.

Finally, in Fig. 8 we combine all of our results to show the complete $T=0$
quasiclassical (Néel and Néel-II) phase diagram of our bilayer model. The
darker (blue and pink) shaded regions show, respectively, where the stable GS
phase has AFM Néel and Néel-II order on each monolayer, while in the lighter
(grey) shaded region all (quasiclassical) collinear magnetic order is absent
in each monolayer.

##  5 Discussion

We have implemented the powerful CCM technique of microscopic quantum many-
body theory (QMBT) to very high orders in a systematic (LSUB$n$) hierarchy of
approximations, to investigate the regions of stability of two quasiclassical
collinear magnetically-ordered phases of a frustrated spin-$\frac{1}{2}$
$J_{1}$–$J_{2}$–$J_{1}^{\perp}$ Heisenberg magnet on a $AA$-stacked honeycomb-
lattice bilayer, in the case where the intralayer NN and NNN couplings are
both antiferromagnetic in nature (i.e., $J_{1}>0$ and $J_{2}>0$), for all
values of the interlayer NN coupling parameter $J_{1}^{\perp}$. The dual
features of the CCM that it satisfies both the Goldstone linked cluster
theorem and the Hellmann-Feynman theorem at all levels of approximation are
unmatched by almost all alternative techniques of ab initio  QMBT, and thus
make it an ideal method for our purposes. An immediate consequence of
satisfying the former theorem is that we have been able to perform all
calculations from the outset in the thermodynamic (infinite-lattice,
$N\rightarrow\infty$) limit. As a result we have never had to resort to any
finite-size scaling of our numerical results, as is required in almost all
other high-accuracy techniques, and which is then often the major cause of
uncertainty in their results.

Our sole approximation has thus been to extrapolate numerically the sequences
of LSUB$n$ approximants for the GS energy and GS magnetic order parameter to
the $n\rightarrow\infty$ limit, where the CCM becomes exact. By implementing
the method computationally with the purpose-built computer-algebra (CCCM)
package [69] to very high orders (viz., with $n\leq 10$), we have been able to
perform the extrapolations with very high accuracy. In order to do so we have
circumvented not only the well-known $(2m-1)/2m$ staggering effect (where
$m\in\mathbb{Z}^{+}$ is a positive integer) in the sequences of $n$th order
approximants for various physical observables, which is always present, but
also the additional (and less well known) $(4m-2)/4m$ staggering effect of the
even $(n=2m)$ subsequences, which the honeycomb lattice also exhibits as a
probable consequence of its non-Bravais nature, by performing all LSUB$n$
extrapolations with the limited data sets $n=\\{2,6,10\\}$.

Given that we have fitted our LSUB$n$ results for the magnetic order parameter
$M$ with the data set $n=\\{2,6,10\\}$ to the three-parameter fit of Eq. (9),
it is not possible in this case to give a precise estimate of any systematic
errors in our results. Nevertheless, many prior applications of the
methodology have shown that the extrapolation scheme is very robust. In
general, in cases where we have more LSUB$n$ data points than three, a least-
squares fit gives a reliable estimate of the errors. Many such estimates have
been given in prior publications involving the CCM, particularly for
applications on comparable models on Bravais lattices for which only the
$(2m-1)/2m$ staggering effect is present, and in which the $(4m-2)/4m$
staggering effect is absent, and hence for which more LSUB$n$ data points are
available to perform the extrapolations. A good example can be found in Ref.
[60] where results have been presented for the analogous model to that
presented here, but on a square-lattice bilayer. The high accuracy of the
extrapolations demonstrated explicitly in that case gives us considerable
confidence in the location of the phase boundaries in our final phase diagram
of Fig. 8.

We note that for the case of the honeycomb-lattice monolayer (i.e., with
$\delta=0$), Néel-II order does not exist classically as the stable GS of the
$J_{1}$–$J_{2}$ model, except precisely at the point $\kappa=\frac{1}{2}$,
which is a point of maximal GS degeneracy. Thus, for the classical
$J_{1}$–$J_{2}$ model on the honeycomb lattice (with $J_{1}>0$, $J_{2}>0$)
Néel order exists for all values $\kappa\leq\frac{1}{6}$. Conversely, for all
values $\kappa>\frac{1}{6}$ the model has an infinitely degenerate family of
incommensurate, coplanar, spiral configurations of spins, described by an
ordering wave vector $\mathbf{Q}$, for which the direction is arbitrary. For
values of $\kappa$ in the range $\frac{1}{6}<\kappa<\frac{1}{2}$ these
classically degenerate solutions for $\mathbf{Q}$ form a closed contour around
the center point $\mathbf{Q}=\mathbf{\Gamma}\equiv(0,0)$ of the hexagonal
first Brillouin zone. By contrast, when $\kappa>\frac{1}{2}$, the solutions
for $\mathbf{Q}$ lie on pairs of closed contours that are now centered on any
two of the inequivalent corners of the first Brillouin zone. The point
$\kappa=\frac{1}{2}$ marks a second classical transition point between two
different forms of GS spiral phases, and it is precisely at this point that
both are also degenerate with the collinear Néel-II phase.

We have seen from Fig. 6 that for the case $\delta=0$ of the
spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the bilayer
honeycomb lattice, quantum fluctuations have now stabilized Néel-II order over
the range $\kappa_{c_{2}}^{<}(\delta=0)\approx
0.45<\kappa<\kappa_{c_{2}}^{>}(\delta=0)\approx 1.49$, in keeping with the
general observation that quantum fluctuations tend to favor collinear over
non-collinear forms of order. Nevertheless, the Néel-II order present is quite
fragile, as may be seen from the rather small values of the order parameter
$M$ in the $\delta=0$ curve in Fig. 6 over the entire region of its existence
(i.e., when $M>0$). By contrast, Fig. 6 also shows that even a small
interlayer coupling enhances the stability of the Néel-II phase considerably.

Although it is out of the scope of the present work to study in detail the
nature of the GS phases of the model in the lighter (grey) shaded areas of
Fig. 8, where collinear magnetic ordering is absent, nevertheless it is
interesting to speculate. Thus, for example, from our discussion above, we are
led to expect that in the region around the $\delta=0$ axis, and between the
two lobes of the region of stability where Néel-II ordering is present on each
monolayer, the GS phase will have spiral ordering. This assertion has been
lent credence, again by the use of the CCM, but now to study excited states
also, by the finding [67] that (at least part of) this region is gapless. From
the same CCM study of the excitation energy to the lowest excited state of the
model, the clear presence of a gapped state has also been shown to exist over
much of the paramagnetic region in Fig. 8 between the two islands of stability
(shown by the darker, blue or, pink, shaded regions) in which Néel or Néel-II
magnetic ordering is present on each monolayer. It is likely that this gapped
paramagnetic region comprises one or more types of VBC phase, including those
of the staggered-dimer (SDVBC) or hexagonal plaquette (PVBC) variety on each
monolayer. Another type, namely the interlayer dimer (IDVBC) variety, also
seems bound to occur, particularly for larger absolute values $|\delta|$ of
the interlayer coupling, where the dimers now occur between NN pairs of spins
across the two $AA$-stacked layers.

In practical applications of the current model, the bilayer may well rest on
top of a suitable substrate. One way to simulate the effects of such a
substrate would be to include disorder on the system. It would then be of
considerable interest to investigate which phases and, particularly, which
phase boundaries, in Fig. 8 are robust against weak on-site potential
disorder.

## ACKNOWLEDGMENTS

We thank the University of Minnesota Supercomputing Institute for the grant of
supercomputing facilities, on which some of the work reported here was
performed. The authors would also like to acknowledge the assistance given by
Research IT and the use of the HPC Pool funded by the Research Lifecycle
Programme at the University of Manchester. One of us (RFB) gratefully
acknowledges the Leverhulme Trust (United Kingdom) for the award of an
Emeritus Fellowship (Grant No. EM-2020-013).

## References

## References

  * [1] L. Savary, L. Balents, Quantum spin liquids: a review, Rep. Prog. Phys. 80 (1) (2017) 016502.  doi:10.1088/0034-4885/80/1/016502.
URL https://doi.org/10.1088%2F0034-4885%2F80%2F1%2F016502

  * [2] E. Lieb, T. Schultz, D. Mattis, Two soluble models of an antiferromagnetic chain, Ann. Phys. (N.Y.) 16 (3) (1961) 407–466.  doi:10.1016/0003-4916(61)90115-4.
URL http://www.sciencedirect.com/science/article/pii/0003491661901154

  * [3] M. B. Hastings, Lieb$-$Schultz$-$Mattis in higher dimensions, Phys. Rev. B 69 (2004) 104431.  doi:10.1103/PhysRevB.69.104431.
URL https://link.aps.org/doi/10.1103/PhysRevB.69.104431

  * [4] I. Affleck, Spin gap and symmetry breaking in CuO2 layers and other antiferromagnets, Phys. Rev. B 37 (1988) 5186–5192.  doi:10.1103/PhysRevB.37.5186.
URL https://link.aps.org/doi/10.1103/PhysRevB.37.5186

  * [5] M. Oshikawa, Commensurability, excitation gap, and topology in quantum many-particle systems on a periodic lattice, Phys. Rev. Lett. 84 (2000) 1535–1538.  doi:10.1103/PhysRevLett.84.1535.
URL https://link.aps.org/doi/10.1103/PhysRevLett.84.1535

  * [6] H. Watanabe, H. C. Po, A. Vishwanath, M. Zaletel, Filling constraints for spin-orbit coupled insulators in symmorphic and nonsymmorphic crystals, Proc. Natl. Acad. Sci. USA 112 (47) (2015) 14551–14556.
URL https://www.jstor.org/stable/26465856

  * [7] E. Rastelli, A. Tassi, L. Reatto, Non-simple magnetic order for simple Hamiltonians, Physica B & C 97 (1) (1979) 1–24.  doi:10.1016/0378-4363(79)90002-0.
URL http://www.sciencedirect.com/science/article/pii/0378436379900020

  * [8] A. Mattsson, P. Fröjdh, T. Einarsson, Frustrated honeycomb Heisenberg antiferromagnet: A Schwinger-boson approach, Phys. Rev. B 49 (1994) 3997–4002.  doi:10.1103/PhysRevB.49.3997.
URL https://link.aps.org/doi/10.1103/PhysRevB.49.3997

  * [9] J. B. Fouet, P. Sindzingre, C. Lhuillier, An investigation of the quantum ${J}_{1}$–${J}_{2}$–${J}_{3}$ model on the honeycomb lattice, Eur. Phys. J. B 20 (2) (2001) 241–254.  doi:10.1007/s100510170273.
URL https://doi.org/10.1007/s100510170273

  * [10] A. Mulder, R. Ganesh, L. Capriotti, A. Paramekanti, Spiral order by disorder and lattice nematic order in a frustrated Heisenberg antiferromagnet on the honeycomb lattice, Phys. Rev. B 81 (2010) 214419.  doi:10.1103/PhysRevB.81.214419.
URL https://link.aps.org/doi/10.1103/PhysRevB.81.214419

  * [11] S. Okumura, H. Kawamura, T. Okubo, Y. Motome, Novel spin-liquid states in the frustrated Heisenberg antiferromagnet on the honeycomb lattice, J. Phys. Soc. Jpn. 79 (11) (2010) 114705.  doi:10.1143/JPSJ.79.114705.
URL https://doi.org/10.1143/JPSJ.79.114705

  * [12] F. Wang, Schwinger boson mean field theories of spin liquid states on a honeycomb lattice: Projective symmetry group analysis and critical field theory, Phys. Rev. B 82 (2010) 024419.  doi:10.1103/PhysRevB.82.024419.
URL https://link.aps.org/doi/10.1103/PhysRevB.82.024419

  * [13] D. C. Cabra, C. A. Lamas, H. D. Rosales, Quantum disordered phase on the frustrated honeycomb lattice, Phys. Rev. B 83 (2011) 094506.  doi:10.1103/PhysRevB.83.094506.
URL https://link.aps.org/doi/10.1103/PhysRevB.83.094506

  * [14] R. Ganesh, D. N. Sheng, Y.-J. Kim, A. Paramekanti, Quantum paramagnetic ground states on the honeycomb lattice and field-induced Néel order, Phys. Rev. B 83 (2011) 144414.  doi:10.1103/PhysRevB.83.144414.
URL https://link.aps.org/doi/10.1103/PhysRevB.83.144414

  * [15] R. Ganesh, D. N. Sheng, Y.-J. Kim, A. Paramekanti, Publisher’s Note: Quantum paramagnetic ground states on the honeycomb lattice and field-induced Néel order [Phys. Rev. B 83, 144414 (2011)], Phys. Rev. B 83 (2011) 219903(E).  doi:10.1103/PhysRevB.83.219903.
URL https://link.aps.org/doi/10.1103/PhysRevB.83.219903

  * [16] B. K. Clark, D. A. Abanin, S. L. Sondhi, Nature of the spin liquid state of the Hubbard model on a honeycomb lattice, Phys. Rev. Lett. 107 (2011) 087204.  doi:10.1103/PhysRevLett.107.087204.
URL https://link.aps.org/doi/10.1103/PhysRevLett.107.087204

  * [17] D. J. J. Farnell, R. F. Bishop, P. H. Y. Li, J. Richter, C. E. Campbell, Frustrated Heisenberg antiferromagnet on the honeycomb lattice: A candidate for deconfined quantum criticality, Phys. Rev. B 84 (2011) 012403.  doi:10.1103/PhysRevB.84.012403.
URL https://link.aps.org/doi/10.1103/PhysRevB.84.012403

  * [18] J. Reuther, D. A. Abanin, R. Thomale, Magnetic order and paramagnetic phases in the quantum ${J}_{1}$-${J}_{2}$-${J}_{3}$ honeycomb model, Phys. Rev. B 84 (2011) 014417.  doi:10.1103/PhysRevB.84.014417.
URL https://link.aps.org/doi/10.1103/PhysRevB.84.014417

  * [19] A. F. Albuquerque, D. Schwandt, B. Hetényi, S. Capponi, M. Mambrini, A. M. Läuchli, Phase diagram of a frustrated quantum antiferromagnet on the honeycomb lattice: Magnetic order versus valence-bond crystal formation, Phys. Rev. B 84 (2011) 024406.  doi:10.1103/PhysRevB.84.024406.
URL https://link.aps.org/doi/10.1103/PhysRevB.84.024406

  * [20] H. Mosadeq, F. Shahbazi, S. A. Jafari, Plaquette valence bond ordering in a ${J}_{1}$–${J}_{2}$ Heisenberg antiferromagnet on a honeycomb lattice, J. Phys.: Condens. Matter 23 (22) (2011) 226006.  doi:10.1088/0953-8984/23/22/226006.
URL http://stacks.iop.org/0953-8984/23/i=22/a=226006

  * [21] J. Oitmaa, R. R. P. Singh, Phase diagram of the ${J}_{1}$-${J}_{2}$-${J}_{3}$ Heisenberg model on the honeycomb lattice: A series expansion study, Phys. Rev. B 84 (2011) 094424.  doi:10.1103/PhysRevB.84.094424.
URL https://link.aps.org/doi/10.1103/PhysRevB.84.094424

  * [22] F. Mezzacapo, M. Boninsegni, Ground-state phase diagram of the quantum ${J}_{1}$–${J}_{2}$ model on the honeycomb lattice, Phys. Rev. B 85 (2012) 060402(R).  doi:10.1103/PhysRevB.85.060402.
URL https://link.aps.org/doi/10.1103/PhysRevB.85.060402

  * [23] R. F. Bishop, P. H. Y. Li, D. J. J. Farnell, C. E. Campbell, The frustrated Heisenberg antiferromagnet on the honeycomb lattice: ${J}_{1}$–${J}_{2}$ model, J. Phys.: Condens. Matter 24 (23) (2012) 236002.  doi:10.1088/0953-8984/24/23/236002.
URL http://stacks.iop.org/0953-8984/24/i=23/a=236002

  * [24] P. H. Y. Li, R. F. Bishop, D. J. J. Farnell, C. E. Campbell, Phase diagram of a frustrated Heisenberg antiferromagnet on the honeycomb lattice: The ${J}_{1}$–${J}_{2}$–${J}_{3}$ model, Phys. Rev. B 86 (2012) 144404.  doi:10.1103/PhysRevB.86.144404.
URL https://link.aps.org/doi/10.1103/PhysRevB.86.144404

  * [25] R. F. Bishop, P. H. Y. Li, C. E. Campbell, Valence-bond crystalline order in the $s=1/2$ ${J}_{1}$–${J}_{2}$ model on the honeycomb lattice, J. Phys.: Condens. Matter 25 (30) (2013) 306002.  doi:10.1088/0953-8984/25/30/306002.
URL http://stacks.iop.org/0953-8984/25/i=30/a=306002

  * [26] H. Zhang, C. A. Lamas, Exotic disordered phases in the quantum ${J}_{1}$-${J}_{2}$ model on the honeycomb lattice, Phys. Rev. B 87 (2013) 024415.  doi:10.1103/PhysRevB.87.024415.
URL https://link.aps.org/doi/10.1103/PhysRevB.87.024415

  * [27] R. Ganesh, S. Nishimoto, J. van den Brink, Plaquette resonating valence bond state in a frustrated honeycomb antiferromagnet, Phys. Rev. B 87 (2013) 054413.  doi:10.1103/PhysRevB.87.054413.
URL https://link.aps.org/doi/10.1103/PhysRevB.87.054413

  * [28] R. Ganesh, J. van den Brink, S. Nishimoto, Deconfined criticality in the frustrated Heisenberg honeycomb antiferromagnet, Phys. Rev. Lett. 110 (2013) 127203.  doi:10.1103/PhysRevLett.110.127203.
URL https://link.aps.org/doi/10.1103/PhysRevLett.110.127203

  * [29] Z. Zhu, D. A. Huse, S. R. White, Weak plaquette valence bond order in the ${S}=1/2$ honeycomb ${J}_{1}$–${J}_{2}$ Heisenberg model, Phys. Rev. Lett. 110 (2013) 127205.  doi:10.1103/PhysRevLett.110.127205.
URL https://link.aps.org/doi/10.1103/PhysRevLett.110.127205

  * [30] S.-S. Gong, D. N. Sheng, O. I. Motrunich, M. P. A. Fisher, Phase diagram of the spin-$\frac{1}{2}$ ${J}_{1}$-${J}_{2}$ Heisenberg model on a honeycomb lattice, Phys. Rev. B 88 (2013) 165138.  doi:10.1103/PhysRevB.88.165138.
URL https://link.aps.org/doi/10.1103/PhysRevB.88.165138

  * [31] X.-L. Yu, D.-Y. Liu, P. Li, L.-J. Zou, Ground-state and finite-temperature properties of spin liquid phase in the ${J}_{1}$–${J}_{2}$ honeycomb model,, Physica E 59 (2014) 41–49.  doi:10.1016/j.physe.2013.12.017.
URL http://www.sciencedirect.com/science/article/pii/S1386947713004542

  * [32] A. Di Ciolo, J. Carrasquilla, F. Becca, M. Rigol, V. Galitski, Spiral antiferromagnets beyond the spin-wave approximation: Frustrated ${XY}$ and Heisenberg models on the honeycomb lattice, Phys. Rev. B 89 (2014) 094413.  doi:10.1103/PhysRevB.89.094413.
URL https://link.aps.org/doi/10.1103/PhysRevB.89.094413

  * [33] C.-M. Jian, M. Zaletel, Existence of featureless paramagnets on the square and the honeycomb lattices in 2+1 dimensions, Phys. Rev. B 93 (2016) 035114.  doi:10.1103/PhysRevB.93.035114.
URL https://link.aps.org/doi/10.1103/PhysRevB.93.035114

  * [34] F. Ferrari, S. Bieri, F. Becca, Competition between spin liquids and valence-bond order in the frustrated spin-$\frac{1}{2}$ Heisenberg model on the honeycomb lattice, Phys. Rev. B 96 (2017) 104401.  doi:10.1103/PhysRevB.96.104401.
URL https://link.aps.org/doi/10.1103/PhysRevB.96.104401

  * [35] A. W. Sandvik, Finite-size scaling and boundary effects in two-dimensional valence-bond solids, Phys. Rev. B 85 (2012) 134407.  doi:10.1103/PhysRevB.85.134407.
URL https://link.aps.org/doi/10.1103/PhysRevB.85.134407

  * [36] F. Coester, Bound states of a many-particle system, Nucl. Phys. 7 (1958) 421–424.  doi:10.1016/0029-5582(58)90280-3.
URL http://www.sciencedirect.com/science/article/pii/0029558258902803

  * [37] F. Coester, H. Kümmel, Short-range correlations in nuclear wave functions, Nucl. Phys. 17 (1960) 477–485.  doi:10.1016/0029-5582(60)90140-1.
URL http://www.sciencedirect.com/science/article/pii/0029558260901401

  * [38] J. Čižek, On the correlation problem in atomic and molecular systems. Calculation of wavefunction components in Ursell-type expansion using quantum-field theoretical methods, J. Chem. Phys. 45 (11) (1966) 4256–4266.  doi:10.1063/1.1727484.
URL https://doi.org/10.1063/1.1727484

  * [39] H. Kümmel, K. H. Lührmann, J. G. Zabolitzky, Many-fermion theory in exp${S}$ (or coupled cluster) form, Phys Rep. 36 (1) (1978) 1–63.  doi:10.1016/0370-1573(78)90081-9.
URL http://www.sciencedirect.com/science/article/pii/0370157378900819

  * [40] R. F. Bishop, K. H. Lührmann, Electron correlations: I. Ground-state results in the high-density regime, Phys. Rev. B 17 (1978) 3757–3780.  doi:10.1103/PhysRevB.17.3757.
URL https://link.aps.org/doi/10.1103/PhysRevB.17.3757

  * [41] R. F. Bishop, K. H. Lührmann, Electron correlations. II. Ground-state results at low and metallic densities, Phys. Rev. B 26 (1982) 5523–5557.  doi:10.1103/PhysRevB.26.5523.
URL https://link.aps.org/doi/10.1103/PhysRevB.26.5523

  * [42] J. Arponen, Variational principles and linked-cluster exp ${S}$ expansions for static and dynamic many-body problems, Ann. Phys. (N.Y.) 151 (2) (1983) 311–382.  doi:10.1016/0003-4916(83)90284-1.
URL http://www.sciencedirect.com/science/article/pii/0003491683902841

  * [43] R. F. Bishop, H. G. Kümmel, The coupled-cluster method, Phys. Today 40(3) (1987) 52–60.  doi:10.1063/1.881103.
  * [44] J. S. Arponen, R. F. Bishop, E. Pajanne, Extended coupled-cluster method. I. Generalized coherent bosonization as a mapping of quantum theory into classical Hamiltonian mechanics, Phys. Rev. A 36 (1987) 2519–2538.  doi:10.1103/PhysRevA.36.2519.
URL https://link.aps.org/doi/10.1103/PhysRevA.36.2519

  * [45] J. S. Arponen, R. F. Bishop, E. Pajanne, Extended coupled-cluster method. II. Excited states and generalized random-phase approximation, Phys. Rev. A 36 (1987) 2539–2549.  doi:10.1103/PhysRevA.36.2539.
URL https://link.aps.org/doi/10.1103/PhysRevA.36.2539

  * [46] R. J. Bartlett, Coupled-cluster approach to molecular structure and spectra: A step toward predictive quantum chemistry, J. Phys. Chem. 93 (5) (1989) 1697–1708.  doi:10.1021/j100342a008.
URL https://doi.org/10.1021/j100342a008

  * [47] J. S. Arponen, R. F. Bishop, Independent-cluster parametrizations of wave functions in model field theories. I. Introduction to their holomorphic representations, Ann. Phys. (N.Y.) 207 (1) (1991) 171–217.  doi:10.1016/0003-4916(91)90183-9.
URL http://www.sciencedirect.com/science/article/pii/0003491691901839

  * [48] R. F. Bishop, An overview of coupled cluster theory and its applications in physics, Theor. Chim. Acta 80 (1991) 95–148.  doi:10.1007/BF01119617.
  * [49] R. F. Bishop, The coupled cluster method, in: J. Navarro, A. Polls (Eds.), Microscopic Quantum Many-Body Theories and Their Applications, Lecture Notes in Physics Vol. 510, Springer-Verlag, Berlin, 1998, pp. 1–70.  doi:10.1007/BFb0104523.
  * [50] C. Zeng, D. J. J. Farnell, R. F. Bishop, An efficient implementation of high-order coupled-cluster techniques applied to quantum magnets, J. Stat. Phys. 90 (1) (1998) 327–361.  doi:10.1023/A:1023220222019.
URL https://doi.org/10.1023/A:1023220222019

  * [51] D. J. J. Farnell, R. F. Bishop, The coupled cluster method applied to quantum magnetism, in: U. Schollwöck, J. Richter, D. J. J. Farnell, R. F. Bishop (Eds.), Quantum Magnetism, Lecture Notes in Physics Vol. 645, Springer-Verlag, Berlin, 2004, pp. 307–348.  doi:10.1007/BFb0119597.
  * [52] R. J. Bartlett, M. Musiał, Coupled-cluster theory in quantum chemistry, Rev. Mod. Phys. 79 (2007) 291–352.  doi:10.1103/RevModPhys.79.291.
URL https://link.aps.org/doi/10.1103/RevModPhys.79.291

  * [53] R. F. Bishop, P. H. Y. Li, C. E. Campbell, Highly frustrated spin-lattice models of magnetism and their quantum phase transitions: A microscopic treatment via the coupled cluster method, AIP Conf. Proc. 1619 (1) (2014) 40–50.  doi:10.1063/1.4899216.
URL https://aip.scitation.org/doi/abs/10.1063/1.4899216

  * [54] D. J. J. Farnell, R. F. Bishop, J. Richter, Non-coplanar model states in quantum magnetism applications of the high-order coupled cluster method, J. Stat. Phys. 176 (1) (2019) 180–213.  doi:10.1007/s10955-019-02297-1.
URL https://doi.org/10.1007/s10955-019-02297-1

  * [55] R. F. Bishop, P. H. Y. Li, C. E. Campbell, Frustrated spin-$\frac{1}{2}$ ${J}_{1}$–${J}_{2}$ isotropic ${XY}$ model on the honeycomb lattice, Phys. Rev. B 89 (2014) 214413.  doi:10.1103/PhysRevB.89.214413.
URL https://link.aps.org/doi/10.1103/PhysRevB.89.214413

  * [56] P. H. Y. Li, R. F. Bishop, C. E. Campbell, Phase diagram of a frustrated spin-$\frac{1}{2}$ ${J}_{1}$–${J}_{2}$ ${XXZ}$ model on the honeycomb lattice, Phys. Rev. B 89 (2014) 220408(R).  doi:10.1103/PhysRevB.89.220408.
URL https://link.aps.org/doi/10.1103/PhysRevB.89.220408

  * [57] R. F. Bishop, P. H. Y. Li, O. Götze, J. Richter, C. E. Campbell, Frustrated Heisenberg antiferromagnet on the honeycomb lattice: Spin gap and low-energy parameters, Phys. Rev. B 92 (2015) 224434.  doi:10.1103/PhysRevB.92.224434.
URL https://link.aps.org/doi/10.1103/PhysRevB.92.224434

  * [58] R. F. Bishop, P. H. Y. Li, Large-$s$ expansions for the low-energy parameters of the honeycomb-lattice Heisenberg antiferromagnet with spin quantum number $s$, J. Magn. Magn. Mater. 407 (2016) 348–357.  doi:10.1016/j.jmmm.2016.01.101.
URL http://www.sciencedirect.com/science/article/pii/S0304885316301019

  * [59] P. H. Y. Li, R. F. Bishop, Ground-state phases of the spin-1 ${J}_{1}$–${J}_{2}$ Heisenberg antiferromagnet on the honeycomb lattice, Phys. Rev. B 93 (2016) 214438.  doi:10.1103/PhysRevB.93.214438.
URL https://link.aps.org/doi/10.1103/PhysRevB.93.214438

  * [60] R. F. Bishop, P. H. Y. Li, O. Götze, J. Richter, Frustrated spin-$\frac{1}{2}$ Heisenberg magnet on a square-lattice bilayer: High-order study of the quantum critical behavior of the ${J}_{1}\text{$-$}{J}_{2}\text{$-$}{J}_{1}^{\perp}$ model, Phys. Rev. B 100 (2019) 024401.  doi:10.1103/PhysRevB.100.024401.
URL https://link.aps.org/doi/10.1103/PhysRevB.100.024401

  * [61] H. Zhang, M. Arlego, C. A. Lamas, Quantum phases in the frustrated Heisenberg model on the bilayer honeycomb lattice, Phys. Rev. B 89 (2014) 024403.  doi:10.1103/PhysRevB.89.024403.
URL https://link.aps.org/doi/10.1103/PhysRevB.89.024403

  * [62] M. Arlego, C. A. Lamas, H. Zhang, Self consistent study of the quantum phases in a frustrated antiferromagnet on the bilayer honeycomb lattice, J. Phys.: Conf. Ser. 568 (4) (2014) 042019.  doi:10.1088/1742-6596/568/4/042019.
URL http://stacks.iop.org/1742-6596/568/i=4/a=042019

  * [63] R. F. Bishop, P. H. Y. Li, Frustrated honeycomb-bilayer Heisenberg antiferromagnet: The spin-$\frac{1}{2}$ ${J}_{1}$–${J}_{2}$–${J}_{1}^{\perp}$ model, Phys. Rev. B 95 (2017) 134414\.  doi:10.1103/PhysRevB.95.134414.
URL https://link.aps.org/doi/10.1103/PhysRevB.95.134414

  * [64] P. H. Y. Li, R. F. Bishop, Transverse magnetic susceptibility of a frustrated spin-$\frac{1}{2}$ ${J}_{1}$–${J}_{2}$–${J}_{1}^{\perp}$ Heisenberg antiferromagnet on a bilayer honeycomb lattice, AIP Conf. Proc. 1912 (1) (2017) 020011.  doi:10.1063/1.5016136.
URL https://aip.scitation.org/doi/abs/10.1063/1.5016136

  * [65] J. Villain, A magnetic analogue of stereoisomerism: application to helimagnetism in two dimensions, J. Phys. (France) 38 (1977) 385–391.  doi:10.1051/jphys:01977003804038500.
URL https://doi.org/10.1051/jphys:01977003804038500

  * [66] J. Villain, R. Bidaux, J.-P. Carton, R. Conte, Order as an effect of disorder, J. Phys. (France) 41 (11) (1980) 1263–1272.  doi:10.1051/jphys:0198000410110126300.
URL https://doi.org/10.1051/jphys:0198000410110126300

  * [67] P. H. Y. Li, R. F. Bishop, Collinear antiferromagnetic phases of a frustrated spin-$\frac{1}{2}$ ${J}_{1}$–${J}_{2}$–${J}_{1}^{\perp}$ Heisenberg model on an ${AA}$-stacked bilayer honeycomb lattice, J. Magn. Magn. Mater. 482 (2019) 262–273.  doi:10.1016/j.jmmm.2019.03.033.
URL http://www.sciencedirect.com/science/article/pii/S0304885318336667

  * [68] S.-S. Gong, W. Zhu, D. N. Sheng, Quantum phase diagram of the spin-1 ${J}_{1}$-${J}_{2}$ Heisenberg model on the honeycomb lattice, Phys. Rev. B 92 (2015) 195110.  doi:10.1103/PhysRevB.92.195110.
URL https://link.aps.org/doi/10.1103/PhysRevB.92.195110

  * [69] We use the program package CCCM of D. J. J. Farnell and J. Schulenburg, see http://www-e.uni-magdeburg.de/jschulen/ccm/index.html.  [link].
URL http://www-e.uni-magdeburg.de/jschulen/ccm/index.html

  * [70] R. F. Bishop, D. J. J. Farnell, S. E. Krüger, J. B. Parkinson, J. Richter, C. Zeng, High-order coupled cluster method calculations for the ground- and excited-state properties of the spin-half ${XXZ}$ model, J. Phys.: Condens. Matter 12 (30) (2000) 6887–6902.  doi:10.1088/0953-8984/12/30/317.
URL http://stacks.iop.org/0953-8984/12/i=30/a=317

  * [71] S. E. Krüger, J. Richter, J. Schulenburg, D. J. J. Farnell, R. F. Bishop, Quantum phase transitions of a square-lattice Heisenberg antiferromagnet with two kinds of nearest-neighbor bonds: A high-order coupled-cluster treatment, Phys. Rev. B 61 (2000) 14607–14615.  doi:10.1103/PhysRevB.61.14607.
URL https://link.aps.org/doi/10.1103/PhysRevB.61.14607

  * [72] D. J. J. Farnell, K. A. Gernoth, R. F. Bishop, High-order coupled-cluster method for general spin-lattice problems: An illustration via the anisotropic Heisenberg model, Phys. Rev. B 64 (2001) 172409.  doi:10.1103/PhysRevB.64.172409.
URL https://link.aps.org/doi/10.1103/PhysRevB.64.172409

  * [73] R. Darradi, J. Richter, D. J. J. Farnell, Coupled cluster treatment of the Shastry-Sutherland antiferromagnet, Phys. Rev. B 72 (2005) 104425.  doi:10.1103/PhysRevB.72.104425.
URL https://link.aps.org/doi/10.1103/PhysRevB.72.104425

  * [74] R. F. Bishop, P. H. Y. Li, R. Darradi, J. Richter, The quantum ${J}_{1}$–${J}_{1}^{\prime}$–${J}_{2}$ spin-$1$ Heisenberg model: Influence of the interchain coupling on the ground-state magnetic ordering in 2D, EPL 83 (4) (2008) 47004.  doi:10.1209/0295-5075/83/47004.
URL http://stacks.iop.org/0295-5075/83/i=4/a=47004

  * [75] R. F. Bishop, P. H. Y. Li, R. Darradi, J. Richter, C. E. Campbell, The effect of anisotropy on the ground-state magnetic ordering of the spin-1 quantum ${J}_{1}^{XXZ}$–${J}_{2}^{XXZ}$ model on the square lattice, J. Phys.: Condens. Matter 20 (41) (2008) 415213.  doi:10.1088/0953-8984/20/41/415213.
URL http://stacks.iop.org/0953-8984/20/i=41/a=415213

  * [76] R. F. Bishop, P. H. Y. Li, D. J. J. Farnell, C. E. Campbell, Magnetic order in a spin-$\frac{1}{2}$ interpolating square-triangle Heisenberg antiferromagnet, Phys. Rev. B 79 (2009) 174405.  doi:10.1103/PhysRevB.79.174405.
URL https://link.aps.org/doi/10.1103/PhysRevB.79.174405

  * [77] R. F. Bishop, P. H. Y. Li, D. J. J. Farnell, C. E. Campbell, Magnetic order on a frustrated spin-$\frac{1}{2}$ Heisenberg antiferromagnet on the Union Jack lattice, Phys. Rev. B 82 (2010) 024416.  doi:10.1103/PhysRevB.82.024416.
URL https://link.aps.org/doi/10.1103/PhysRevB.82.024416

  * [78] R. F. Bishop, P. H. Y. Li, D. J. J. Farnell, C. E. Campbell, Magnetic order in a spin-$\frac{1}{2}$ interpolating kagome/square Heisenberg antiferromagnet, Phys. Rev. B 82 (2010) 104406.  doi:10.1103/PhysRevB.82.104406.
URL https://link.aps.org/doi/10.1103/PhysRevB.82.104406

  * [79] R. F. Bishop, P. H. Y. Li, A frustrated quantum spin-${\boldmath s}$ model on the Union Jack lattice with spins ${\boldmath s>\frac{1}{2}}$, Eur. Phys. J. B 81 (1) (2011) 37–48.  doi:10.1140/epjb/e2011-10817-7.
URL https://doi.org/10.1140/epjb/e2011-10817-7

  * [80] P. H. Y. Li, R. F. Bishop, D. J. J. Farnell, J. Richter, C. E. Campbell, Ground-state phases of the frustrated spin-$\frac{1}{2}$ ${J}_{1}$–${J}_{2}$–${J}_{3}$ Heisenberg ferromagnet (${J}_{1}<0$) on the honeycomb lattice with ${J}_{3}={J}_{2}>0$, Phys. Rev. B 85 (2012) 085115.  doi:10.1103/PhysRevB.85.085115.
URL https://link.aps.org/doi/10.1103/PhysRevB.85.085115

  * [81] R. F. Bishop, P. H. Y. Li, Complete phase diagram of the spin-$\frac{1}{2}$ ${J}_{1}$–${J}_{2}$–${J}_{3}$ model (with ${J}_{3}={J}_{2}$) on the honeycomb lattice, Phys. Rev. B 85 (2012) 155135.  doi:10.1103/PhysRevB.85.155135.
URL https://link.aps.org/doi/10.1103/PhysRevB.85.155135

  * [82] P. H. Y. Li, R. F. Bishop, Magnetic order in spin-1 and spin-$\frac{3}{2}$ interpolating square-triangle Heisenberg antiferromagnets, Eur. Phys. J. B 85 (2012) 25.  doi:10.1140/epjb/e2011-20766-8.
URL https://doi-org.manchester.idm.oclc.org/10.1140/epjb/e2011-20766-8

  * [83] P. H. Y. Li, R. F. Bishop, C. E. Campbell, D. J. J. Farnell, O. Götze, J. Richter, Spin-$\frac{1}{2}$ Heisenberg antiferromagnet on an anisotropic kagome lattice, Phys. Rev. B 86 (2012) 214403.  doi:10.1103/PhysRevB.86.214403.
URL https://link.aps.org/doi/10.1103/PhysRevB.86.214403

  * [84] D. J. J. Farnell, O. Götze, J. Richter, R. F. Bishop, P. H. Y. Li, Quantum $s=\frac{1}{2}$ antiferromagnets on Archimedean lattices: The route from semiclassical magnetic order to nonmagnetic quantum states, Phys. Rev. B 89 (2014) 184407.  doi:10.1103/PhysRevB.89.184407.
URL https://link.aps.org/doi/10.1103/PhysRevB.89.184407

  * [85] P. H. Y. Li, R. F. Bishop, C. E. Campbell, Frustrated Heisenberg antiferromagnet on the honeycomb lattice with spin quantum number $s\geq 1$, J. Phys.: Conf. Ser. 702 (1) (2016) 012001.  doi:10.1088/1742-6596/702/1/012001.
URL http://stacks.iop.org/1742-6596/702/i=1/a=012001

  * [86] R. F. Bishop, P. H. Y. Li, R. Darradi, J. Richter, The quantum ${J}_{1}$–${J}_{1}^{\prime}$–${J}_{2}$ spin-$1/2$ Heisenberg model: influence of the interchain coupling on the ground-state magnetic ordering in two dimensions, J. Phys.: Condens. Matter 20 (25) (2008) 255251.  doi:10.1088/0953-8984/20/25/255251.
URL http://stacks.iop.org/0953-8984/20/i=25/a=255251

  * [87] R. F. Bishop, P. H. Y. Li, R. Darradi, J. Schulenburg, J. Richter, Effect of anisotropy on the ground-state magnetic ordering of the spin-half quantum ${J}_{1}^{XXZ}$–${J}_{2}^{XXZ}$ model on the square lattice, Phys. Rev. B 78 (2008) 054412.  doi:10.1103/PhysRevB.78.054412.
URL https://link.aps.org/doi/10.1103/PhysRevB.78.054412

  * [88] R. Darradi, O. Derzhko, R. Zinke, J. Schulenburg, S. E. Krüger, J. Richter, Ground state phases of the spin-1/2 ${J}_{1}$–${J}_{2}$ Heisenberg antiferromagnet on the square lattice: A high-order coupled cluster treatment, Phys. Rev. B 78 (2008) 214415.  doi:10.1103/PhysRevB.78.214415.
URL https://link.aps.org/doi/10.1103/PhysRevB.78.214415

  * [89] J. Reuther, P. Wölfle, R. Darradi, W. Brenig, M. Arlego, J. Richter, Quantum phases of the planar antiferromagnetic ${J}_{1}$–${J}_{2}$–${J}_{3}$ Heisenberg model, Phys. Rev. B 83 (2011) 064416.  doi:10.1103/PhysRevB.83.064416.
URL https://link.aps.org/doi/10.1103/PhysRevB.83.064416

  * [90] O. Götze, S. E. Krüger, F. Fleck, J. Schulenburg, J. Richter, Ground-state phase diagram of the spin-$\frac{1}{2}$ square-lattice ${J}_{1}$-${J}_{2}$ model with plaquette structure, Phys. Rev. B 85 (2012) 224424.  doi:10.1103/PhysRevB.85.224424.
URL https://link.aps.org/doi/10.1103/PhysRevB.85.224424

  * [91] R. F. Bishop, P. H. Y. Li, High-order study of the quantum critical behavior of a frustrated spin-$\frac{1}{2}$ antiferromagnet on a stacked honeycomb bilayer, Phys. Rev. B 96 (2017) 224416\.  doi:10.1103/PhysRevB.96.224416.
URL https://link.aps.org/doi/10.1103/PhysRevB.96.224416