text:
# Quantum Phase Transitions of a Square-Lattice Heisenberg Antiferromagnet with Two Kinds of Nearest-Neighbour Bonds: A High-Order Coupled-Cluster Treatment

## I Introduction

The subject of quantum spin-half antiferromagnetism in low-dimensional systems has attracted a great deal of interest in recent times in connection with the magnetic properties of the cuprate high-temperature superconductors. However, low-dimensional quantum spin systems are of interest in their own right as examples of strongly interacting quantum many-body systems. Although we know from the Mermin-Wagner theorem  that thermal fluctuations are strong enough to destroy magnetic long-range order at any finite temperature, the role of quantum fluctuations is less understood. As a result of intensive work in the late eighties, it is now well-established that the ground-state of the Heisenberg antiferromagnet on the square lattice with nearest-neighbour interactions is long-range ordered (see for example the review in Ref. ). However, Anderson’s and Fazekas’ investigations  of the triangular lattice led to a conjecture that quantum fluctuations plus frustration may be sufficient to destroy the Néel-like long-range order in two dimensions. Another specific area of recent research is the spin-half $`J_1`$$`J_2`$ antiferromagnet on the square-lattice where the frustrating diagonal $`J_2`$ bonds plus quantum fluctuations are able to realize a second-order transition from Néel ordering to a disordered quantum spin liquid (see for example Refs.  and references therein). On the other hand, there are cases in which frustration causes a first-order transition in quantum spin systems in contrast to a second-order transition in the corresponding classical model (see for example Refs. ).

Besides frustration, there is another mechanism to realize the “melting” of Néel ordering in the ground states of unfrustrated Heisenberg antiferromagnets, namely the formation of local singlet pairs of two coupled spins. This mechanism may be relevant for the quantum disordered state in bilayer systems  as well as in CaV4O9 (see for example Refs.  and references therein). The formation of local singlets is connected with a gap in the excitation spectrum. By contrast, the opening of a gap in the excitation spectrum of frustrated systems seems to be less clear and might be dependent on details of the exchange interactions.

In the present paper, we study a model which contains both mechanisms, frustration and singlet formation, in different parameter regions. We mainly use in this article the coupled cluster method, which has become widely recognized as one of the most powerful and most universal techniques in quantum many-body theory. In recent years there has been increasing success in applying the CCM to quantum spin systems, especially with the advent of high-order approximations which utilize computer algebra. Subsequently, high-order CCM approximations have been applied to the $`XXZ`$ model, the anisotropic $`XY`$ model and the $`J_1`$$`J_2`$ model. In addition to the CCM results we also present variational, spin-wave theory (SWT) and exact diagonalization (ED) results for the sake of comparison.

## II The model

We consider a spin-half Heisenberg model on a square lattice with nearest neighbor bonds $`J`$ and $`J^{}`$ in a regular zigzag pattern as shown in Fig. 1. The Hamiltonian is given by

$`H`$ $`=`$ $`J{\displaystyle \underset{<ij>_1}{}}\text{s}_i\text{s}_j+J^{}{\displaystyle \underset{<ij>_2}{}}\text{s}_i\text{s}_j`$ (1)
$`\stackrel{J=1}{=}`$ $`{\displaystyle \underset{iA}{}}{\displaystyle \underset{p}{}}(1+\delta _{p,p_J^{}}(J^{}1))\text{s}_i\text{s}_{i+p}.`$ (2)
The sums over $`<ij>_1`$, and $`<ij>_2`$ represent sums over the nearest-neighbour bonds shown in Fig. 1 by dashed and solid lines, respectively. Throughout the paper we fix the $`J`$ bond to be antiferromagnetic ($`J>0`$) and henceforth scale it to the value $`J=1`$, and consider $`J^{}`$ as the free parameter of the model. We also split the square lattice into the equivalent $`A`$ and $`B`$ sublattices shown in Fig. 1. In Eq. (1) the sum over $`i`$ runs over the sites of the sublattice $`A`$, with vectors $`p=\{(0,\pm 1),(\pm 1,0)\}`$ connecting nearest neighbours. In particular, $`p_J^{}=(1,0)`$ pertains to the coupling with $`J^{}`$ bonds.

Each square lattice plaquette consists of three $`J=1`$ bonds and one $`J^{}`$ bond. In the case of ferromagnetic $`J^{}`$ bonds (i.e., $`J^{}<0`$), the plaquettes are frustrated. Conversely, for antiferromagnetic $`J^{}`$ bonds (i.e., $`J^{}>0`$) there is no frustration in the system, although the difference of the coupling strengths $`J`$ and $`J^{}`$ leads to quantum competition. This model has been treated previously using perturbation theory, renormalized spin wave theory (RSWT) and exact diagonalization (ED). It allows us to study the influence of local singlet formation ($`J^{}>1`$) and frustration ($`J^{}<0`$) on the stability of the Néel order within a single model.

Ferromagnetic bonds in an antiferromagnetic matrix have been discussed in recent times  in connection with the proposal by Aharony and coworkers  to model localized oxygen holes in the Cu-O-planes by local ferromagnetic bonds between the copper spins. It was argued that random ferromagnetic bonds may influence the antiferromagnetic order drastically and may support the realization of a quantum spin-liquid state.

On the other hand, the case of antiferromagnetic $`J^{}`$ bonds with $`J^{}>1`$ resembles the situation in bilayer systems and in the depleted square-lattice antiferromagnet CaV4O9, in which the competition between two different antiferromagnetic bonds leads to a phase transition from antiferromagnetic long-range order to quantum disorder with a finite gap. It is seen in this article that the transition point obtained for the model of Eq. (1) is quite close to that obtained for the bilayer model.

There are some special cases of the model Hamiltonian of Eq. (1):
(i) $`J^{}=1`$: square-lattice antiferromagnet, for which the ground state is long-range ordered;
(ii) $`J^{}=0`$: honeycomb-lattice antiferromagnet, for which the ground state is long-range ordered;
(iii) $`J^{}=\mathrm{}`$: spin-$`1`$ triangular lattice, for which the ground state is long-range ordered; and
(iv) $`J^{}=+\mathrm{}`$: valence-bond solid, for which the ground state is a rotationally invariant quantum dimer state with an excitation gap.

Classical ground state. For $`J^{}>1/3`$ the Néel state is the classical ground state of the Hamiltonian of Eq. (1). At $`J_c^{}=1/3`$ there is classically a second-order phase transition to a ground state of helical nature (see Fig. 1), with a characteristic pitch angle $`\mathrm{\Phi }=\pm |\mathrm{\Phi }_{\mathrm{cl}}|`$ given by

$$|\mathrm{\Phi }_{\mathrm{cl}}|=\{\begin{array}{cc}0\hfill & J^{}>\frac{1}{3}\hfill \\ \mathrm{arccos}\left(\frac{1}{2}\sqrt{1\frac{1}{J^{}}}\right)\hfill & J^{}\frac{1}{3}\hfill \end{array}$$
(3)
where the different signs correspond to the two chiralities  of this helical state. Note that for $`\mathrm{\Phi }=0`$ this is just the Néel state. More generally, the pitch angle varies with $`J^{}`$ from $`|\mathrm{\Phi }_{\mathrm{cl}}|=0`$ for $`J^{}>1/3`$ to $`|\mathrm{\Phi }_{\mathrm{cl}}|=\pi /3`$ for $`J^{}=\mathrm{}`$. Note that $`|\mathrm{\Phi }_{\mathrm{cl}}|=\pi /3`$ (realized at $`J^{}=\mathrm{}`$) corresponds to the ground state of the spin-1 triangular lattice. We describe the directions of the spins $`\text{s}_A`$ and $`\text{s}_B`$, belonging to the $`A`$ and $`B`$ sublattices respectively, for the classical helical state with a characteristic angle $`\mathrm{\Phi }`$ as follows (and see Fig. (1)),

$`\text{s}_A(𝐑)`$ $`=`$ $`\widehat{𝐮}\mathrm{cos}𝐐𝐑+\widehat{𝐯}\mathrm{sin}𝐐𝐑,`$ (4)
$`\text{s}_B(𝐑+\widehat{x})`$ $`=`$ $`\widehat{𝐮}\mathrm{cos}(𝐐𝐑+\pi +3\mathrm{\Phi })+\widehat{𝐯}\mathrm{sin}(𝐐𝐑+\pi +3\mathrm{\Phi }),`$ (5)
where $`\widehat{𝐮}`$ and $`\widehat{𝐯}`$ are perpendicular unit vectors in the spin space, $`𝐑`$ runs over the sites of the sublattice $`A`$, and we have $`𝐐=(2\mathrm{\Phi },0)`$ for the pitch vector $`𝐐`$. We note that this general helical state does not have a periodicity in the $`x`$-direction because $`\mathrm{\Phi }`$ is in general not of the form $`m\pi /n`$ with $`m`$ and $`n`$ integral. We also note that we have only three different angles between nearest-neighbour spins, namely $`\pm (\pi \mathrm{\Phi })`$ for the $`J=1`$ couplings and $`\pi 3\mathrm{\Phi }`$ for the coupling with $`J^{}`$.

The maximum frustration is in the region around $`J^{}1`$. Bearing in mind the situation for the $`J_1`$$`J_2`$ model, one might expect that for the extreme quantum case (spin-half) quantum fluctuations might be able to open the window to a spin-liquid phase for a finite range of parameters around this region of maximum frustration. On the other hand, for strong antiferromagnetic bonds ($`J^{}1`$) there is, of course, no indication in the classical model for the breakdown of the Néel order.

Simple variational ansatz for the quantum ground-state. In the quantum case, the region of strong antiferromagnetic $`J^{}`$ bonds ($`J^{}1`$) is characterized by a tendency to singlet pairing of the two spins corresponding to a $`J^{}`$ bond. Using a high-order series expansion  the Néel order was found to be stable up to a critical value $`J_s^{}2.56`$.

A comparable value can be obtained using a simple variational wave function similar to that used  for bilayer systems, namely

$$|\mathrm{\Psi }_{\mathrm{var}}=\underset{iA}{}\frac{1}{\sqrt{1+t^2}}[|_i_{i+\widehat{x}}t|_i_{i+\widehat{x}}],$$
(6)
where the lattice sites $`i`$ and $`i+\widehat{x}`$ correspond to a $`J^{}`$ bond, and where the product in Eq. (6) is thus effectively taken over the $`J^{}`$ bonds of the lattice of Eq. (1). The trial function depends on the variational parameter $`t`$ and interpolates between a valence-bond state realized for $`t=1`$ and the Néel state for $`t=0`$. For $`t=1`$, the singlet pairing is complete and $`|\mathrm{\Psi }_{\mathrm{var}}`$ represents an eigenstate of the model of Eq. (1) in the limit $`J^{}\mathrm{}`$ (dimer state). By minimizing $`\mathrm{\Psi }_{\mathrm{var}}|H|\mathrm{\Psi }_{\mathrm{var}}`$ with respect to the variational parameter $`t`$ we get an upper bound for the ground-state energy per spin of the model of Eq. (1),

$$E_{\mathrm{var}}/N=\{\begin{array}{cc}(J^2+3J^{}+9)/24\hfill & J^{}3\hfill \\ 3J^{}/8\hfill & J^{}>3\hfill \end{array}.$$
(7)
The relevant order parameter describing the Néel order is

$$M_{\mathrm{var}}=\mathrm{\Psi }_{\mathrm{var}}|s_i^z|\mathrm{\Psi }_{\mathrm{var}}=\{\begin{array}{cc}1/2\sqrt{1J^2/9}\hfill & J^{}3\hfill \\ 0\hfill & J^{}>3\hfill \end{array},$$
(8)
showing a breakdown of the Néel order at a critical value $`J_s^{}=3`$.

## III Coupled cluster calculations

### A The ground-state formalism

The starting point for any CCM calculation (see overwiev in Ref. ) is the choice of a normalized model or reference state $`|\mathrm{\Phi }`$, together with a set of mutually commuting multispin creation operators $`C_I^+`$ which are defined over a complete set of many-body configurations $`I`$. The operators $`C_I`$ are the multispin destruction operators and are defined to be the Hermitian adjoints of the $`C_I^+`$. We choose $`\{|\mathrm{\Phi };C_I^+\}`$ in such a way that we have $`\mathrm{\Phi }|C_I^+=0=C_I|\mathrm{\Phi }`$, $`I0`$, where, by definiton, $`C_0^+=𝟙`$, the identity operator.

For spin systems, an appropriate choice for the CCM model state $`|\mathrm{\Phi }`$ is often a classical spin state, in which the most general situation is that each spin can point in an arbitrary direction. For the case of the Hamiltonian of Eq. (1), we choose the helical state illustrated in Fig. 1 to be our model state. Although the classical ground state of Eq. (1) is precisely of this form, we do not choose the classical result for the pitch angle $`\mathrm{\Phi }`$ but we consider it rather as a free parameter in the CCM calculation.

In order to perform a CCM calculation, we would like to treat each site equivalently and we do this by performing a rotation of the local spin axes at each site about the $`y`$-axis such that all spins in the model state align in the same direction, say down (along the negative $`z`$-axis). After this transformation we have

$$|\mathrm{\Phi }=|\mathrm{}\mathrm{};C_I^+=s_r^+,s_r^+s_r^{}^+,s_r^+s_r^{}^+s_{r^{\prime \prime }}^+,\mathrm{},$$
(9)
(where the indices $`r,r^{},r^{\prime \prime },\mathrm{}`$ denote any lattice site) respectively, for the model state and the multispin creation operators, which now consist of spin-raising operators only.

In order to make the spin $`\text{s}_i`$ point down let us suppose we need to perform such a rotation of the spin axes by an angle $`\delta _i`$. This is equivalent to the transformation

$`s_i^x`$ $``$ $`\mathrm{cos}\delta _is_i^x+\mathrm{sin}\delta _is_i^z`$ (10)
$`s_i^y`$ $``$ $`s_i^y`$ (11)
$`s_i^z`$ $``$ $`\mathrm{sin}\delta _is_i^x+\mathrm{cos}\delta _is_i^z.`$ (12)
A similar rotation about the $`y`$-axis by an angle $`\delta _j`$ is performed for the spin $`\text{s}_j`$. Thus we get for the transformation of the scalar product of the two spins, $`\text{s}_i\text{s}_j(\text{s}_i\text{s}_j)_\phi `$, where

$`(\text{s}_i\text{s}_j)_\phi `$ $``$ $`\mathrm{sin}\phi [s_i^xs_j^zs_i^zs_j^x]+\mathrm{cos}\phi [s_i^xs_j^x+s_i^zs_j^z]+s_i^ys_j^y`$ (13)
$`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{sin}\phi [s_i^+s_j^zs_i^zs_j^++s_i^{}s_j^zs_i^zs_j^{}]+\mathrm{cos}\phi s_i^zs_j^z`$ (16)
$`+{\displaystyle \frac{1}{4}}(\mathrm{cos}\phi +1)[s_i^+s_j^{}+s_i^{}s_j^+]`$
$`+{\displaystyle \frac{1}{4}}(\mathrm{cos}\phi 1)[s_i^+s_j^++s_i^{}s_j^{}].`$
The angle $`\phi \delta _j\delta _i`$ is the angle between the two spins, and $`s^\pm s^x\pm \mathrm{i}s^y`$ are the spin-raising and spin-lowering operators. Note that this product of two spins after the rotation depends not only on the angle between them, but also on the sign of this angle. In case of the Néel model state ($`\mathrm{\Phi }=0`$), the angle between any neighbouring spins is $`\pi `$, and hence Eq. (13) becomes $`\text{s}_i\text{s}_j\frac{1}{2}[s_i^+s_j^++s_i^{}s_j^{}]s_i^zs_j^z`$.

Using the helical state of Eq. (4) with the characteristic angle $`\mathrm{\Phi }`$, the Hamiltonian of Eq. (1) is now rewritten in the local coordinates as,

$$H=\underset{iA}{}\underset{p}{}(1+\delta _{p,p_J^{}}(J^{}1))(\text{s}_i\text{s}_{i+p})_{\phi _p},$$
(17)
where the angles between neighbouring spins are $`\phi _{\pm \widehat{y}}=\pi +\mathrm{\Phi }`$, $`\phi _{\widehat{x}}=\pi \mathrm{\Phi }`$ and $`\phi _{\widehat{x}}=\pi +3\mathrm{\Phi }`$. While the general helical state (see Fig. 1) does not have translational symmetry in the $`x`$-direction, the transformed Hamiltonian of Eq. (17) does have this symmetry since it depends only on the angles between neighbouring spins.

Having defined an appropriate model state $`|\mathrm{\Phi }`$ with creation operators $`C_I^+`$, the CCM parameterizations of the ket and bra ground states are given by

$$|\mathrm{\Psi }=e^S|\mathrm{\Phi },S=\underset{I0}{}𝒮_IC_I^+,$$
(18)
$$\stackrel{~}{\mathrm{\Psi }}|=\mathrm{\Phi }|\stackrel{~}{S}e^S,\stackrel{~}{S}=1+\underset{I0}{}\stackrel{~}{𝒮}_IC_I.$$
(19)
The correlation operator $`S`$ is expressed in terms of the creation operators $`C_I^+`$ and the ket-state correlation coefficients $`𝒮_I`$. We can now write the ground-state energy as,

$$E=\mathrm{\Phi }|e^SHe^S|\mathrm{\Phi }.$$
(20)
To describe the magnetic order of the system, we use a simple order parameter which is expressed in terms of the local, rotated spin axes, and which is given by

$$M\stackrel{~}{\mathrm{\Psi }}|s_i^z|\mathrm{\Psi },$$
(21)
such that the order parameter represents the on-site magnetization. Note that $`M`$ is the usual sublattice magnetization for the case of the Néel state as the CCM model state.

To find the ket-state and bra-state correlation coefficients we have to require that the expectation value $`\overline{H}=\stackrel{~}{\mathrm{\Psi }}|H|\mathrm{\Psi }`$ is a minimum with respect to the bra-state and ket-state correlation coefficients. This formalism is exact if we include all possible multispin configurations in the correlation operators $`S`$ and $`\stackrel{~}{S}`$, which is usually impossible in a practical situation. We use the LSUB$`n`$ approximation scheme to truncate the expansion of $`S`$ and $`\stackrel{~}{S}`$ in the Eqs. (18) and (19).

Using the lattice symmetries, we have now to find all different possible configurations with respect to the point and space group symmetries of both the lattice and Hamiltonian with up to $`n`$ spins spanning a range of no more than $`n`$ adjacent lattice sites (LSUB$`n`$ approximation) and these are referred to as the fundamental configurations.

The Hamiltonian of Eq. (1) has four lattice point-group symmetries namely two rotational operations ($`0^{}`$, $`180^{}`$) and two reflections (along the $`x`$\- and $`y`$-axes), defined by:

$`xx,yy;`$ $`x(x+1),yy,`$ (22)
$`xx,yy;`$ $`x(x+1),yy.`$ (23)
The rotation of $`180^{}`$ and the reflection along the $`y`$-axis are connected by a shift of $`\widehat{x}=(1,0)`$. The translational operator $`T`$ is defined by

$$T=(n+m)\widehat{x}+(mn)\widehat{y},n\text{}m\text{ integral},$$
(24)
such that translational symmetry is preserved.

The Néel model state also contains these symmetries, and so for this model state we can directly apply all these symmetries in finding the fundamental configurations. On the other hand the general helical model state ($`\mathrm{\Phi }0`$) has only two of the above four lattice point-group symmetries, namely $`xx,yy`$ and $`xx,yy`$, and so this reduced symmetry yields a larger number of fundamental configurations.

In the case of the Néel model state ($`\mathrm{\Phi }=0`$), the number of fundamental configurations can further be reduced by explicitly conserving the total uniform magnetization $`s_T^z_ks_k^z`$ (the sum on $`k`$ runs over all lattice sites) because the ground state is known to lie in the $`s_T^z=0`$ subspace. This means that we exclude configurations with an odd number of spins, and therefore we do not use LSUB3, LSUB5, etc. approximations. The helical state is not an eigenstate of $`s_T^z`$ and we cannot apply this property when using the helical model state. The fundamental configurations can now be calculated computationally, and the resulting numbers of LSUB$`n`$ configurations for $`n8`$ are given in Table I.

The Ket-state and bra-state equations are calculated computationally. For the Néel model state, we are able to carry out the CCM up to the LSUB8 level (where we need to solve 4986 coupled equations), whereas for the helical state we could do this only up to the LSUB6 level (where we need to solve 1638 coupled equations).

### B The excited state formalism

We use the excited-state formalism of Emrich to approximate the excited-state wavefunctions. We apply an excitation operator $`X^e`$ linearly to the ket state wavefunction (18), such that

$$|\mathrm{\Psi }_e=X^ee^S|\mathrm{\Phi };X^e=\underset{I0}{}𝒳_I^eC_I^+.$$
(25)
Using the Schrödinger equation, $`H|\mathrm{\Psi }_e=E_e|\mathrm{\Psi }_e`$, we find that

$$ϵ_eX^e|\mathrm{\Phi }=e^S[H,X^e]_{}e^S|\mathrm{\Phi },$$
(26)
where $`ϵ_e`$ ($`E_eE`$) is the difference between the excited-state energy ($`E_e`$) and the ground-state energy ($`E`$). Applying $`\mathrm{\Phi }|C_I`$ to Eq. (26) we find that,

$$ϵ_e𝒳_I^e=\mathrm{\Phi }|C_Ie^S[H,X^e]_{}e^S|\mathrm{\Phi },$$
(27)
which is an eigenvalue equation with eigenvalues $`ϵ_e`$ and corresponding eigenvectors $`𝒳_I^e`$.

As for the ground state, we must use an approximation scheme for $`X^e`$ in Eq. (25). Although it is not necessary to use the same approximation for the excited state as for the ground state, we in fact do so to keep the CCM calculations as systematic and self-consistent as possible. We define the fundamental configurations for LSUB$`n`$ (for the Néel state) as previously, but we now restrict the choice of configurations to contain only those which produce a change of $`s_T^z`$ of $`\pm 1`$ with respect to the model state. Since we are only interested in the lowest-lying excitations, the restriction to these single-magnon spin-wave-like excitations is the correct choice. The number of fundamental excited-state configurations for LSUB$`n`$ is given in Table I.

To calculate the terms of the right hand side of Eq. (27) we use the same computational algorithm as for the calculation of the ground-state ket equation. The terms contain the ground ket-state correlation coefficients $`𝒮_I`$, so once these coefficients have been determined the eigenvalue equation (27) can be solved (numerically). We furthermore choose the lowest energy eigenvalue of Eq. (27) in order to calculate the excitation energy gap, $`\mathrm{\Delta }`$. We note that the eigenvalues of Eq. (27) are not guaranteed to be real, since as a generalized eigenvalue equation it is not symmetric. However, over the entire regime of interest, the values of $`\mathrm{\Delta }`$ so obtained are found to be real. We have performed these calculations for the excited state up to the LSUB6 level of approximation.

### C Extrapolation of the CCM-LSUB$`n`$ results

Although no scaling theory for results of LSUB$`n`$ approximations has yet been proven, there are empirical indications of scaling laws for the energy, the magnetization, and the excited-state energy gap for various spin models.

These scaling laws can be justified by the observations that they fit the results well (i.e., with low mean-square deviation), and that the extrapolated results are in good agreement with results of other methods (e.g., Green function Monte Carlo or series expansion for the 2D $`XXZ`$ model) or with exact results (e.g., 1D $`XY`$ model). In accordance with those previous results we use the following scaling laws: for the ground-state energy,

$$E=a_0+a_1(1/n^2)+a_2(1/n^2)^2;$$
(28)
for the ground-state magnetization,

$$M=b_0+b_1(1/n)+b_2(1/n)^2;$$
(29)
and for the gap of the lowest-lying excitations,

$$\mathrm{\Delta }=c_0+c_1(1/n)+c_2(1/n)^2;$$
(30)
where $`n`$ is the LSUB$`n`$ approximation level.

### D Choice of the CCM model state

As stated previously, we use the helical state of Eq. (4) with the characteristic angle $`\mathrm{\Phi }`$, illustrated in Fig. 1, as the model state for the CCM. We must therefore make a selection of an appropriate value for $`\mathrm{\Phi }`$. A possible choice would be the classical ground state of the Hamiltonian of (1) (i.e., $`\mathrm{\Phi }=\mathrm{\Phi }_{\mathrm{cl}}`$ as given by Eq. (3)).

Another possibility is to perform a CCM-LSUB$`n`$ approximation calculation and then to minimize the corresponding LSUB$`n`$ approximation to the energy with respect to $`\mathrm{\Phi }`$,

$$E_{\mathrm{LSUB}n}(\mathrm{\Phi })\mathrm{min}\mathrm{\Phi }=\mathrm{\Phi }_{\mathrm{LSUB}n}.$$
(31)
The results for $`\mathrm{\Phi }_{\mathrm{LSUB}n}`$ will be given later (Fig. 7). However, we note now that although the CCM does not yield a strict upper bound for the ground-state energy, using $`\mathrm{\Phi }=\mathrm{\Phi }_{\mathrm{LSUB}n}`$ (i.e., using the CCM with a variational parameter) has been found to be a reasonable assumption. There are several additional arguments to suggest that $`\mathrm{\Phi }=\mathrm{\Phi }_{\mathrm{LSUB}n}`$ is indeed a better choice than $`\mathrm{\Phi }=\mathrm{\Phi }_{\mathrm{cl}}`$, as indicated below.

In the first place we note that we cannot find solutions for the LSUB6 equations using $`\mathrm{\Phi }=\mathrm{\Phi }_{\mathrm{cl}}(J^{})`$ in the region $`0.7J^{}0.47`$ insofar as the Newton method used to solve these equations does not converge in that region. This is a clear indication that this model state is not a good one. By contrast, such behaviour is not found for $`\mathrm{\Phi }=\mathrm{\Phi }_{\mathrm{LSUB}n}`$.

Secondly, it is generally known that quantum fluctuations tend to prefer collinear order (e.g., Néel order). We will indeed find (and see Fig. 7 below) that the Néel ordering ($`\mathrm{\Phi }=0`$) seems to survive for some $`J^{}<1/3`$, in which region it has already broken down in the classical case. This is also in agreement with results of exact diagonalizations for our model (and see Fig. 5 below).

Thirdly we find better agreement of the CCM results for the energy compared to exact diagonalization results by using the helical state as the model state with the value $`\mathrm{\Phi }=\mathrm{\Phi }_{\mathrm{LSUB}n}`$ rather than with the classical value $`\mathrm{\Phi }=\mathrm{\Phi }_{\mathrm{cl}}`$. We find that CCM results for the ground-state energy usually agree well with the corresponding ED results (and with results of other methods), provided that a good CCM model state is chosen.

We therefore use the helical state with $`\mathrm{\Phi }=\mathrm{\Phi }_{\mathrm{LSUB}n}`$ as the CCM model state throughout this paper. Note that for $`J^{}1/3`$ this model state is identical to the classical ground state of Eq. (1) but that for $`J^{}<1/3`$ it is not.

## IV Results

Using the CCM scheme described above, we calculate the approximate ground state and the low-lying excitations of the Hamiltonian of Eq. (1). For comparison we also exactly diagonalize finite sized lattices of square shape. We use periodic boundary conditions with $`N=16,18,20,26`$ and $`32`$ spins, and we extrapolate to the infinite system using standard finite-size scaling laws. We present results for the ground-state energy, the order parameter and the excitation gap. We examine the formation of local singlets (for $`J^{}>1`$), the effects of frustration (for $`J^{}<0`$), and the special case of the honeycomb lattice ($`J^{}=0`$).

$`J^{}>1`$: Formation of local singlets. Using the CCM we obtain clear indications of a second-order phase transition to a disordered dimer-like phase at a certain critical value of $`J^{}`$, namely $`J_s^{}`$. For $`J^{}>J_s^{}`$, the Néel-like long range order melts, (i.e., the sublattice magnetization, $`M`$, given by Eq. (21) becomes zero). Our estimate for $`J_s^{}`$ using the four extrapolated LSUB$`n`$ results for $`M`$ with $`n=2,4,6,8`$ (see Fig. 2) is $`J_s^{}3.41`$. However, using only the three CCM LSUB$`n`$ approximations with $`n=4,6,8`$ for the extrapolation, we obtain a value $`J_s^{}3.16`$, which indicates that the true value could be even somewhat smaller. This is in agreement with our corresponding result using exact diagonalizations of small systems. By using the extrapolation scheme of Ref. , we find a critical value $`J_s^{}2.45`$ for the magnetization. Note however that better accuracy requires larger systems because of the exact diagonalization (ED) extrapolation ansatz for $`M`$ (i.e., $`M=M_{\mathrm{}}+\mathrm{const}\times N^{1/2}`$). Therefore, we cannot consider the ED results for the magnetization (and see Fig. 2) as quantitatively correct. Our two results for the critical $`J_s^{}`$ also agree with the estimate $`J_s^{}2.56`$ from series expansion, and even the result $`J_s^{}=3`$ from the simple variational ansatz of Eq. (6) agrees surprisingly well with these values. By contrast, the second-order renormalized spin wave theory (RSWT) gives the larger result $`J_s^{}5.0`$, indicating that the standard spin-wave approach is insufficient to describe this type of transition.

Another indication of a dimerized phase is the appearance of a gap $`\mathrm{\Delta }`$ between the ground state and the lowest-lying excited state. We clearly expect a spectrum with gapless Goldstone modes if the ground state is Néel long-range ordered, whereas for a disordered singlet ground-state the formation of triplet excitations may cost a finite amount of energy. This behaviour is reflected by our results using both CCM and exact diagonalization (see Fig. 3), which agree well with each other. For $`J^{}J_s^{}`$, there is a gap proportional to $`J^{}`$, corresponding to the dimer-like nature of the ground state. The gap obviously opens in the range $`2.5J_s^{}3.0`$ in both the ED and CCM calculations. This is in good agreement with the corresponding estimates for the critical point using the order parameter. Note, that the standard linear spin wave theory fails in calculating the gap (i.e., gives gapless modes for all values of $`J^{}`$ with $`J^{}>0`$).

By comparing the results for the ground-state energy, we find excellent agreement between the CCM results and the results from exact diagonalization (see Fig. 4) for $`J^{}>0`$. By contrast, spin-wave theory (SWT) calculations show a significant deviation from these results for larger $`J^{}1`$. These spin-wave results are obviously poor since the simple upper bound for the energy given by Eq. (7) (e.g., $`E_0=1.5`$ for $`J^{}=4`$) is smaller than the corresponding SWT results (e.g., $`E_0=1.42`$ from 2nd-order RSWT). By contrast, CCM and ED results (both are about $`E_0=1.54`$ for $`J^{}=4`$) are slightly smaller than the variational result. While both CCM and SWT calculations have the Néel state as starting point, we find the CCM is much better able than SWT to describe the transition to the rotationally invariant disordered state and to the completely dimerized state (represented by the variational function of Eq. (6) with $`t=1`$). Note that even the simplest CCM approximation (LSUB2) gives the correct asymptotic result for the energy (i.e., Eq. (7)) for very large values of $`J^{}`$, whereas SWT does not.

For the case of the pure square-lattice Heisenberg antiferromagnet (i.e., $`J^{}=1`$), we reproduce the CCM results of Refs. , which have already been demonstrated to agree well with those from other methods.

$`J^{}=0`$: honeycomb lattice. For the special case of $`J^{}=0`$ (which is equivalent to the honeycomb lattice), we find that the CCM and the ED results are in good agreement (see Table II). However, the magnetization $`M`$ for ED is found to be smaller then the CCM result. We note, however, that the CCM result for $`M`$ at this point agrees with the result of high-order SWT ($`M=0.28`$) as well as with the result of series expansion ($`M=0.26`$), although it does not agree so well with the result of Monte Carlo calculations ($`M=0.22`$). We note too that our CCM results here agree perfectly with previous lower-order CCM calculations.

$`J^{}<0`$: Frustration. For $`J^{}2`$, we find that the extrapolated ED results for the energy lie appreciably above the CCM (and SWT) results (see Fig. 4). This is because the energies for the small lattices considered do not fit well to the finite-size scaling law ($`E/N=E_{\mathrm{}}/N+\mathrm{const}\times N^{3/2}`$) in this region. The finite-size effects for systems with an incommensurate helical structure are found to be larger than for systems with, for example, Néel order or with dimerized spin pairs. However, we find that our best ED result (with 32 spins) shows only very small deviations from the CCM result, even in the frustrated region.

While classically we have a second-order phase transition (from Néel order to helical order) at $`J_c^{}=1/3`$, using the CCM we find indications for a shift of this critical point to a value $`J_c^{}1.35`$ (see Fig. 7) in the quantum case. The ED data of the structure factors (see Fig. 5) also show a shift of the transition to stronger ferromagnetic $`J^{}`$ bonds. Both correspond to the general picture that quantum fluctuations prefer a collinear ordering (such as Néel order). Hence this ordered state can survive for the quantum case into a region where classically it is already unstable. The Néel model state ($`\mathrm{\Phi }=0`$) gives the minimum ground-state energy for all values $`J^{}>J_c^{}`$, where $`J_c^{}`$ is also dependent on the level of LSUB$`n`$ approximation level. For $`J^{}<J_c^{}`$ another minimum in the energy for $`\mathrm{\Phi }0`$ is found to lie lower than the minimum at $`\mathrm{\Phi }=0`$ (see Fig. 6). The state for $`\mathrm{\Phi }0`$ is believed to be a quantum analogue of the classical ground state in two dimensions. Furthermore, the crossover from one minimum solution to the other is not smooth but is abrupt at this point (see Fig. 6 and Fig. 7). This behaviour is assumed here to be an indication of a phase transition. Furthermore, it might also indicate that this is a first-order phase transition and, consequently, that due to quantum fluctuations the nature of this phase transition is changed from the classical second-order type to a first-order type.

The behaviour of the order parameter (see Fig. 2) in the region around $`J_c^{}1.35`$ where we expect the above-mentioned quantum phase transition is also quite marked. We cannot extrapolate the LSUB$`n`$ results directly, because the phase transition points shift with the order of the LSUB$`n`$ approximations (Fig. 7). We thus find a large statistical deviation of the extrapolated results in the region $`1.4J^{}1.0`$. Hence, we use the minima for $`M`$ in that region to extrapolate an estimation of the order parameter. We find that minimum ($`M0.05`$) to be at a value of $`J^{}1.2`$. The extrapolated ED results do not agree very well with this result, since these give $`M`$ to be zero at $`J^{}0.8`$. However, these results are also very poor in that region because of the strong influence of the boundary conditions and large statistical errors. As a result of these difficulties we are not able to decide whether or not quantum fluctuations and frustration are able to form a disordered quantum spin liquid phase (i.e., with $`M=0`$) between the Néel state and the helical state for some finite frustrating $`J^{}<0`$ regime. However, the CCM results suggest that there is either no quantum spin liquid phase or that, if it does exist, it does so only in a very small region.

## V Summary

Using the CCM we have studied the influence of quantum spin fluctuations on both the ground-state phase diagram and the excited states of a spin-half square-lattice Heisenberg antiferromagnet with two kinds of nearest-neighbour exchange bonds. The phase diagram is found to contain a quantum helical phase, a Néel ordered phase, and a finite-gap quantum disordered phase. While we have clearly a second-order transition from the Néel phase to the finite-gap quantum disordered phase, we also found indications of a quantum-induced first-order transition from the Néel phase to the helical phase, for which classically we have a second-order transition. While our CCM results were in general in good agreement with the ED data, we found the CCM particularly good at describing the dimerized phase. By contrast, spin wave theory fails in that region due to enhanced longitudinal spin fluctuations. Accurate high-order CCM results for the antiferromagnet on the honeycomb lattice were also presented.

###### Acknowledgements.

We would like to thank N.B. Ivanov for his stimulating discussions. This work has also been supported in part by the Deutsche Forschungsgemeinschaft (GRK 14, Graduiertenkolleg on “Classification of Phase Transitions in Crystalline Materials”), and by a research grant (GR/M45429) from the Engineering and Physical Sciences Research Council (EPSRC) of Great Britain.