Frustrated spin-$\frac{1}{2}$ Heisenberg magnet on a square-lattice bilayer: High-order study of the quantum critical behavior of the $J_{1}$--$J_{2}$--$J_{1}^{\perp}$ model

TL;DR

This study uses high-order coupled cluster calculations to map the quantum phase diagram of a frustrated spin-1/2 bilayer Heisenberg model, identifying phase boundaries and a quantum triple point with high precision.

Contribution

The paper provides the first high-order coupled cluster analysis of the $J_1$--$J_2$--$J_1^ot$ model, accurately locating phase boundaries and the quantum triple point.

Findings

Identified the quantum triple point at approximately ($rac{0.547}{rac{0.45}$) in the $rac{ ext{J}_2}{ ext{J}_1}$--$rac{ ext{J}_1^ot}{ ext{J}_1}$ plane.

Mapped the phase boundaries of collinear antiferromagnetic phases with high precision.

Showed the absence of a quantum triple point in the $rac{ ext{J}_2}{ ext{J}_1}$--$rac{ ext{J}_1^ot}{ ext{J}_1}$ plane for $rac{ ext{J}_1^ot}{ ext{J}_1} > 0$.

Abstract

The zero-temperature phase diagram of the spin-$\frac{1}{2}$ $J_{1}$--$J_{2}$--$J_{1}^{\perp}$ model on an $AA$-stacked square-lattice bilayer is studied using the coupled cluster method implemented to very high orders. Both nearest-neighbor (NN) and frustrating next-nearest-neighbor Heisenberg exchange interactions, of strengths $J_{1}>0$ and $J_{2} \equiv \kappa J_{1}>0$, respectively, are included in each layer. The two layers are coupled via a NN interlayer Heisenberg exchange interaction with a strength $J_{1}^{\perp} \equiv \delta J_{1}$. The magnetic order parameter $M$ (viz., the sublattice magnetization) is calculated directly in the thermodynamic (infinite-lattice) limit for the two cases when both layers have antiferromagnetic ordering of either the N\'{e}el or the striped kind, and with the layers coupled so that NN spins between them are either parallel (when $\delta < 0$) or antiparallel (when $\delta > 0$) to one another. Calculations are performed at $n$th order in a well-defined sequence of approximations, which exactly preserve both the Goldstone linked cluster theorem and the Hellmann-Feynman theorem, with $n \leq 10$. The sole approximation made is to extrapolate such sequences of $n$th-order results for $M$ to the exact limit, $n \to \infty$. By thus locating the points where $M$ vanishes, we calculate the full phase boundaries of the two collinear AFM phases in the $\kappa$--$\delta$ half-plane with $\kappa > 0$. In particular, we provide the accurate estimate, ($\kappa \approx 0.547,\delta \approx -0.45$), for the position of the quantum triple point (QTP) in the region $\delta < 0$. We also show that there is no counterpart of such a QTP in the region $\delta > 0$, where the two quasiclassical phase boundaries show instead an ``avoided crossing'' behavior, such that the entire region that contains the nonclassical paramagnetic phases is singly connected.

Equations

$$H$$

$$+J_1^{\perp }\sum _i{\mathbf{s}}_{i,1}\cdot {\mathbf{s}}_{i,2}$$

$$\equiv J_{1}h(\kappa ,\delta );\kappa \equiv J_2/J_1,\delta \equiv J_1^{\perp }/J_1,$$

$$\sum _I{C}_I^{+}|\Phi ⟩⟨\Phi |C_I^{-}=1,$$

$$⟨\Phi |C_I^{+}=0=C_I^{-}|\Phi ⟩,\forall I\ne 0;C_0^{+}\equiv 1,$$

$$[C_I^{+},C_J^{+}]=0=[C_I^{-},C_J^{-}],$$

$$⟨\Phi |C_I^{-}{C}_J^{+}|\Phi ⟩={\delta }_{I,J},$$

$$H|\Psi ⟩=E|\Psi ⟩;⟨\tilde{\Psi }|H=E⟨\tilde{\Psi }|,$$

$$|\Psi ⟩={\mathrm{e}}^S|\Phi ⟩;S=\sum _{I\ne 0}{\mathcal{S}}_I{C}_I^{+},$$

$$⟨\tilde{\Psi }|=⟨\Phi |\tilde{S}{\mathrm{e}}^{-S};\tilde{S}=1+\sum _{I\ne 0}{\tilde{\mathcal{S}}}_I{C}_I^{-},$$

$$⟨\Phi |\tilde{S}=\frac{⟨\Phi |{\mathrm{e}}^{S^{†}}{\mathrm{e}}^S}{⟨\Phi |{\mathrm{e}}^{S^{†}}{\mathrm{e}}^S|\Phi ⟩},$$

$$\bar{H}\equiv ⟨\tilde{\Psi }|H|\Psi ⟩=⟨\Phi |\tilde{S}{\mathrm{e}}^{-S}H{\mathrm{e}}^S|\Phi ⟩,$$

$$⟨\Phi |C_I^{-}{\mathrm{e}}^{-S}H{\mathrm{e}}^S|\Phi ⟩=0,\forall I\ne 0,$$

$$⟨\Phi |\tilde{S}({\mathrm{e}}^{-S}[H,C_I^{+}]{\mathrm{e}}^S|\Phi ⟩=0,\forall I\ne 0,$$

$$E=E({\mathcal{S}}_I)=⟨\Phi |{\mathrm{e}}^{-S}H{\mathrm{e}}^S|\Phi ⟩.$$

$$⟨\Phi |\tilde{S}({\mathrm{e}}^{-S}H{\mathrm{e}}^S-E)C_I^{+}|\Phi ⟩=0,\forall I\ne 0.$$

$$\bar{A}=\bar{A}({\mathcal{S}}_I,{\tilde{\mathcal{S}}}_I)=⟨\Phi |\tilde{S}{\mathrm{e}}^{-S}A{\mathrm{e}}^S|\Phi ⟩.$$

$${\mathrm{e}}^{-S}H{\mathrm{e}}^S=\sum _{n=0}^{\infty }\frac{1}{n!}[H,S{]}_n,$$

$$[H,S{]}_n\equiv [[H,S{]}_{n-1},S];[H,S{]}_0=H.$$

$$M=-\frac{1}{N}\sum _{i=1}^N⟨\Phi |\tilde{S}{\mathrm{e}}^{-S}{s}_{l_i}^z{\mathrm{e}}^S|\Phi ⟩.$$

$$M(n)={\mu }_0+{\mu }_1{n}^{-1/2}+{\mu }_2{n}^{-3/2}.$$

$$e(n)=e_0+e_1{n}^{-2}+e_2{n}^{-4}.$$

Full text

I INTRODUCTION

The frustrated spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$ Heisenberg antiferromagnet on the square lattice, which contains isotropic Heisenberg exchange interactions with strengths $J_{1}>0$ between all nearest-neighbor (NN) pairs of spins and $J_{2}>0$ between all next-nearest-neighbor (NNN) pairs, has become a paradigmatic model of quantum magnetism. It has received enormous attention over the last thirty or so years Chandra and Doucot (1988); Dagotto and Moreo (1989); Gelfand et al. (1989); Sachdev and Bhatt (1990); Chubukov and Jolicoeur (1991); Read and Sachdev (1991); Richter (1993); Richter et al. (1994); Ivanov and Richter (1994); Schulz et al. (1996); Oitmaa and Weihong (1996); Zhitomirsky and Ueda (1996); Trumper et al. (1997); Bishop et al. (1998); Singh et al. (1999); Kotov et al. (1999); Capriotti and Sorella (2000); Capriotti et al. (2001); Takano et al. (2003); Roscilde et al. (2004); Lante and Parola (2006); Sirker et al. (2006); Schmalfuß et al. (2006); Mambrini et al. (2006); Bishop et al. (2008a, b); Darradi et al. (2008); Isaev et al. (2009); Murg et al. (2009); Ralko et al. (2009); Richter and Schulenburg (2010); Reuther and Wölfle (2010); Reuther et al. (2011); Yu and Kao (2012); Götze et al. (2012); Jiang et al. (2012); Mezzacapo (2012); Li et al. (2012); Wang et al. (2013); Zhang and Beach (2013); Hu et al. (2013); Gong et al. (2014); Doretto (2014); Qi and Gu (2014); Metavitsiadis et al. (2014); Ren et al. (2014); Wang (2014); Chou and Chen (2014); Morita et al. (2015); Richter et al. (2015); Wang et al. (2016); Poilblanc and Mambrini (2017); Haghshenas and Sheng (2018); Yu et al. (2018); Wang and Sandvik (2018); Liu et al. (2018), starting with its proposed relationship to the disappearance of antiferromagnetic (AFM) long-range order (LRO) in the high-$T_{c}$ cuprate superconductors. The conjecture here was that frustrated AFM exchange couplings might lead to a quantum spin liquid (QSL) state in which preformed pairs, or resonating valence bonds, could become superconducting upon doping Anderson (1987); Lee et al. (2006). More recently, as frustrated quantum magnets have emerged as an active research field in their own right, the model has become recognized as one of the most challenging quantum spin-lattice systems. Accordingly, it has been widely studied Chandra and Doucot (1988); Dagotto and Moreo (1989); Gelfand et al. (1989); Sachdev and Bhatt (1990); Chubukov and Jolicoeur (1991); Read and Sachdev (1991); Richter (1993); Richter et al. (1994); Ivanov and Richter (1994); Schulz et al. (1996); Oitmaa and Weihong (1996); Zhitomirsky and Ueda (1996); Trumper et al. (1997); Bishop et al. (1998); Singh et al. (1999); Kotov et al. (1999); Capriotti and Sorella (2000); Capriotti et al. (2001); Takano et al. (2003); Roscilde et al. (2004); Lante and Parola (2006); Sirker et al. (2006); Schmalfuß et al. (2006); Mambrini et al. (2006); Bishop et al. (2008a, b); Darradi et al. (2008); Isaev et al. (2009); Murg et al. (2009); Ralko et al. (2009); Richter and Schulenburg (2010); Reuther and Wölfle (2010); Reuther et al. (2011); Yu and Kao (2012); Götze et al. (2012); Jiang et al. (2012); Mezzacapo (2012); Li et al. (2012); Wang et al. (2013); Zhang and Beach (2013); Hu et al. (2013); Gong et al. (2014); Doretto (2014); Qi and Gu (2014); Metavitsiadis et al. (2014); Ren et al. (2014); Wang (2014); Chou and Chen (2014); Morita et al. (2015); Richter et al. (2015); Wang et al. (2016); Poilblanc and Mambrini (2017); Haghshenas and Sheng (2018); Yu et al. (2018); Wang and Sandvik (2018); Liu et al. (2018), by a large number of theoretical techniques, as a prototypical system in which to examine quantum phase transition (QPTs) between quasiclassical ground-state (GS) phases with magnetic LRO and magnetically disordered (paramagnetic) quantum phases that are driven by frustration.

In addition to this extensive theoretical interest, it is also worth noting that several good experimental realizations of spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$ models on a quasi-two-dimensional square lattice exist with $J_{1}>0$ and $J_{2}>0$. Examples include the vanadium-layered oxide materials Li2VO(Si,Ge)O4 Melzi et al. (2000); Melzi:2001_sqLatt_J1J2mod_merge; Rosner:2003_sqLatt_J1J2mod_merge and the $B$-site ordered double-perovskite oxides Ba2CuWO6 Todate et al. (2007), Sr2CuMoO6 Vasala et al. (2014a), Sr2CuWO6 Vasala et al. (2014a, b), and Sr2CuTeO6 Koga et al. (2016).

Despite the intense interest in this model from both theorists and experimentalists, as outlined above, the nature of its GS phase around the value $\kappa=\frac{1}{2}$ of the frustration parameter, $\kappa\equiv J_{2}/J_{1}$, which represents the point of maximum frustration in the classical version of the model, still remains largely unresolved. Thus, if the spins on the square-lattice sites carry spin quantum number $s$, the model becomes classical in the limit $s\to\infty$. In this classical limit the GS phase is simple in the two limiting cases $\kappa=0$ and $\kappa\to\infty$. Clearly, when $\kappa=0$, the model has Néel AFM order [i.e., with a magnetic wave vector $\mathbf{Q}=(\pi,\pi)$]. For nonzero values of $\kappa$ the Néel-ordered state has a GS energy per spin given by $E_{{\rm cl}}/N=2J_{1}(-1+\kappa)s^{2}$. By contrast, when $\kappa\to\infty$, the classical ordering is such that on each of the two equivalent interpenetrating sublattices (i.e., which comprise NNN sites on the original lattice connected by $J_{2}$ bonds) the spins are separately Néel-ordered, and with a relative angle $\theta$ between the ordering directions on the two sublattices. The GS energy per spin in this case is given by $E_{{\rm cl}}/N=-2J_{1}\kappa s^{2}$, independent of $\theta$, for any value of $\kappa$.

Clearly, the classical $J_{1}$–$J_{2}$ model on the square lattice thus has a first-order phase transition at $\kappa=\kappa_{{\rm cl}}=\frac{1}{2}$ between two AFM states, viz., the Néel state for $\kappa<\frac{1}{2}$ and an infinitely degenerate family of ground states specified by the relative angle $\theta$ between the ordering directions on the two interpenetrating sublattices, for $\kappa>\frac{1}{2}$. Thus, for $\kappa>\frac{1}{2}$, the classical GS manifold has $\text{SU(2)}\times\text{SU(2)}$ symmetry, which is larger than the SU(2) symmetry of the Hamiltonian. In this latter case, although the effects of the exchange fields ($J_{1}$) between the two sublattices cancel out, the zero-point quantum fluctuations, as well as the thermal fluctuations, will depend on the angle $\theta$ between the two sublattice spin orientations. This leads to a prototypical example Chandra et al. (1990) of the phenomenon of order by disorder Villain (1977); Villain:1980_ordByDisord_merge; Shender (1982), whereby the GS degeneracy is lifted by quantum fluctuations with the angle $\theta$ now selected to be 0 or $\pi$. The AFM GS ordering is now collinear, and the corresponding GS phase is a striped one consisting of successive alternating columns (or rows) of parallel spins [i.e., with a magnetic wave vector $\mathbf{Q}=(\pi,0)$ or $\mathbf{Q}=(0,\pi)$, respectively]. The GS symmetry is thereby reduced from $\text{SU(2)}\times\text{SU(2)}$ to $\text{SU(2)}\times\text{Z}_{2}$, and the collinear striped state breaks the invariance of the Heisenberg Hamiltonian under both spin rotations [SU(2)] and rotations by 90∘* of the square lattice [$\text{Z}_{2}$].

In the classical, $s\to\infty$, limit of the model, lowest-order spin-wave theory, wherein the effects of quantum fluctuations are taken into account perturbatively at $O(s^{-1})$, thus shows Chandra and Doucot (1988) that the critical coupling $\kappa_{{\rm cl}}=\frac{1}{2}$ marks a first-order transition between the Néel and striped collinear AFM phases. In the extreme quantum case, $s=\frac{1}{2}$, in which we are interested here, where quantum fluctuations now have to be taken fully into account beyond perturbation theory, it may be anticipated that these two quasiclassical AFM phases persist, but are now separated by one or more intermediate paramagnetic phases with no classical counterparts (i.e., without magnetic LRO). While there is essentially complete consensus that this scenario is realized in the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$ model on the square lattice, the nature of both the phase (or phases) in the intermediate regime and their associated QPTs, as well as the precise critical values of $\kappa$ at which the latter occur, are still not completely resolved, despite many calculations over the last thirty or so years. These have included investigations of the model using a wide diversity of modern theoretical techniques and numerical tools of ever increasing sophistication. Examples include those based on mean-field theories of various (e.g., cluster, hierarchical) types Gelfand et al. (1989); Isaev et al. (2009); Reuther and Wölfle (2010); Ren et al. (2014), the exact diagonalization (ED) of finite-sized clusters Dagotto and Moreo (1989); Ivanov and Richter (1994); Schulz et al. (1996); Capriotti and Sorella (2000); Roscilde et al. (2004); Mambrini et al. (2006); Richter and Schulenburg (2010); Götze et al. (2012), linked-cluster series expansions Gelfand et al. (1989); Oitmaa and Weihong (1996); Singh et al. (1999); Sirker et al. (2006); Reuther et al. (2011), the bond-operator formalism Sachdev and Bhatt (1990); Zhitomirsky and Ueda (1996); Doretto (2014), resonating valence bond (RVB) approaches Capriotti et al. (2001); Li et al. (2012); Wang et al. (2013); Zhang and Beach (2013); Hu et al. (2013); Wang (2014); Chou and Chen (2014), variational Monte Carlo (VMC) approaches based on various families of trial GS wave functions (e.g., RVB states, entangled plaquette states) Capriotti et al. (2001); Mezzacapo (2012); Li et al. (2012); Wang et al. (2013); Zhang and Beach (2013); Hu et al. (2013); Wang (2014); Chou and Chen (2014); Morita et al. (2015); Yu et al. (2018), various quantum field-theoretical approaches Takano et al. (2003); Lante and Parola (2006); Ralko et al. (2009); Chandra et al. (1990) including the dynamic functional renormalization group Reuther and Wölfle (2010); Reuther et al. (2011), the density-matrix renormalization group (DMRG) Jiang et al. (2012); Gong et al. (2014); Wang and Sandvik (2018), matrix-product or tensor-network approaches Murg et al. (2009); Yu and Kao (2012); Wang et al. (2013, 2016); Poilblanc and Mambrini (2017); Haghshenas and Sheng (2018); Liu et al. (2018), and the coupled cluster method (CCM) Schmalfuß et al. (2006); Bishop et al. (2008a, b); Darradi et al. (2008); Reuther et al. (2011); Götze et al. (2012); Richter et al. (2015).

While, in the intermediate region, where the GS phase or phases are non-magnetic, the SU(2) spin symmetry is not broken, various symmetries of the lattice still may or may not be broken. In the former case one can have various valence-bond crystalline (VBC) phases where the lattice symmetries are broken by the formation of some pattern of spin singlets. Examples include the columnar dimer VBC phase, which breaks both translational and rotational lattice symmetries, and the plaquette VBC phase, which breaks only the translational symmetry. Alternatively, one could have a QSL phase that conserves all lattice symmetries. Such a QSL phase could be either gapped or gapless (e.g., of the $\text{Z}_{2}$ type).

Each of these phases has been proposed to form the stable GS in part of all of the paramagnetic intermediate regime of the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$ model on the square lattice by various of the above-cited references, with no overall consensus having yet emerged. Part of the reason for this uncertainty undoubtedly must lie in the fact that of the various methods discussed above that have high potential accuracy and/or are capable of systematic improvement via some well-defined hierarchical approximation scheme, almost all are either intrinsically biased in favor of some particular GS phase and/or are not directly performed in the thermodynamic (infinite-lattice) limit of interest. In the latter regard, for example, the great majority of the techniques employed are performed on lattices of a finite size ($N$ spins), and some form of finite-size scaling is then used to extrapolate to the thermodynamic ($N\to\infty$) limit.

As has been very rigorously and authoritatively demonstrated in a recent study Sandvik (2012) of the spin-$\frac{1}{2}$ $J$–$Q$ model on the square lattice, for which the infamous quantum Monte Carlo (QMC) minus-sign problem is absent, and hence where large-scale QMC calculations can be undertaken, by contrast with the corresponding $J_{1}$–$J_{2}$ model of interest here, such extrapolations to the thermodynamic limit can have great uncertainties. This is specially true in cases where is little or no analytic guidance from theoretical considerations, as is often the case, but can also even hold when such guidance is present. In this context it is particularly noteworthy that the CCM Coester (1958); Coester and Kümmel (1960); Čižek (1966); Kümmel et al. (1978); Bishop and Lührmann (1978, 1982); Arponen (1983); Bishop and Kümmel (1987); Arponen et al. (1987a, b); Bartlett (1989); Arponen and Bishop (1991); Bishop (1991, 1998); Zeng et al. (1998); Farnell and Bishop (2004); Bartlett and Musiał (2007); Bishop et al. (2014) provides a rather singular example of a theoretical quantum many-body technique that can and does study arbitrary spin-lattice models directly in the thermodynamic limit. It is precisely for that reason that we employ it here.

Furthermore, in view of the still puzzling nature of the phase or phases present in the intermediate paramagnetic regime of the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$ model on the square lattice around the value $\kappa=\frac{1}{2}$ of the frustration parameter, it is also potentially useful to examine a larger class of systems for which this model reduces to a special case. Thus, we are strongly motivated to consider the corresponding spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on a square-lattice bilayer. Each of the two monolayers is just a frustrated $J_{1}$–$J_{2}$ system, but the two layers are now connected by Heisenberg exchange bonds of strength $J_{1}^{\perp}\equiv\delta J_{1}$ between NN interlayer pairs of spins, with the two layer arranged in $AA$ stacking [i.e., with each site of one (horizontal) monolayer placed immediately above its counterpart on the other monolayer]. The original $J_{1}$–$J_{2}$ model is then just the special case $\delta=0$ of the larger $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model. In this paper we use the CCM to study the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on a square-lattice bilayer, for both signs of the interlayer coupling parameter $\delta$. In particular, we will concentrate our efforts on examining the complete phase boundaries of the two quasiclassical collinear (AFM) phases (viz., the phases with Néel and striped AFM order on each of the coupled monolayers) in the $\kappa$–$\delta$ half-plane with $\kappa>0$ (and $J_{1}>0$), and specifically in the window $0\leq\kappa\leq 1$ that contains the intermediate paramagnetic regime in the case $\delta=0$.

In this context it is interesting to note too that the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model has also been previously studied on a stacked square lattice (i.e., where the number of layers $n\to\infty$, rather than the case $n=2$ studied here) Schmalfuß et al. (2006), where use was also made of the CCM. Thus, the bilayer model we study here lies, in some sense, between the strictly two-dimensional square-lattice $J_{1}$–$J_{2}$ model (i.e., where $\delta=0$) and the strictly three-dimensional $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the stacked square lattice with an infinite number of layers. Furthermore, unlike the latter case, the bilayer case also exhibits the additional physical phenomenon of dimerization between NN interlayer pairs, as discussed more fully in Sec. II. These features thus provide considerable additional motivation to study the bilayer model.

The plan for the remainder of this paper is as follows. The $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model is itself first described in Sec. II, where we also discuss more fully the main features of the limiting case, $J_{1}^{\perp}=0$, of the monolayer model. We also give there some discussion of what we might expect to be some of the main features of the phase boundaries of the two quasiclassical AFM phases as the interlayer coupling parameter, $\delta\equiv J_{1}^{\perp}/J_{1}$, is introduced. The main features of the CCM as applied to quantum spin-lattice problems are then reviewed in Sec. III before our numerical results are presented in Sec. IV. Finally, our findings are summarized and discussed in Sec. V, where we also make comparisons with the results of others.

II THE MODEL

The Hamiltonian of the $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on a square-lattice bilayer is specified as

[TABLE]

such that the sites on each (horizontal) monolayer are labelled by the index $i$ (i.e., with the two layers in $AA$ stacking such that sites $i$ on the top layer lie vertically above those on the bottom layer), and the two layers are labelled by the index $\alpha=1,2$. Every site $(i,\alpha)$ is occupied by a spin-$s$ particle described in terms of the usual SU(2) operators ${\bf s}_{i,\alpha}\equiv(s^{x}_{i,\alpha},s^{y}_{i,\alpha},s^{z}_{i,\alpha})$, with ${\bf s}^{2}_{i,\alpha}=s(s+1)\mathbb{1}$, and where we restrict discussion here to the case $s=\frac{1}{2}$. The first two sums over $\langle i,j\rangle$ and $\langle\langle i,k\rangle\rangle$ in Eq. (1) run over all NN and NNN intralayer pairs of spins, respectively, with each Heisenberg bond (with respective strengths $J_{1}$ and $J_{2}$) counted once and once only. The third sum in Eq. (1) over the index $i$ counts all corresponding interlayer NN Heisenberg bonds of strength $J_{1}^{\perp}$. We shall be interested here in the case when both intralayer bonds are AFM in nature (i.e., $J_{1}>0$ and $J_{2}\equiv\kappa J_{1}>0$), such that frustration is present in each monolayer, but where the interlayer coupling parameter, $J_{1}^{\perp}\equiv\delta J_{1}$, may be either AFM ($\delta>0$) or ferromagnetic (FM) ($\delta<0$) in nature. Since the parameter $J_{1}$ merely sets the overall energy scale, the Hamiltonian may be expressed as in the last line of Eq. (1), such that the relevant parameters of the model are $\kappa$ and $\delta$.

Our main interest here will thus be to investigate the regions of stability of the two collinear AFM phases in each monolayer (i.e., the quasiclassical Néel and striped phases) in the $\kappa$–$\delta$ half-plane with $\kappa>0$, as the interlayer coupling, $\delta$, is turned on. The square-lattice bilayer is illustrated in Fig. 1(a), while the patterns of spins of the two quasiclassical AFM phases on each monolayer are shown in Fig. 1(b) and 1(c), respectively.

For both AFM phases the original square lattice (with spacing $d$) is decomposed into two equivalent sublattices. For the Néel state each sublattice is itself square (i.e., with spacing $\sqrt{2}d\times\sqrt{2}d$), such that each site on one sublattice has its 4 NN sites on the original lattice on the other sublattice. By contrast, for the striped states, the original lattice is decomposed into equivalent sublattices chosen either as alternating columns (each with spacing $2d\times d$) or as alternating rows (each with spacing $d\times 2d$). Each of the corresponding classical AFM phases then has its spins on one sublattice pointing in a given, arbitrary (say, down) direction, and those on the other sublattice pointing in the opposite (say, up) direction. The classical Néel state is thus as shown in Fig. 1(b), while the classical columnar striped state is as shown in Fig. 1(c).

For the case of the spin-$\frac{1}{2}$ square-lattice monolayer ($\delta=0$) there is essentially complete agreement that Néel order persists for $\kappa<\kappa_{c_{1}}$, striped order persists for $\kappa>\kappa_{c_{2}}$, and some paramagnetic phase (or phases) form the stable GS phase in the intermediate regime $\kappa_{c_{1}}<\kappa<\kappa_{c_{2}}$. While the nature of the phase(s) in the intermediate regime remains unresolved, as noted in Sec. I, modern high-quality calculations do seem to be converging on values for the two critical points of $\kappa_{c_{1}}\approx 0.43(3)$ and $\kappa_{c_{2}}\approx 0.605(15).$

Thus, for example, three recent independent DMRG calculations yielded the values $\kappa_{c_{1}}\approx 0.41$, $\kappa_{c_{2}}\approx 0.62$ Jiang et al. (2012), $\kappa_{c_{1}}\approx 0.44$, $\kappa_{c_{2}}\approx 0.61$ Gong et al. (2014), and $\kappa_{c_{1}}\approx 0.46$, $\kappa_{c_{2}}\approx 0.62$ Wang and Sandvik (2018), while two high-order CCM calculations yielded the values $\kappa_{c_{1}}\approx 0.447$, $\kappa_{c_{2}}\approx 0.586$ Götze et al. (2012) and $\kappa_{c_{1}}\approx 0.454$, $\kappa_{c_{2}}\approx 0.588$ Richter et al. (2015). Similar results have also been found, for example, from a plaquette-renormalized tensor-network study Yu and Kao (2012) that gave values $\kappa_{c_{1}}\approx 0.40$, $\kappa_{c_{2}}\approx 0.62$; a renormalization group (RG) approach Metavitsiadis et al. (2014), in which the RG flows were numerically integrated, that gave values $\kappa_{c_{1}}\approx 0.416$, $\kappa_{c_{2}}\approx 0.606$; a cluster mean-field theory approach Ren et al. (2014) that gave values $\kappa_{c_{1}}\approx 0.42$, $\kappa_{c_{2}}\approx 0.59$; a VMC calculation using an AFM fermionic RVB class of trial wave functions Chou and Chen (2014) that gave values $\kappa_{c_{1}}\approx 0.45$, $\kappa_{c_{2}}\approx 0.6$; and a separate many-variable VMC calculation combined with quantum-number projections Morita et al. (2015) that gave values $\kappa_{c_{1}}\approx 0.4$, $\kappa_{c_{2}}\approx 0.6$.

We should note, however, that while there is broad agreement on the value for $\kappa_{c_{2}}$, there are still outlier calculations for $\kappa_{c_{1}}$. For example, a bond-operator formalism approach that included cubic and quartic interactions beyond the harmonic approximation Doretto (2014) yielded a lower value of $\kappa_{c_{1}}\approx 0.34$ (and $\kappa_{c_{2}}\approx 0.59$), while a recent approach using the cluster update algorithm for tensor product states Wang et al. (2016) yielded the much higher value of $\kappa_{c_{1}}\approx 0.572$. In this context it is interesting to note too that a large-scale ED calculation using finite-size scaling to the thermodynamic limit ($N\to\infty$) on finite square lattices of up to $N=40$ sites Richter and Schulenburg (2010) yielded values $\kappa_{c_{1}}\approx 0.35$, $\kappa_{c_{2}}\approx 0.66$ based on the points where the Néel and striped order parameters vanish, respectively, but also gave values $\kappa_{c_{1}}\approx 0.46$, $\kappa_{c_{2}}\approx 0.60$ based on points where the respective zero-field transverse (uniform) magnetic susceptibility vanishes. The latter estimates are clearly in much better agreement with the modern consensual values. On the other hand we should note that the vanishing of the magnetic susceptibility only denotes the opening up of a new gapped phase. Any non-magnetic gapless state (e.g., of the QSL variety) would not, of course, be seen by calculations of the susceptibility alone.

Turning our attention now to the bilayer, it is clear that the interlayer $J_{1}^{\perp}$ bonds have no additional frustrating effect on the intralayer magnetic LRO. Indeed, in the classical limit ($s\to\infty$) they have zero effect. However, for finite spin quantum numbers $s$, if we consider first the case of zero frustration ($\kappa=0$), the $J_{1}$ and $J_{1}^{\perp}$ bonds do still compete with one another since the $J_{1}^{\perp}$ bonds by themselves promote the formation of NN interlayer dimers. For the present spin-$\frac{1}{2}$ case when $J_{1}^{\perp}>0$ these are spin-singlet pairs, while for $J_{1}^{\perp}<0$ they are spin-triplet pairs. Thus, even with zero frustration ($J_{2}=0$), the introduction of AFM $J_{1}^{\perp}$ bonds induces a competition between a GS magnetic phase with Néel LRO and a nonclassical paramagnetic phase of the VBC kind, which is formed of interlayer dimers. The resulting spin-$\frac{1}{2}$ $J_{1}$–$J_{1}^{\perp}$ model on a square-lattice bilayer has been studied previously Hida (1990, 1992); Millis and Monien (1993, 1994); Sandvik and Scalapino (1994); Sandvik et al. (1995); Chubukov and Morr (1995); Weihong (1997); Shevchenko and Sushkov (1999); Shevchenko et al. (2000); Wang et al. (2006); Collins and Hamer (2008); Fritz et al. (2011); Ganesh et al. (2011); Helmes and Wessel (2014); Devakul and Singh (2014); Lohöfer et al. (2015). Since QMC calculations can be performed in this case (i.e., when $\kappa=0$), the position $\delta_{c_{1}}^{>}(\kappa=0)$ of the QPT between the Néel-ordered state and the quantum disordered interlayer-dimer VBC (IDVBC) state can be ascertained with high accuracy. For example, a finite-size scaling of QMC results on lattices with $2L^{2}$ spins with $L\leq 10$ Sandvik and Scalapino (1994) gave a value $\delta_{c_{1}}^{>}(\kappa=0)=2.51(2)$, while a more recent QMC calculation of Wang et al. Wang et al. (2006) using the improved stochastic series-expansion algorithm with operator-loop updates and finite-size scaling on $L\times L\times 2$ lattices with $L\leq 42$ gave the very precise value $\delta_{c_{1}}^{>}(\kappa=0)=2.5220(1)$. An exponent-biased SE analysis of Zheng Weihong (1997) gave the comparable result $\delta_{c_{1}}^{>}(\kappa=0)=2.537(5)$.

We turn now finally to the case of interest here where we also introduce intralayer frustration via the NNN AFM $J_{2}$ bonds. The resulting spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on a square-lattice bilayer has received much less attention Hida (1996, 1998) than either of the limiting cases $\delta=0$ Chandra and Doucot (1988); Dagotto and Moreo (1989); Gelfand et al. (1989); Sachdev and Bhatt (1990); Chubukov and Jolicoeur (1991); Read and Sachdev (1991); Richter (1993); Richter et al. (1994); Ivanov and Richter (1994); Schulz et al. (1996); Oitmaa and Weihong (1996); Zhitomirsky and Ueda (1996); Trumper et al. (1997); Bishop et al. (1998); Singh et al. (1999); Kotov et al. (1999); Capriotti and Sorella (2000); Capriotti et al. (2001); Takano et al. (2003); Roscilde et al. (2004); Lante and Parola (2006); Sirker et al. (2006); Schmalfuß et al. (2006); Mambrini et al. (2006); Bishop et al. (2008a, b); Darradi et al. (2008); Isaev et al. (2009); Murg et al. (2009); Ralko et al. (2009); Richter and Schulenburg (2010); Reuther and Wölfle (2010); Reuther et al. (2011); Yu and Kao (2012); Götze et al. (2012); Jiang et al. (2012); Mezzacapo (2012); Li et al. (2012); Wang et al. (2013); Zhang and Beach (2013); Hu et al. (2013); Gong et al. (2014); Doretto (2014); Qi and Gu (2014); Metavitsiadis et al. (2014); Ren et al. (2014); Wang (2014); Chou and Chen (2014); Morita et al. (2015); Richter et al. (2015); Wang et al. (2016); Poilblanc and Mambrini (2017); Haghshenas and Sheng (2018); Yu et al. (2018); Wang and Sandvik (2018); Liu et al. (2018) or $\kappa=0$ Hida (1990, 1992); Millis and Monien (1993, 1994); Sandvik and Scalapino (1994); Sandvik et al. (1995); Chubukov and Morr (1995); Weihong (1997); Shevchenko and Sushkov (1999); Shevchenko et al. (2000); Wang et al. (2006); Collins and Hamer (2008); Fritz et al. (2011); Ganesh et al. (2011); Helmes and Wessel (2014); Devakul and Singh (2014); Lohöfer et al. (2015) discussed above. We note, however, that a very recent paper Stapmanns et al. (2018) studied the case where frustration is introduced instead via an interlayer NNN AFM $J_{2}^{\perp}$ bond, resulting in a $J_{1}$–$J_{1}^{\perp}$–$J_{2}^{\perp}$ model, with very different properties and behavior (and see also Ref. Alet et al. (2016)).

Before presenting our results in Sec. IV for the $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model it may be worthwhile to outline first our aims and expectations. As the interlayer coupling parameter $\delta$ is turned on we anticipate that its initial effect at a fixed value of $\kappa$ at which magnetic order of either the Néel or striped sort exists, will be to enhance the stability of the corresponding quasiclassical state, since the effect of the $J_{1}^{\perp}$ bonds is to increase the number of NN bonds and hence to take a step towards three-dimensionality. A priori, one can expect that this effect is roughly symmetric with respect to small positive and negative values of $\delta$. Hence, we anticipate that on a $\kappa\delta$ plot the Néel and striped phase boundaries will both show a cusp at $\delta=0$. Thus, we expect that $\kappa_{c_{1}}(\delta)$ will initially increase and $\kappa_{c_{2}}(\delta)$ will initially decrease as $\delta$ is either increased to small positive values or decreased to small negative values. Accordingly, we are particularly interested in what then happens as $|\delta|$ is increased further in both the FM ($\delta<0$) and AFM ($\delta>0$) regimes of interlayer coupling. The situation is expected to be very different in the two cases.

Thus, firstly, in the region where $\delta<0$ it is evident that in the limit $\delta\to-\infty$ the system will simply behave as a spin-1 $J_{1}$–$J_{2}$ model on a square-lattice monolayer. Unlike the corresponding spin-$\frac{1}{2}$ case the spin-1 $J_{1}$–$J_{2}$ model on the square-lattice seems to show Bishop et al. (2008c, d); Haghshenas et al. (2018) a direct transition between the Néel and striped phases at a critical value $\kappa\approx 0.55$, although an early DMRG calculationJiang et al. (2009) indicated a disordered paramagnetic phase in the narrow region $0.525\lesssim\kappa\lesssim 0.555$. Interestingly, a later and very recent DMRG calculation Haghshenas et al. (2018) using larger finite lattices showed that if such an intermediate region did exist it could do so only in the much smaller regime $0.545\lesssim\kappa\lesssim 0.550$. While the system sizes in the DMRG calculations Haghshenas et al. (2018); Jiang et al. (2009) were too small for a critical analysis, both the CCM analysis Bishop et al. (2008c, d) and an infinite projected entangled-pair state analysis Haghshenas et al. (2018), have shown that the direct transition between the Néel and striped phases for the spin-1 case is a first-order transition. The best estimate for the critical coupling of the transition is $\kappa\approx 0.549$ Haghshenas et al. (2018).

Returning to our bilayer model, let us denote by $\delta_{c_{1}}^{{\mathrm{F}}}(\kappa)$ and $\delta_{c_{2}}^{{\mathrm{F}}}(\kappa)$ the critical values of $\delta$ (for a given value of $\kappa$) at which Néel order and striped order, respectively, melt in the regime of FM interlayer coupling ($\delta<0$). Equivalently, these phase boundaries, $\delta=\delta_{c_{1}}^{{\mathrm{F}}}(\kappa)$ and $\delta=\delta_{c_{2}}^{{\mathrm{F}}}(\kappa)$, are also denoted, respectively as $\kappa=\kappa_{c_{1}}^{{\mathrm{F}}}(\delta)$ and $\kappa=\kappa_{c_{2}}^{{\mathrm{F}}}(\delta)$. In the light of the above discussion it seems clear that in the half-plane $\delta<0$ there must exist a quantum triple point (QTP) that occurs at a value $\delta=\delta_{{\mathrm{T}}}^{{\mathrm{F}}}$ such that $\kappa_{c_{1}}^{{\mathrm{F}}}(\delta_{{\mathrm{T}}}^{{\mathrm{F}}})=\kappa_{c_{2}}^{{\mathrm{F}}}(\delta_{{\mathrm{T}}}^{{\mathrm{F}}})$ or, equivalently, when $\delta_{c_{1}}^{{\mathrm{F}}}(\kappa_{{\mathrm{T}}}^{{\mathrm{F}}})=\delta_{c_{2}}^{{\mathrm{F}}}(\kappa_{{\mathrm{T}}}^{{\mathrm{F}}})$. Thus, if the position of this QTP is ($\kappa_{{\mathrm{T}}}^{{\mathrm{F}}},\delta_{{\mathrm{T}}}^{{\mathrm{F}}}$), then for all values $\delta<\delta_{{\mathrm{T}}}^{{\mathrm{F}}}$ there will be a direct transition between the Néel and striped phases at a value $\kappa^{{\mathrm{F}}}(\delta)$, where we expect $\lim_{\delta\to-\infty}\kappa^{{\mathrm{F}}}(\delta)\approx 0.549$. One of our aims will be to evaluate accurately the position ($\kappa_{{\mathrm{T}}}^{{\mathrm{F}}},\delta_{{\mathrm{T}}}^{{\mathrm{F}}}$) of the QTP where the Néel, striped, and disordered paramagnetic phases meet in the half-plane $\delta<0$. It seems almost certain that $\kappa_{{\mathrm{T}}}^{{\mathrm{F}}}$ will lie between the monolayer values $\kappa_{c_{1}}(\delta=0)\approx 0.43(3)$ and $\kappa_{c_{2}}(\delta=0)\approx 0.605(15)$.

The possible scenarios in the half-plane $\delta>0$ are even more interesting. One possibility is the obvious analog to that discussed above, with another QTP between the Néel, striped, and intermediate paramagnetic phases. However, in this scenario, such a QTP would presumably have to be accompanied by another QTP, at a larger value of $\delta$, now between the Néel and striped phases together with the gapped IDVBC state that we know must physically occur for large enough values of $\delta$ at any fixed value of $\kappa$. Such a scenario (at least as far as the first QTP is concerned) was obtained in an earlier CCM calculation Schmalfuß et al. (2006) of the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the stacked square lattice (i.e., the same model as considered here but with an infinite number of layers in which all NN interlayer pairs are connected via $J_{1}^{\perp}$ bonds). In this case, of course, no IDVBC state occurs, and the Néel and striped AFM phases simply undergo a first-order transition for all values of the coupling parameter $\delta$ beyond the first (and only) QTP in this case.

An alternative, perhaps more intriguing, scenario in the half-plane $\delta>0$ is one in which the boundaries of the two quasiclassical AFM phases turn back on themselves, in a reentrant fashion, sufficiently rapidly as $\delta$ is increased so that they avoid crossing. One of our major aims here is to perform sufficiently accurate calculations as to be able to distinguish with confidence between such different scenarios. Before we present our findings in Sec. IV, however, we first briefly discuss in Sec. III the most important features of the CCM that we use to obtain them.

III THE COUPLED CLUSTER METHOD

The CCM Coester (1958); Coester and Kümmel (1960); Čižek (1966); Kümmel et al. (1978); Bishop and Lührmann (1978, 1982); Arponen (1983); Bishop and Kümmel (1987); Arponen et al. (1987a, b); Bartlett (1989); Arponen and Bishop (1991); Bishop (1991, 1998); Zeng et al. (1998); Farnell and Bishop (2004); Bartlett and Musiał (2007); Bishop et al. (2014) is widely recognized as providing one of the most flexible, most widely utilized, and most accurate of all available ab initio techniques in modern microscopic quantum many-body theory. One of the keys to its success is the fact that it preserves size-extensivity and size-consistency at every level of approximation, thereby enabling it to be implemented from the very outset in the thermodynamic ($N\to\infty$) limit. Hence any errors associated with finite-size scaling, as needs to be performed in almost all competing methods, are obviated. A second key to the success of the CCM lies in the fact that it also exactly preserves at all levels of approximation the very important Hellmann-Feynman theorem as well as the Goldstone linked-cluster theorem. A third key to its success is that there exist well-defined, systematic, and very widely tested hierarchies of truncations within which the method can be computationally implemented to very high orders of approximation, as will be done here. Since the CCM becomes exact within such a truncation hierarchy as the order $n$ of the approximation tends to infinity ($n\to\infty$), the only approximation ever made is in the extrapolation of such a sequence of approximants for any physical parameter calculated for the system under study. The combination of these features ensures that the CCM yields accurate and self-consistent sets of results for all GS and excited-state (ES) quantities calculated.

Amongst many applications to quantum many-body problems in fields as diverse as nuclear physics, subnuclear physics, quantum chemistry, atomic and molecular physics, quantum optics, and condensed matter physics, the CCM has, in particular, by now been applied to a wide variety of spin-lattice systems of interest in quantum magnetism (see, e.g., Refs. Bishop et al. (1998); Schmalfuß et al. (2006); Bishop et al. (2008a, b); Darradi et al. (2008); Reuther et al. (2011); Götze et al. (2012); Richter et al. (2015); Zeng et al. (1998); Farnell and Bishop (2004); Bishop et al. (2014, 2008c, 2008d); Li and Bishop (2019) and references therein). Since its application to such systems has already been widely described in the literature, therefore we content ourselves here with presenting a brief overview of only those features that are most relevant to us now.

The first step in any implementation of the CCM is to choose a suitable model (or reference) state $|\Phi\rangle$ for the $N$-body system with ($N\to\infty$) under consideration, together with a complete set of mutually commuting, multiconfigurational creation operators, $C^{+}_{I}\equiv(C_{I}^{-})^{\dagger}$. The main requirement on $|\Phi\rangle$ is that it should be a cyclic vector (or, equivalently, a generalized vacuum state) with respect to the set of operators $\{C_{I}^{+}\}$. The set-index $I$ here is used to indicate a complete labelling of the many-particle configuration created in the state $C_{I}^{+}|\Phi\rangle$. We thus require the set $\{|\Phi\rangle;C_{I}^{+}\}$ to obey the conditions,

[TABLE]

[TABLE]

[TABLE]

where $\mathbb{1}$ is the unit vector in the $N$-particle Hilbert space. It is also convenient to choose the states $\{C_{I}^{+}|\Phi\rangle\}$ that so span the $N$-body Hilbert space to be an orthonormal set,

[TABLE]

with $\delta_{I,J}$ a suitably generalized Kronecker symbol.

The exact many-body GS ket and bra states, $|\Psi\rangle$ and $\langle\tilde{\Psi}|$ ($=\langle\Psi|/\langle\Psi|\Psi\rangle$), respectively, which satisfy the respective GS Schrödinger equations,

[TABLE]

with $|\Psi\rangle$ now satisfying the intermediate normalization condition, $\langle\Phi|\Psi\rangle=1=\langle\Phi|\Phi\rangle$, together with $\langle\tilde{\Psi}|\Psi\rangle=1$, are now parametrized within the CCM with respect to the model state $|\Phi\rangle$ via the distinctive exponentiated forms of correlation operators,

[TABLE]

[TABLE]

that are one of the distinguishing features of the method. Although Hermiticity implies that the destruction correlation operator $\tilde{S}$ is formally related to its creation counterpart $S$ via the relation

[TABLE]

the CCM treats $S$ and $\tilde{S}$ as independent operators. They will clearly satisfy Eq. (9) when no approximations are made, but may violate it when truncations are made in the sums over the index $I$ in Eqs. (7) and (8), as described below in practical implementations. The compensation paid for this loss of explicit Hermiticity is the huge advantage that in all such truncations the Hellmann-Feynman theorem is now manifestly maintained.

All GS properties of the system may thus be calculated in terms of the set of real $c$-number correlation coefficients $\{{\cal S}_{I},\tilde{{\cal S}}_{I}\}$. In turn, these may be found by insertion of the parametrizations of Eqs. (7) and (8) into the respective Schrödinger equations (6), followed by projection onto the complete sets of states $\langle\Phi|C_{I}^{-}$ and $C_{I}^{+}|\Phi\rangle$, respectively. As a completely equivalent alternative procedure we may derive $\{{\cal S}_{I},\tilde{{\cal S}}_{I}\}$ by demanding that the GS energy expectation value functional, $\bar{H}=\bar{H}({\cal S}_{I},\tilde{{\cal S}}_{I})$, defined as

[TABLE]

be an extremum with respect to the entire set of parameters $\{{\cal S}_{I},{\cal\tilde{S}}_{I}\}$. By either method we may readily derive the sets of equations,

[TABLE]

[TABLE]

By using Eq. (11) we may also show that the GS energy at the stationary point may be expressed purely in terms of the set of creation coefficients $\{{\cal S}_{I}\}$ as

[TABLE]

Correspondingly, Eq. (12) may be written in the equivalent form,

[TABLE]

By contrast, the GS expectation value, $\bar{A}\equiv\langle\tilde{\Psi}|A|\Psi\rangle$, of any other operator $A$ requires both sets of GS CCM coefficients for its evaluation,

[TABLE]

We note that the characteristic CCM exponentiated operators ${\mathrm{e}}^{\pm S}$ only enter into Eqs. (11) and (12), which need to be solved for the GS coefficients $\{{\cal S}_{I},\tilde{{\cal S}}_{I}\}$, in the form of the associated similarity transform of the Hamiltonian, ${\mathrm{e}}^{-{S}}H{\mathrm{e}}^{S}$. In order to solve Eqs. (11) and (12) in practice we utilize the nested commutator expansion,

[TABLE]

where the $n$-fold nested commutators $[H,S]_{n}$ are defined iteratively as

[TABLE]

Another key feature of the CCM parametrizations of Eqs. (7) and (8) is now that the otherwise infinite sum in Eq. (16) terminates in practice for all Hamiltonians that contain only finite-order multinomials in the appropriate single-particle operators, as in the present case. The reason for this is simple, namely that all components in the expansion of Eq. (7) mutually commute by construction, as in Eq. (4). For the present Hamiltonian of Eq. (1), which is bilinear in the basic one-body operators $(s^{x}_{i,\alpha},s^{y}_{i,\alpha},s^{z}_{i,\alpha})$, the sum in Eq. (16) will terminate with the term $n=2$ for the choices for $\{|\Phi\rangle;C_{I}^{+}\}$ that we describe below, due to the basic SU(2) commutation relations. Thus, all nested commutators with $n>2$ simply vanish identically.

For the same reason, all terms in the expansion of $\bar{H}$ are linked, and it is this fact that leads to the CCM satisfying the Goldstone linked-cluster theorem (and consequently being size-extensive) at all levels of truncation in the expansions of Eqs. (7) and (8). Thus, in the solutions of Eqs. (11) and (12) the sole approximation made is in what set of configurations $\{I\}$ will be retained in the expansions of Eqs. (7) and (8), as described below.

We turn now to the choice of model state $|\Phi\rangle$ and the associated set of multiconfigurational creation operators $\{C_{I}^{+}\}$ for the present case. In any spin-lattice application of the CCM a convenient (but not the only) choice for $|\Phi\rangle$ is always any quasiclassical state with perfect magnetic LRO, i.e., one for which the spin on every lattice site is specified independently via its given spin projection onto some specified spin quantization axis. Here we will thus use both the Néel and striped AFM states as our independent CCM model states. It is highly convenient to treat all such states in the same way, so that all sites may be treated equivalently. A simple means of doing so is to choose a local spin quantization axis independently on each site (i.e., equivalently, by making a suitable passive rotation of each spin separately) so that in these local axes the reference state is a tensor product of spin-down states, $|\Phi\rangle=|\downarrow\downarrow\downarrow\cdots\downarrow\rangle$, such that all spins point along the negative $z_{s}$ direction in these local sets of axes. A beneficial effect of choosing such rotations is that all cases may henceforward be treated on an equal footing and by a universal computational code. The cases are distinguished only by that the spin Hamiltonian needs to be rewritten in terms of the particular local axes needed for each specific model state.

Such passive rotations are unitary transformations that leave the basic $\text{SU}(2)$ algebra unchanged, but also have the other beneficial effect of allowing us to choose the operators $C_{I}^{+}$ as products of single-spin raising operators $s_{k,\alpha}^{+}\equiv s_{k,\alpha}^{x}+is_{k,\alpha}^{y}$. Thus, we have that the set-index $I$ becomes a set of lattice indices, $\{I\}\to\{l_{1},l_{2},\cdots,l_{n}\,;n=1,2,\cdots,2sN\}$, where $l_{i}\equiv(k_{i},\alpha)$, and in which any given lattice site index $l_{i}$ may be repeated so that it appears no more than $2s$ times, where $s$ is the spin quantum number of the spins in the general case. Thus, we have $C_{I}^{+}\to s_{l_{1}}^{+}s_{l_{2}}^{+}\cdots s_{l_{n}}^{+}$, with $n=1,2,\cdots,2sN$. Here, we take $s=\frac{1}{2}$, so that in each configuration $I$ no lattice site may appear more than once.

In these local rotated spin axes the order parameter (i.e., the sublattice magnetization) takes the universal form,

[TABLE]

Once again, the expression ${\mathrm{e}}^{-S}s_{l_{i}}^{z}{\mathrm{e}}^{S}$ may be evaluated exactly via a nested commutator expansion akin to Eq. (16) that now terminates at the term with $n=1$.

As we have indicated above, the only approximation that we make in implementing the CCM is the truncation of the expansion of Eqs. (7) and (8) for the correlation operators $S$ and $\tilde{S}$. We shall use here the well-established localised (lattice-animal-based subsystem) LSUB$n$ scheme wherein we retain at $n$th order all multispin-flip correlations on the lattice over no more than $n$ contiguous sites. A set of sites is contiguous in this sense if every site in the set is NN to at least one other in the set, in some specified geometry. As the truncation index $n$ tends to infinity $(n\to\infty)$, the corresponding LSUB$n$ approximation becomes exact.

One may use the space- and point-group symmetries of the lattice and the particular CCM model state $|\Phi\rangle$ being used, as well as any pertinent conservation laws, to reduce the number of independent configurations retained at any order $n$ of approximation. For example, for each of the model states considered in Sec. IV (i.e., the Néel and striped states on each monolayer), the Hamiltonian of Eq. (1) conserves the total $z$ component of spin, $s_{T}^{z}\equiv\sum_{i=1}^{N}s_{l_{i}}^{z}$ (where global spin axes are now assumed), to the sector $s_{T}^{z}=0$. We denote by $N_{f}(n)$ the minimal number of distinct (and nonzero) fundamental multispin-slip configurations that are retained at a given LSUB$n$ level of approximation after all such symmetries and conservation laws are take into account. The number $N_{f}(n)$ typically still grows rapidly with $n$, and available computational resources then determine the maximum order $n$ that can be computed.

In the present case we are able to perform LSUB$n$ calculations up to the very high order $n=10$. Thus, for the spin-$\frac{1}{2}$ square-lattice bilayer model under consideration we have $N_{f}(10)=239\,021$ $(443\,813)$ when the CCM model state is chosen so that each monolayer has Néel (striped) AFM order. To derive and then to solve such larger sets of coupled nonlinear multinomial equations (11) for $S$ and linear equations (12) for $\tilde{S}$ we use both massive parallelization and large-scale supercomputer resources. In order to derive the equations (and see Ref. Zeng et al. (1998)) we also use a purpose-built and customized computer algebra package ccm , without which it would not be possible to go to such large orders $n$ of LSUB$n$ truncation.

Since no approximations have been made in the evaluation of any finite-order LSUB$n$ truncation of our basic CCM equations, nor in the subsequent evaluation of any GS parameter of the system, our only approximation is now made at the last step where we extrapolate an LSUB$n$ sequence of approximants to the (exact) $n\to\infty$ limit. By now there is a great deal of empirical evidence on how to do so. For example, for the LSUB$n$ approximants $M(n)$ to the magnetic order parameter $M$ of Eq. (18), a well-tested scheme for systems with strong frustration, and/or for which the system has a QPT between states with and without magnetic LRO, has been found (and see, e.g., Refs. Bishop et al. (2008a, b); Darradi et al. (2008); Reuther et al. (2011); Götze et al. (2012); Bishop et al. (2014, 2008c, 2008d); Li and Bishop (2019) and references cited therein) to be given by

[TABLE]

By fitting a sequence of LSUB$n$ approximants $M(n)$ to Eq. (19) we thus extract the corresponding extrapolated (LSUB$\infty$) value $\mu_{0}$ for $M$.

IV RESULTS

We show first in Fig. 2 our CCM results for the magnetic order parameter $M$ based on a model state in which each monolayer has Néel order. For values of the interlayer coupling parameter $\delta>0$ ($\delta<0$) the two layers are coupled so that NN spins are antiparallel (parallel). Results are shown as functions of $\delta$ in Figs. 2(a) and 2(b) respectively for two fixed values of the intralayer frustration parameter, $\kappa=0.3$ and $\kappa=0.5$. These values are chosen to lie on either side of the critical value $\kappa_{c_{1}}(\delta=0)\approx 0.45$, at which Néel order melts in the monolayer. In both cases results are shown for even-order LSUB$n$ approximations with $4\leq n\leq 10$, as well as for the (LSUB$\infty$) estimates $\mu_{0}$ obtained from fitting these data to the extrapolation scheme of Eq. (19). The LSUB2 data are omitted since, in principle, they are expected to be of too low order to fit well to such a three-term extrapolation scheme. Nevertheless, when separate fits are made that include them (i.e., to the LSUB$n$ data sets $n=\{2,4,6,8,10\}$), the extrapolated values for $\mu_{0}$ hardly change at all, thereby demonstrating the robustness of the fits. To enable readers to judge by eye for themselves how robust and accurate are our extrapolations, we also show in Fig. 2(a) the corresponding extrapolation (labelled LSUB$\infty$’) based on the restricted LSUB$n$ data set $n=\{6,8,10\}$. Although a fit based on only three data points to a three-term extrapolation scheme such as that of Eq. (19) is not, a priori, expected to be as robust as one based on more data points, the agreement between the LSUB$\infty$ and LSUB$\infty$’ curves can be clearly seen. Similar agreement is observed for all other values of the intralayer frustration parameter $\kappa$.

Turning first to Fig. 2(a), the LSUB$\infty$ extrapolated curve exhibits all of the features we expect from our discussion in Sec. II. Thus, firstly, the cusp at $\delta=0$ (where $M\approx 0.17$ for the value $\kappa=0.3$ shown) is exactly as expected from the observation that as $|\delta|$ is increased from zero in either direction, the order is first enhanced due to the increase in the number of NN bonds and the consequent step towards increasing the dimensionality of the system. Secondly, we observe that as $\delta$ is increased further, in the regime $\delta>0$ of AFM coupling between the two layers, $M$ attains a maximum value of about 0.29 at a value $\delta\approx 0.57$ before the effects of interlayer dimerization become sufficiently strong to start to weaken the intralayer Néel order as $\delta$ is increased further beyond that point. This continues up to an upper critical value, $\delta_{c_{1}}^{>}=\delta_{c_{1}}^{>}(\kappa)$, above which Néel order disappears entirely. For the value $\kappa=0.3$ shown in Fig. 2(a), this upper critical value is seen to be at $\delta_{c_{1}}^{>}(0.3)\approx 1.66$.

As the value of $\kappa$ is increased beyond 0.3 the cusp at $\delta=0$ in the LSUB$\infty$ curve of Fig. 2(a) is lowered until, at the value $\kappa\approx 0.447$, it reaches the $\delta=0$ axis. This is precisely the value we obtain for $\kappa_{c_{1}}(\delta=0)$ within this same extrapolation scheme, using the LSUB$n$ data set $n=\{4,6,8,10\}$. We note in passing that for the monolayer case ($\delta=0$) it has also been possible to perform LSUB12 calculations Richter et al. (2015). The corresponding value obtained from fitting to the LSUB$n$ data set $n=\{4,6,8,10,12\}$ is $\kappa_{c_{1}}(\delta=0)\approx 0.454$ Richter et al. (2015), which again demonstrates the accuracy and robustness of our extrapolations.

If we now slightly increase $\kappa$ beyond this value $\kappa_{c_{1}}(0)$, we obtain curves such as those shown in Fig. 2(b). Thus, for a certain range of values above $\kappa_{c_{1}}(0)$, for which Néel order is absent for the monolayer $(\delta=0)$, as $\delta$ is either decreased or increased from this value, Néel order becomes re-established at certain critical values, $\delta_{c_{1}}^{{\mathrm{F}}}(\kappa)$ for $\delta<0$ and $\delta_{c_{1}}^{<}(\kappa)$ for $\delta>0$. For the value $\kappa=0.5$ shown in Fig. 2(b) these values are seen to be $\delta_{c_{1}}^{{\mathrm{F}}}(0.5)\approx-0.12$ and $\delta_{c_{1}}^{<}(0.5)\approx 0.12$. Once again, as $\delta$ is now increased beyond the value $\delta_{c_{1}}^{<}(\kappa)$, Néel order is at first enhanced until $M$ attains a maximal value. For the value $\kappa=0.5$ shown in Fig. 2(b) this maximum value for $M$, for the case of AFM interlayer coupling, is about 0.11 at a value $\delta\approx 0.42$. Further increase in $\delta$ then reduce the magnetic order until it again melts entirely at a value $\delta_{c_{1}}^{>}(\kappa)$. For $\kappa=0.5$ this upper critical value is seen from Fig. 2(b) to be at $\delta_{c_{1}}^{>}(0.5)\approx 0.67$.

If we continue to increase $\kappa$ slowly beyond the value $\kappa=0.5$ shown in Fig. 2(b), the lower and upper critical values, $\delta_{c_{1}}^{<}(\kappa)$ and $\delta_{c_{1}}^{>}(\kappa)$, move towards one another until at some value $\kappa_{1}^{\mathrm{max}}\approx 0.535$ they merge, $\delta_{c_{1}}^{<}(\kappa_{1}^{\mathrm{max}})=\delta_{c_{1}}^{>}(\kappa_{1}^{\mathrm{max}})\approx 0.3$. Néel order is then wholly absent for any value $\kappa>\kappa_{1}^{\mathrm{max}}$ of the intralayer frustration parameter and for any value $\delta>0$ of the AFM interlayer coupling.

Corresponding results to those shown in Fig. 2, which are based on a CCM model state with Néel order on each monolayer, as in Fig. 1(b), are now shown in Fig. 3 based on a corresponding model state with striped AFM order on each monolayer, as in Fig. 1(c). The two layers are again coupled so that NN interlayer spins are antiparallel for $\delta>0$ and parallel for $\delta<0$. The two values of $\kappa$ shown, viz., $\kappa=0.7$ in Fig. 3(a) and $\kappa=0.57$ in Fig. 3(b) are now chosen to lie on either side of the critical value $\kappa_{c_{2}}(\delta=0)\approx 0.59$ for which striped order melts in the monolayer. In Fig. 3(a) we also show the corresponding LSUB$\infty$’ extrapolated curve based on the restricted LSUB$n$ data set $n=\{6,8,10\}$. The agreement with the LSUB$\infty$ curve based on the full LSUB$n$ data set $n=\{4,6,8,10\}$ is again observed to be excellent. Similar levels of agreement are found for all other values of $\kappa$. Once again, for the special case of the monolayer LSUB12 calculations have also been performed Richter et al. (2015) based on the striped model state. In this case we find $\kappa_{c_{2}}(0)\approx 0.587$ based on the extrapolation of Eq. (19) with the LSUB$n$ data set $n=\{4,6,8,10\}$ as input, while the corresponding result based on the set $n=\{4,6,8,10,12\}$ yielded the almost identical value $\kappa_{c_{2}}(0)\approx 0.588$ Richter et al. (2015).

The similarity between Fig. 2(a) and 3(a), and also between Fig. 2(b) and 3(b), is self-evident. In this case we find that the lower and upper critical values, $\delta_{c_{2}}^{<}(\kappa)$ and $\delta_{c_{2}}^{>}(\kappa)$, as seen in Fig. 3(b) in the region $\delta>0$ for the value $\kappa=0.57<\kappa_{c_{2}}(0)$ again move towards one another as $\kappa$ is now slowly decreased beyond this value, until they merge at a value $\kappa_{2}^{\mathrm{min}}\approx 0.555$, where $\delta_{c_{2}}^{<}(\kappa_{2}^{\mathrm{min}})=\delta_{c_{2}}^{{\mathrm{>}}}(\kappa_{2}^{\mathrm{min}})\approx 0.3$. For all values $\kappa<\kappa_{2}^{\mathrm{min}}$ of the intralayer frustration parameter striped order is then absent, whatever the value $\delta>0$ of the AFM interlayer coupling.

In Figs. 4(a) and 4(b) we now show sets of extrapolated (LSUB$\infty$) curves for the magnetic order parameter $M$ as a function of the interlayer coupling parameter $\delta$ for each layer with Néel order and striped order, respectively, for various fixed values of the intralayer frustration $\kappa$. These correspond, respectively, to curves such as those shown in Figs. 2 and 3. From Fig. 4(b) it is clear that the position of the cusp at $\delta=0$ for the striped-ordered phase rather rapidly approaches a limiting value for $M$ as $\kappa$ is increased. This corresponds to the limit ($\kappa\to\infty$) of the model where each layer corresponds to two independent, and equivalent, interpenetrating square sublattices, each of which is Néel-ordered. A recent CCM calculation Farnell et al. (2018) for the spin-$\frac{1}{2}$ Heisenberg antiferromagnet on the square lattice utilized an LSUB$n$ data set $n=\{4,6,8,10,12\}$ to yield the extrapolated value $M\approx 0.3093$. By comparison, we find here the remarkably close value $M(\delta=0)\approx 0.3094$ from using Eq. (19) with the LSUB$n$ data set $n=\{4,6,8,10\}$, and as shown in Fig. 4(b), for the curve $\kappa=1.0$, which value is itself already extremely close to the limiting value obtained as $\kappa\to\infty$.

We also note from Fig. 4(a) that for the case of zero intralayer frustration ($\kappa=0$) our extrapolation of Eq. (19) leads to an upper critical value, $\delta_{c_{1}}^{>}(\kappa=0)\approx 2.84$, of the interlayer coupling parameter, beyond which Néel order melts and a phase with IDVBC order is stabilized. This may be compared with the corresponding value of 2.5220 obtained from a large-scale QMC simulation Wang et al. (2006) of the spin-$\frac{1}{2}$ $J_{1}$–$J_{1}^{\perp}$ model discussed in Sec. II. In this context it is worth noting that the location of the phase boundary is less accurately determined in our CCM calculations for the region where the order-disorder transition is essentially driven by singlet-dimerization than in the region where it is essentially due to frustration.

Particularly in the half-plane $\delta<0$, corresponding to FM interlayer coupling, it is also convenient to investigate the magnetic order parameter $M$ for the two quasiclassical phases with AFM orderings on each monolayer as functions of the intralayer frustration parameter $\kappa$, for various fixed values of the interlayer coupling parameter $\delta$. In Fig. 5 we show a set of such curves for both AFM monolayer orderings. The curves for $\delta=0$ again exhibit the corresponding critical points, $\kappa_{c_{1}}(\delta=0)\approx 0.447$ and $\kappa_{c_{2}}(\delta=0)\approx 0.587$, for the melting of Néel and striped order, respectively, for the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$ model on the square-lattice monolayer. We also show the corresponding curves for $\delta=0.3$ in the region of AFM coupling between the layers, which is approximately the value of $\delta$ for which the phase boundaries of the two quasiclassical phases are closest together in this region $\delta>0$. Thus, at $\delta=0.3$, the paramagnet state exists only in the very narrow regime $0.535\lesssim\kappa\lesssim 0.555$ of the frustration parameter.

We also show in Fig. 5 similar curves for various values of $\delta$ in the half-plane $\delta<0$, which corresponds to FM interlayer coupling.

One sees clearly that for all values $\delta<\delta_{{\mathrm{T}}}^{{\mathrm{F}}}\approx-0.45$ the respective curves $M=M(\kappa)$ for the two quasiclassical phases cross one another at a value $\kappa^{{\mathrm{F}}}(\delta)>0$, indicating a direct first-order transition between them. The value $\kappa^{{\mathrm{F}}}(\delta)$ is seen to be almost independent of $\delta$ in this regime. For example $\kappa^{{\mathrm{F}}}(-3)\approx 0.539$, which may be compared with the expected value ${\mathrm{lim}}_{\delta\to-\infty}\kappa^{{\mathrm{F}}}(\delta)\approx 0.549$, viz., the value that corresponds to the critical coupling for the direct transition between the two states in the spin-1 $J_{1}$–$J_{2}$ model on the square-lattice monolayer Haghshenas et al. (2018), as discussed in Sec. II. From Fig. 5 we see that the QTP where the Néel, striped, and paramagnetic phases meet is situated at $(\kappa_{{\mathrm{T}}}^{{\mathrm{F}}},\delta_{{\mathrm{T}}}^{{\mathrm{F}}})\approx(0.547,-0.45)$.

Of course, in view of the observation that the direct transition between the two quasiclassical phases for $\delta\lesssim-0.45$ is of first-order type, we can also corroborate our results in this regime by using the fact that the GS energies of the two states should also cross at the phase boundary. We use the appropriate extrapolation scheme for the LSUB$n$ approximants $e(n)$ to the GS energy per spin, $e\equiv E/N$, which is well known to be given by

[TABLE]

We may thus make use of Eq. (20) with the same LSUB$n$ input data sets $n=\{4,6,8,10\}$ as used in Fig. 5 to corroborate the points shown in Fig. 5 by cross ($\times$) symbols, which denote the points where the respective $M=M(\kappa)$ curves for the two quasiclassical phases cross one another for a given value of $\delta$ at the value $\kappa=\kappa^{\mathrm{F}}(\delta)$. We find, for example, at the value $\delta=-3$, the two corresponding LSUB$\infty$ extrapolated curves $e=e(\kappa)$ cross one another at the value $\kappa\approx 0.558$, which may be compared with the value $\kappa^{\mathrm{F}}(-3)\approx 0.539$ cited above, where the respective $M=M(\kappa)$ curves cross. The agreement between these two essentially independent calculations is good. This clearly also provides internal confirmation of the robustness of the extrapolation schemes that form the sole approximation made in our results

By combining results from curves such as those shown in Figs. 2–5 to extract the points where the extrapolated (LSUB$\infty$) GS magnetic order parameter $M$ vanishes for the two quasiclassical AFM phases, we may finally construct the zero-temperature ($T=0$) quantum phase diagram of our spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$–$J_{1}^{\perp}$ model on the bilayer square lattice. It is shown in Fig. 6 in the $\kappa$–$\delta$ half-plane with $\kappa>0$.

We note that different symbols are used in Fig. 6 to distinguish between points on the phase boundaries that have been extracted from calculations done at fixed values of $\kappa$ and those extracted from calculations done at fixed values of $\delta$. It is clear by visual inspection that these two sets of critical points lie very accurately on a smooth boundary curve for each collinear AFM state. This again provides good internal evidence that our extrapolations are robust and accurate.

Our results are now discussed and summarized in Sec. V.

V DISCUSSION AND SUMMARY

We see from Fig. 6 that the phase boundary for each of the AFM quasiclassical states is rather accurately linear in the $\kappa$–$\delta$ plane for large enough values of AFM interlayer coupling (viz., for $\delta\gtrsim 1.5$). For the Néel phase boundary the slope of the curve at large values of $\delta$ is ${\mathrm{d}}\delta/{\mathrm{d}}\kappa\approx-3.9$, while the corresponding value for the striped phase boundary is ${\mathrm{d}}\delta/{\mathrm{d}}\kappa\approx 3.7$. If this linear behavior would continue unchanged to smaller values of $\delta$ the two curves would cross at a point $(\kappa,\delta)\approx(0.54,0.7)$. Instead, as we see from Fig. 6, the two curves turn against each other before this point, in a typical “avoided crossing” manner, although they do approach one another rather closely in the region around $\delta\approx 0.3$. This has the effect that the entire disordered paramagnetic regime in Fig. 6 is singly connected.

At small values of $\delta$ both phase boundaries exhibit the reentrant behavior expected from our discussion in Sec. II, with both displaying cusps at $\delta=0$. However, a closer inspection of Fig. 6 shows that the nature of the cusps is quite different for the two phase boundaries, with that on the Néel side much “sharper” than its counterpart on the striped side. More quantitatively, on the striped side the slopes of the curves on each side of the cusp at the point $\kappa=\kappa_{c_{2}}(0),\delta=0$ are clearly nonzero. By contrast, on the Néel side the corresponding slopes of the curves of the cusp at the point $\kappa=\kappa_{c_{1}}(0),\delta=0$ appear to be zero (or very close) to zero. This difference would certainly explain why the critical parameter $\kappa_{c_{1}}(0)$ is more difficult to calculate accurately than the corresponding parameter $\kappa_{c_{2}}(0)$ for the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$ model on the square lattice, as has been discussed previously in Sec. II. The difference is surely also a direct reflection of the different natures of the two QPTs in that model. Thus, our results are in clear accord with the consensual view that, while the transition at $\kappa_{c_{2}}(0)$ between the striped and paramagnetic phases is a first-order one, that at $\kappa_{c_{1}}(0)$ between the Néel and paramagnetic phases is continuous.

In the half-plane $\delta<0$ the phase boundaries of the two quasiclassical AFM phases end at a QTP where they meet a line of direct first-order transitions between the two phases. Starting at the QTP, which we have calculated as being located at $(\kappa_{\mathrm{T}}^{\mathrm{F}}\approx 0.547,\delta_{\mathrm{T}}^{\mathrm{F}}\approx-0.45)$, this line of first-order transitions is very nearly a vertical straight line in the $\kappa$–$\delta$ plane. At large negative values of $\delta$ it approaches the value $\kappa^{\mathrm{F}}(\delta\to-\infty)\approx 0.539$, which is itself an accurate estimate of the QCP for the direct transition between the two quasiclassical phases of the spin-1 $J_{1}$–$J_{2}$ model on the square lattice.

It is completely beyond the scope of the present paper to investigate in detail the nature of the phases in the (grey shaded) paramagnetic regime in Fig. 6, outside the respective regimes in which we have calculated that the Néel state or the striped state on each monolayer forms the stable GS phase. Nevertheless we conclude with a few comments on this issue.

We have already discussed in Sec. I the lack of overall consensus for the nature of the GS phase or phases in the region $\kappa_{c_{1}}(0)<\kappa<\kappa_{c_{2}}(0)$ for the spin-$\frac{1}{2}$ $J_{1}$–$J_{2}$ model on the square-lattice monolayer. However, the results of two independent high-order techniques, viz., the CCM Darradi et al. (2008); Richter et al. (2015) and the DMRG method Gong et al. (2014), applied to the model, are both compatible with the existence of two phases in this paramagnetic region. Both CCM and DMRG calculations also agree on the critical values $\kappa_{c_{1}}(0)\approx 0.45(1)$ and $\kappa_{c_{2}}(0)\approx 0.60(1)$. Furthermore, both can be consistently interpreted with the hypothesis of a gapped plaquette-ordered VBC (PVBC) ground state in the region $0.5\lesssim\kappa<\kappa_{c_{2}}(0)$, and a ground state in the region $\kappa_{c_{1}}(0)<\kappa\lesssim 0.5$ that could be a gapless QSL state. In view of the single-connectedness of the entire paramagnetic regime that we have found, it is clear that, based on the above scenario being true for $\delta=0$, this paramagnetic regime for the bilayer should include at least three phases, viz., QSL, PVBC, and IDVBC. It will clearly be of great future interest to study the boundaries of these phases in detail.

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.