# Various series expansions for a Heisenberg antiferromagnet model for SrCu2(BO3)2

## I INTRODUCTION

The discovery of high-temperature superconductivity has stimulated an enormous amount of activity in the study of two-dimensional antiferromagnetism, which may be connected with the superconductivity phenomenon. The pseudo spin-gap behaviour observed in the high T<sub>c</sub> cuprates has also stimulated intense interest in systems with spin gaps. Several new spin gap systems have been found experimentally. Among them some of the compounds which have two dimensional character include the coupled spin ladder systems, SrCu<sub>2</sub>O<sub>3</sub>, CaV<sub>2</sub>O<sub>5</sub>, (VO<sub>2</sub>)P<sub>2</sub>O<sub>7</sub>, Cu<sub>2</sub>(C<sub>5</sub>H<sub>12</sub>N<sub>2</sub>)<sub>2</sub>Cl<sub>4</sub>, and the plaquette RVB system, CaV<sub>4</sub>O<sub>9</sub>.

Recently a new two dimensional spin gap system SrCu<sub>2</sub>(BO<sub>3</sub>)<sub>2</sub> has been found by Kageyama et al.. It has a spin-singlet ground state with a finite spin gap $`20`$ K. They also found that the peak of the susceptibility is much suppressed compared with standard dimer models, and observed two plateaus in the magnetization at 1/4 and 1/8 of the full moment.

Miyahara and Ueda showed that these observations could be understood on the basis of a simple two-dimensional Heisenberg antiferromagnet model with nearest-neighbor and next-nearest-neighbor couplings. The copper ions in the SrCu<sub>2</sub>(BO<sub>3</sub>)<sub>2</sub> compound are all located at crystallographically equivalent sites, forming a distinctive pattern. Miyahara and Ueda show that, remarkably enough, a singlet dimer state forms an exact eigenstate of this Hamiltonian at all couplings, and is the ground-state in a region where the nearest-neighbour coupling dominates. There is a first order phase transition to a Néel ordered state at $`J^{}/J=0.7\pm 0.01`$. They find that the SrCu<sub>2</sub>(BO<sub>3</sub>)<sub>2</sub> system lies close to this transition, which explains the unusual temperature dependence of the magnetization. The plateaus observed in the magnetization curve can also be understood, on the basis that the triplet excitations from the ground state are almost localized. Their conclusions were reached on the basis of exact diagonalization calculations for lattices of up to 20 sites, and low order dimer expansions in $`(J^{}/J)`$, and high temperature expansions in $`J/T`$ and $`J^{}/T`$.

In this paper, we reinforce these conclusions by carrying out a more extensive series study of the model. These include high-temperature expansions, Ising expansions at zero temperature, and dimer expansions both at zero and finite temperature. The model is presented in Section II, and the framework for the various series expansions is outlined. The results are presented in Section III. It is shown that the temperature dependence of the susceptibility can be fitted accurately and in detail by this method.

## II Series Expansions

The magnetic properties of SrCu<sub>2</sub>(BO<sub>3</sub>)<sub>2</sub> may be described by the two-dimensional spin $`S=1/2`$ Heisenberg antiferromagnetic model with nearest-neighbor (n.n.) and next-nearest-neighbor (n.n.n.) interactions:

$$H=J\underset{\mathrm{n}.\mathrm{n}.}{}𝐒_i𝐒_j+J^{}\underset{\mathrm{n}.\mathrm{n}.\mathrm{n}.}{}𝐒_i𝐒_j.$$
(1)   
The system is illustrated in Fig. 1(a). We denote the ratio of the couplings as $`x`$, that is, $`xJ^{}/J`$. In the present paper, we study only the case of antiferromagnetic coupling, where both $`J`$ and $`J^{}`$ are positive. In the large $`J^{}/J`$ limit, the model is topologically equivalent to the two dimensional nearest-neighbor square lattice Heisenberg model. In the small $`J^{}/J`$ limit, every pair of spins along nearest-neighbor bonds interact only weakly with each other, and the dominant configuration in the ground state is the product state with the spins along nearest-neighbor bonds forming a spin singlet. In fact, it has been proved that this state is an eigenstate of the system for any $`x`$, and it is the ground state for small enough values of $`x`$. We can also prove that this system still has this perfect dimer state as an eigenstate even if we include the coupling $`J_2`$ denoted by dotted lines shown in Fig. 1(b).

We have studied this system by using various linked-cluster expansion methods including dimer expansions at both zero-temperature and finite-temperature, Ising expansions at zero-temperature, and high temperature expansions. The linked-cluster expansion method has been previously reviewed in several articles, and will not be repeated here. Here we only summarize the expansion methods used, and the results derived from them are presented in the next Section.

### A Dimer Expansions at $`T=0`$

At temperature $`T=0`$, we can construct an expansion in $`x`$ by taking the first term of $`H`$ as the unperturbed Hamiltonian and the second term in $`H`$ as a perturbation. That is, the Hamiltonian of Eq.(1) can be rewritten as

$$H=H_0+xV,$$
(2)   
where

$`H_0`$ $`=`$ $`{\displaystyle \underset{\mathrm{n}.\mathrm{n}.}{}}𝐒_i𝐒_j,`$ (3)   
$`V`$ $`=`$ $`{\displaystyle \underset{\mathrm{n}.\mathrm{n}.\mathrm{n}.}{}}𝐒_i𝐒_j,`$ (5)   
and where we have set $`J=1`$ for convenience. The unperturbed ground state is then a product state of nearest neighbor singlet dimers and the perturbation couples these among themselves and with the pair triplet states. As mentioned above, the unperturbed ground state is also an eigenstate of the full Hamiltonian, but is not the true ground state for $`x>x_c`$.

Dimer expansions can be developed for all ground state properties as well as for the triplet excitation spectrum. Here, because of the trivial nature of the ground state in the dimer phase, we concentrate on the lowest triplet excitations. We have calculated the dispersion relation $`\mathrm{\Delta }(k_x,k_y)`$ to order $`x^{15}`$, extending the calculation of Miyahara and Ueda by 11 terms. This calculation involves 11586 linked clusters up to 8 sites. The resulting series coefficients are given in Table I.

### B Ising Expansions at $`T=0`$

In the limit that $`J^{}J`$, the model is topologically equivalent to the two dimensional square lattice Heisenberg model, so we expect that the system has Neél order: an order in which every pair of spins along a horizontal nearest-neighbor bond (denoted as A) has spin up, while every pair of spins along a vertical nearest-neighbor bond (denoted as B) has spin down.

To construct a $`T=0`$ expansion about the Ising limit for this system, one has to introduce an anisotropy parameter $`\lambda `$, and write the Hamiltonian for the Heisenberg-Ising model as

$$H=H_0^{}+\lambda V^{},$$
(6)   
where

$`H_0^{}`$ $`=`$ $`{\displaystyle \underset{\mathrm{n}.\mathrm{n}.}{}}S_i^zS_j^z+x{\displaystyle \underset{\mathrm{n}.\mathrm{n}.\mathrm{n}.}{}}S_i^zS_j^z+t[{\displaystyle \underset{iB}{}}S_i^z{\displaystyle \underset{iA}{}}S_i^z],`$ (7)   
$`V^{}`$ $`=`$ $`{\displaystyle \underset{\mathrm{n}.\mathrm{n}.}{}}(S_i^xS_j^x+S_i^yS_j^y)+x{\displaystyle \underset{\mathrm{n}.\mathrm{n}.\mathrm{n}.}{}}(S_i^xS_j^x+S_i^yS_j^y)t[{\displaystyle \underset{iB}{}}S_i^z{\displaystyle \underset{iA}{}}S_i^z],`$ (8)   
The last term in both $`H_0^{}`$ and $`V^{}`$ is a local staggered field term, which can be included to improve convergence. The limits $`\lambda =0`$ and $`\lambda =1`$ correspond to the Ising model and the isotropic Heisenberg model, respectively. The operator $`H_0^{}`$ is taken as the unperturbed Hamiltonian, with the unperturbed ground state being the usual Néel state. The operator $`V^{}`$ is treated as a perturbation. It flips a pair of spins on neighboring sites.

Ising series have been calculated for the ground state energy per site, $`E_0/N`$ and the staggered magnetization $`M`$, for several ratios of couplings $`x`$ and (simultaneously) for several values of $`t`$ up to order $`\lambda ^9`$. The resulting series for $`x=0.7,0.8,1,2`$ and $`t=0`$ are listed in Table II; the series for other value of $`x`$ are available upon request.

At the next stage of the analysis, we try to extrapolate the series to the isotropic point ($`\lambda =1`$) for those values of the exchange coupling parameters which lie within the Néel-ordered phase at $`\lambda =1`$. For this purpose, we first transform the series to a new variable

$$\delta =1(1\lambda )^{1/2},$$
(9)   
to remove the singularity at $`\lambda =1`$ predicted by the spin-wave theory. This was first proposed by Huse and was also used in our earlier work on the square lattice case. We then use both integrated first-order inhomogeneous differential approximants and Padé approximants to extrapolate the series to the isotropic point $`\delta =1`$ ($`\lambda =1`$). The results of the Ising expansions will be presented in the next section.

### C High-temperature series expansions

We now turn to the finite-$`T`$ thermodynamic properties. We have developed high-temperature series expansions $`\beta =1/(k_BT)`$ for the uniform magnetic susceptibility $`\chi (T)`$ and the specific heat $`C(T)`$, for the system with arbitrary $`x`$,

$`T\chi (T)`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}{\displaystyle \underset{j}{}}{\displaystyle \frac{\mathrm{Tr}S_i^zS_j^ze^{\beta H}}{\mathrm{Tr}e^{\beta H}}},`$ (10)   
$`C(T)`$ $`=`$ $`{\displaystyle \frac{U}{T}},`$ (12)   
where $`N`$ is the number of sites, and the internal energy $`U`$ is defined by

$$U=\frac{\mathrm{Tr}He^{\beta H}}{\mathrm{Tr}e^{\beta H}}.$$
(13)   
The series for $`\chi (T)`$ and $`C(T)`$ were computed to order $`\beta ^7`$ for arbitrary $`x`$. The expansions have the following general form:

$$\underset{i=0}{\overset{\mathrm{}}{}}\left[\underset{j=0}{\overset{i}{}}c_{i,j}x^j\right]\beta ^i,$$
(14)   
where $`c_{i,j}`$ are the numerical coefficients. The series coefficients are not presented here, but are available upon request. Some of the coefficients can be obtained by using the dimer expansion series at finite temperature discussed in the next subsection. The high temperature series expansions were first carried out by Miyahara and Ueda up to order $`\beta ^2`$. Our results agree with these previous results, and extend the series by 5 terms.

### D Dimer expansions at finite temperature

The study of finite temperature properties via series expansions is usually done by the high temperature expansion method as mentioned above, where we expand in powers of $`\beta `$ for given ratio of exchange couplings. This method often performs poorly at low temperatures for many parameter regimes of interest. To overcome this difficulty, one can develop the dimer expansion at finite temperature, where we expand the thermodynamic quantities in power of $`x`$ for arbitrary temperatures. This method has been used by one of us previously for the Hubbard model  and more recently also by Elstner and Singh for spin models. This method has shown excellent convergence for a wide range of coupling constants at all temperatures for the bilayer Heisenberg model and for alternating spin-chains and spin ladders.

To get an expansion for thermodynamic quantities, such as the susceptibility and the specific heat, in $`x`$ at arbitrary temperature for the Hamiltonian in Eq.(2), one basically needs to expand

$$e^{\beta (H_0+xV)}$$
(15)   
in powers of $`x`$. This can be obtained by using the following relation

$$e^{\beta (H_0+xV)}=e^{\beta H_0}\underset{n}{}(x)^nI_n,$$
(16)   
where $`I_n`$ are the integrals given by

$$I_n=_0^\beta 𝑑t_1_0^{t_1}𝑑t_2\mathrm{}_0^{t_{n1}}𝑑t_n\stackrel{~}{V}(t_1)\stackrel{~}{V}(t_2)\mathrm{}\stackrel{~}{V}(t_n),$$
(17)   
with

$$\stackrel{~}{V}(t)=e^{tH_0}Ve^{tH_0}.$$
(18)   
The integrations needed in this expansion are of the type:

$$t^ke^{lt},$$
(19)   
and can easily be integrated analytically.

The expansions for the susceptibility $`\chi `$ and the logarithm of the partition function $`\mathrm{ln}Z`$ per site take the following general form:

$`T\chi (x,\beta )`$ $`=`$ $`{\displaystyle \underset{n=0}{}}{\displaystyle \frac{f_n^{(\chi )}(\beta )}{12n!}}\left({\displaystyle \frac{x}{12}}\right)^n,`$ (20)   
$`\mathrm{ln}Z(x,\beta )`$ $`=`$ $`3\beta /4+\mathrm{ln}Z_0+{\displaystyle \underset{n=1}{}}{\displaystyle \frac{f_n^{(Z)}(\beta )}{12n!}}\left({\displaystyle \frac{x}{12}}\right)^n,`$ (22)   
where $`Z_0=1+3\mathrm{exp}(\beta )`$, and the coefficient $`f_n(\beta )`$ for both $`\chi `$ and $`\mathrm{ln}Z`$ have the following general form

$$f_n(\beta )=\underset{k=0}{\overset{n}{}}\underset{l=0}{\overset{n+1}{}}c_{k,l}^{(n)}\beta ^kZ_0^l,$$
(23)   
where $`c_{k,l}^{(n)}`$ are expansion coefficients. The results up to order $`n=6`$ are given in Table III. From this expansion, one can recover the results of the high-temperature expansion and dimer expansions if we reexpand in powers of $`\beta `$, or $`x`$ at $`T=0`$.

## III Results

Having obtained the series for the various expansions above we present in this section the results of series analysis. We use integrated first-order inhomogeneous differential approximants and Padé approximants to extrapolate the series.

### A Phase diagram

The ground state energy per site $`E_0/N`$ is shown in Figure 2. The full points in the large $`J^{}/J`$ region are obtained from the Ising expansion, and the horizontal line corresponds to the eigenenergy of the perfect dimer state (which is $`E_0=3NJ/8`$ exactly). These curves cross at a transition point $`x_c`$, which corresponds to a first order ground state phase transition, resulting from a level crossing. The numerical estimate of $`x_c=0.691(6)`$ is a more precise estimate of the result $`x_c=0.70(1)`$ of Miyahara and Ueda discussed above.

The staggered magnetization $`M`$ for those values of the exchange coupling parameters which lie within the Néel-ordered phase at $`\lambda =1`$ is shown in Fig. 3. We can see that $`M`$ decreases as we turn on $`J`$, and appears to vanish at around $`J/J^{}1.4`$, near the transition point determined from the ground state energies. The errors are too large to determine whether or not where is a finite discontinuity at the transition.

### B Triplet excitation spectrum

From the dimer expansions, one can estimate the triplet excitation spectrum for those values of the exchange coupling parameters which lie within the dimer phase. The triplet excitation spectra for $`x=0.6,0.65`$ and 0.678 are shown in Fig. 4. We can see that the minimum energy gap is at $`(k_x,k_y)=(0,0)`$ (and the equivalent point $`(\pi ,\pi )`$). The bandwidth is quite small, indicating that the triplet excitations are highly localized, but increases as $`x_c`$ is approached.

Fig. 5 shows the triplet excitation gap $`\mathrm{\Delta }=\mathrm{\Delta }(0,0)`$ as a function of $`x`$. At $`J^{}=0`$, $`\mathrm{\Delta }/J`$ is equal to 1 exactly, corresponding to a single dimer excited to a triplet state. At the first order transition point $`x_c`$, the gap is $`\mathrm{\Delta }=0.14(5)`$, where the uncertainty is largely associated with the uncertainty of $`x_c`$.

### C Thermodynamic properties

Let us now discuss the thermodynamic properties. Fig. 6 shows the internal energy per site $`E/NJ`$ vs. temperature $`k_BT/J`$ obtained from the high temperature expansion and from the dimer expansion at finite temperature for $`x=0.5`$. We can see that the direct sum to the dimer expansion series at finite temperature converges extremely well down to $`T=0`$, and it can recover the exact ground state energy $`E_0/NJ=3/8`$ at $`T=0`$. The results obtained from the integrated differential approximants to the high-temperature series expansion only converge well down to $`k_BT/J0.3`$.

The results for the specific heat at $`x=0.4`$ obtained from the high-temperature expansion and from the dimer expansion at finite temperature are shown at Fig. 7. We choose $`x=0.4`$ rather than $`x=0.5`$ because the specific heat series converges poorly for larger $`x`$.

Finally, the results for the susceptibility $`\chi `$ at $`x=0.5`$ are shown in Fig. 8. Again, we can see that the results from a direct sum of the dimer expansion series at finite temperature converges well down to $`T=0`$, (it converges less well near the peak position), while the integrated differential approximants to the high-temperature expansion series only converge to the position of the peak.

It is clear that the dimer expansion at finite temperature gives much better results than the high-temperature expansion for those values of the exchange coupling parameters which lie within the dimer phase. It is interesting to explore this also for ratios of the couplings lying within the Neél-ordered phase at $`T=0`$ ($`x>0.691`$). In Fig. 9 and Fig. 10, we present the results for internal energy and susceptibility for $`x=1.5`$, well outside the dimer phase at $`T=0`$, and one can see that dimer expansions at finite temperature still give good convergence in the high temperature region. Evidently at high temperature the system is highly disordered and both methods include contributions from all states.

### D Compound SrCu<sub>2</sub>(BO<sub>3</sub>)<sub>2</sub>

Finally we compare the experimental results for the compound SrCu<sub>2</sub>(BO<sub>3</sub>)<sub>2</sub> with the theoretical calculations to get an estimate of the exchange constants. Kageyama et al. determined the excitation gap $`\mathrm{\Delta }=25`$ K from the NMR relaxation rate. They also measure the temperature dependence of the magnetic susceptibility, and this implies a slightly smaller gap: $`\mathrm{\Delta }=(19\pm 1)`$ K. Since there is some uncertainty in the energy gap, to determine the exchange constants we begin with the magnetic susceptibility. To convert the theoretical $`\chi (T)`$ into units emu/Cu mol in experiment, we multiply our $`\chi `$ (calculated with $`J=1`$) by $`N_Ag^2\mu _B^2/k_B`$, with $`\mu _B`$ the Bohr magneton, $`k_B`$ the Boltzmann constant, $`N_A`$ Avogadro’s number, and $`g`$ the effective $`g`$-factor. There are three fitting parameters: $`J`$ and $`x`$ (or $`J^{}`$) and $`g`$. Our goal is try to find a proper parameter set $`(J,x,g)`$ which gives the minimum value of

$$P=\underset{T_i}{}[\chi ^{\mathrm{exp}}(T_i)\chi ^{\mathrm{theo}}(T_i)]^2,$$
(24)   
where $`\chi ^{\mathrm{exp}}`$, $`\chi ^{\mathrm{theo}}(T_i)`$ are the experimental and theoretical susceptibilities, respectively, and the summation is over all experimental points $`T_i`$. As the theoretical susceptibility from the 6th order dimer expansion at finite temperature converges well above temperatures corresponding to a physical temperature $`T_0100`$ K, we restrict the minimisation function to temperatures $`T>T_0`$. We then scan the values of $`P`$ around the expected parameter region for SrCu<sub>2</sub>(BO<sub>3</sub>)<sub>2</sub> to locate the minimum $`P_{\mathrm{min}}=1.14\times 10^9(\mathrm{emu}/\mathrm{CuMol})^2`$ at $`J=92.0`$ K, $`x=0.563`$ and $`g=2.104`$. This is shown as the full point in Fig. 11

Since there is some error in both the experimental and theoretical $`\chi `$, we cannot expect that the true parameters for SrCu<sub>2</sub>(BO<sub>3</sub>)<sub>2</sub> are located exactly at the minimum point of $`P`$, (in fact, the location of the minimum point depends on the chosen value of $`T_0`$), we need to consider the energy gap also. In the 3-dimensional parameter space $`(J,x,g)`$, there is a long tube-like region which has $`P5P_{\mathrm{min}}`$. The intersection with the plane $`g=2.108`$ is shown as the dashed lines in Fig. 11(a). If we fix the value of $`J`$ and vary $`x`$ and $`g`$, we obtain curves $`x(J)`$ and $`g(J)`$ which give near minimum values of $`P`$. These curves $`x(J)`$ and $`g(J)`$ and the corresponding $`P`$ are given in Fig. 11. We can see that along $`x(J)`$, $`P`$ only changes slightly by about 10%, and $`g`$ is about $`2.11`$ almost independent of $`J`$. To further determine the parameters for SrCu<sub>2</sub>(BO<sub>3</sub>)<sub>2</sub>, we need to consider the energy gap (which relates to the low-temperature behavior of $`\chi (T)`$). From Fig. 5, we can get two curves in $`(J,x)`$ space where the energy gap is 19 K and 25 K, respectively, as shown in Fig. 11. We see that they cross with the optimal $`x(J)`$ curve obtained above at

$$J=82.0\mathrm{K},x=0.678,g=2.108$$

if the energy gap is 19 K, or at

$$J=83.2\mathrm{K},x=0.664,g=2.108$$

if the energy gap is 25 K. We believe this should be the best fit to both $`\chi (T)`$ and the energy gap.

If we take the energy gap to be 19 K, the comparison of the experimental data with theoretical calculations is shown in Fig. 12 and Fig. 13. We can see the fit is extremely good in the high temperature region. We note that the dimer expansion sums tend to show sharp peaks, which are not visible in either the experimental data or the high-temperature expansions. We suspect that these are an artefact, as the results will become more sensitive to the truncation of the dimer expansion at finite order near the transition point in $`x`$. The ratio of couplings $`x`$ obtained here is similar to that obtained by Miyahara and Ueda, but they found $`J=100`$ K which is nearly 20% larger than our estimate. They obtained their results by considering the paramagnetic Curie-Weiss constant $`\theta `$. The Curie-Weiss constant given by the high-temperature expansions is $`\theta =(J+4J^{})/4`$. From experimental data, the susceptibility at high temperature can be fitted with a Curie-Weiss constant $`\theta =92.5`$ K and an effective $`g`$-factor $`g=2.14`$. The result curve for $`(J+4J^{})/4=92.5`$ K is shown as a dotted line in Fig. 11(a). The discrepancy between the two estimates is perhaps due to uncertainty in extracting the Curie-Weiss constant from experimental data at finite temperature. Our estimates imply the Curie-Weiss constant to be $`\theta =76.1`$ K.

## IV Conclusions

Several different series expansions have been calculated to high order for this model: high-temperature expansions, an Ising expansion at zero temperature, and dimer expansion at both zero and finite temperature. The first-order transition from the dimer phase to a Néel-ordered phase has been found to occur at $`J^{}/J=0.691(6)`$, in good agreement with the original estimate 0.70(1) of Miyahara and Ueda. The ground-state energy shows a sharp and distinct break in slope at that point, indicative of a first-order transition, and the triplet spin gap undergoes a small but definite discontinuity: the triplet spin gap is $`\mathrm{\Delta }=0.14(5)`$ at the transition point. A discontinuity in the Néel phase magnetization is less certain, but is not excluded by our results.

The model has been fitted to the experimental susceptibility data for SrCu<sub>2</sub>(BO<sub>3</sub>)<sub>2</sub>, with parameters $`g=2.108`$, $`J=82.0`$ K and $`x=0.678`$ if the energy gap is 19 K or $`g=2.108`$, $`J=83.2`$ K and $`x=0.664`$ if the energy gap is 25 K. A detailed and accurate fit can be obtained at high temperature, and a reasonable fit at low temperatures is also obtained. The fitted value of $`x=(J^{}/J)`$ is about 0.67, which is indeed very close to the transition point 0.691.

We should summarize what is new in this work. The Ising expansion has been developed for this model for the first time, and this allows a more accurate determination of the phase transition point. We have also estimated the spontaneous magnetization in the Neél phase. The finite temperature dimer expansion has not been computed previously. This provides the most reliable theoretical estimate of the high temperature susceptibility and is used to estimate the experimental parameters of the real system. Apart from this we have also substantially increased the length of the dimer series at $`T=0`$ and the conventional high temperature expansion. The former has allowed us to present the first calculation of the full triplet excitation spectrum for this system.

###### Acknowledgements.

This work forms part of a research project supported by a grant from the Australian Research Council. The computation has been performed on Silicon Graphics Power Challenge and Convex machines, and we thank the New South Wales Centre for Parallel Computing for facilities and assistance with the calculations. We would like to thank Prof. Ueda and Dr. Kageyama for providing the experimental data on SrCu<sub>2</sub>(BO<sub>3</sub>)<sub>2</sub>.