# Ab Initio Calculation of Spin Gap Behavior in CaV4O9

## Abstract

Second neighbor dominated exchange coupling in CaV<sub>4</sub>O<sub>9</sub> has been obtained from ab initio density functional (DF) calculations. A DF-based self-consistent atomic deformation model reveals that the nearest neighbor coupling is small due to strong cancellation among the various superexchange processes. Exact diagonalization of the predicted Heisenberg model yields spin-gap behavior in good agreement with experiment. The model is refined by fitting to the experimental susceptibility. The resulting model agrees very well with the experimental susceptibility and triplet dispersion.

CaV<sub>4</sub>O<sub>9</sub> was the first two-dimensional system observed to enter a low-temperature quantum-disordered phase with a spin gap $`\mathrm{\Delta }110`$K. The gap was first apparent in its susceptibility, which vanishes at low temperatures as $`\chi (T0)\mathrm{exp}(\mathrm{\Delta }/kT)`$ , and was observed directly in the dispersion of triplet spin excitations ($`\mathrm{\Omega }_Q`$) measured by neutron scattering. This unexpected behavior has stimulated considerable theoretical study of the exchange couplings between S=$`\frac{1}{2}`$ spins on the V lattice using Heisenberg models .

CaV<sub>4</sub>O<sub>9</sub> is a layered compound—the interlayer distance is sufficiently large to make interlayer V-V coupling negligible. Within a layer, the V atoms form a $`\frac{1}{5}`$-depleted square lattice shown as the circles in Fig. 1 . The lattice was originally viewed as an array of square “plaquettes” of V ions (e.g., 1-2-3-4 in Fig. 1) tending toward singlet formation since isolated plaquettes have a singlet ground state. Examination of the structure however suggests intra- and inter-plaquette nearest neighbor V-V coupling should be similar, so the limit of isolated plaquettes is not realistic.

Self-consistent electronic structure work identified the V<sup>4+</sup> spin orbital as $`d_{xy}`$, which implied that it was a larger square of V ions, the “metaplaquette,” where singlet formation arises. Fitting Heisenberg Hamiltonians to the measured dispersion of the triplet excitations confirmed that the dominant second neighbor exchange coupling is crucial to account for the shape of $`\mathrm{\Omega }_Q`$.

The complete Heisenberg Hamiltonian for CaV<sub>4</sub>O<sub>9</sub> has four different coupling constants: nearest-neighbor ($`nn`$) and next-nearest-neighbor ($`nnn`$) couplings and, for each of these, intra- and inter-plaquette couplings. In notation of Gelfand et al., the Hamiltonian is given by

$`H`$ $`=`$ $`J_1{\displaystyle \underset{nn}{}}𝐒_i𝐒_j+J_1^{}{\displaystyle \underset{nn^{}}{}}𝐒_i𝐒_j`$ (1)   
$`+`$ $`J_2{\displaystyle \underset{nnn}{}}𝐒_i𝐒_j+J_2^{}{\displaystyle \underset{nnn^{}}{}}𝐒_i𝐒_j,`$ (2)   
where $`𝐒_i`$ denotes the spin $`\frac{1}{2}`$ operator in site $`i`$. The $`nn`$ sums run over nearest-neighbor bonds and the $`nnn`$ sums run over next-nearest-neighbor bonds. Unprimed sums connect V’s in the same plaquette, while primed sums connect V’s in different plaquettes. The four couplings are drawn in different line styles in Fig. 1.

In this Letter we show that the spin gap behavior of CaV<sub>4</sub>O<sub>9</sub>, even considering its complex structure with eight very low symmetry V<sup>4+</sup> ions in the primitive cell, can be calculated in ab initio fashion. Our work has three separate aspects. 1) Local spin density approximation (LSDA) calculations are used to obtain energies for various magnetic configurations. The resultant exchange interactions are obtained by fitting these energies to the mean-field Heisenberg model as described below. 2) An approximate but physically motivated local orbital method called the self-consistent atomic deformation (SCAD) method  is used to provide explicit local orbitals, eigenvalues, and hopping integrals for calculating the exchange interactions from perturbation theory. This method reveals that the $`nn`$ interactions are not intrinsically small, but the net value of the superexchange coupling is small due to cancellations among various fourth-order processes. It also indicates that direct V-V exchange coupling is important. 3) The Heisenberg Hamiltonian is solved using exact diagonalization techniques on finite periodic clusters. Spin gap behavior is obtained, and $`\chi `$(T) is similar to the data. The Heisenberg couplings are refined by fitting to $`\chi (T)`$. The resulting Hamiltonian agrees well with $`\chi (T)`$ and with the triplet dispersion determined from neutron scattering.

The LSDA calculations of the energy for various magnetic configurations were more precise extensions of previous work on CaV<sub>4</sub>O<sub>9</sub> . The magnetism of the V ion is found to be robust, allowing us to break the spin symmetry in any manner we choose and obtain the energy from a self-consistent calculation. The symmetry of the non-magnetic state is initially broken as desired by applying the necessary local magnetic fields to the V ions. The seven configurations we have chosen include the ferromagnetic (FM) state, one ferrimagnetic (FiM) state, and five antiferromagnetic (AF) states with zero net spin. These AF states include the Néel state, a state in which FM plaquettes are antialigned (FMPL), and a state in which the metaplaquettes are aligned antiferromagnetically (AFMP). The configurations, given explicitly in Table I, were chosen either because of their physical relevance (AFMP was anticipated to be lowest in energy, as found) or computational considerations such as retaining inversion symmetry.

The resulting energies were fit to the mean-field Heisenberg model, which contains simply the $`\mathrm{S}_i^z`$ or Ising terms of the full Hamiltonian (2), to determine the four coupling constants. The six energy differences lead to six conditions on the four $`J`$s, and a least-squares fit gives the values listed as LSDA in Table II, each with a fitting uncertainty of about 1 meV. Since both nearest and next nearest couplings are AF in sign, there is a great deal of frustration in the magnetic system. The large value of $`J_2^{}`$ indicates that singlet formation on the metaplaquette is the driving force for the spin gap.

To understand how these values of the exchange parameters arise, we evaluate the fourth-order expressions for the exchange constants, using an approximate but parameter-free method based on the SCAD method. For each coupling constant in CaV<sub>4</sub>O<sub>9</sub>, we focus on the relevant clusters for each coupling. The $`nnn`$ interactions require a V<sub>2</sub>O cluster with two V ions (each with one relevant orbital) and one O in between. The $`nn`$ exchange interactions require a V<sub>2</sub>O<sub>2</sub> cluster. All three $`2p`$ orbitals in each O are relevant, since the low symmetry makes them non-degenerate and oriented in directions determined not by symmetry but by electronic interactions.

We neglect the Hubbard $`U`$ and Hund’s rule coupling on the O ions. In what follows, $`U`$ is the V on-site repulsion, $`ϵ_V`$ and $`ϵ_\alpha `$ are site energies of the V and $`\alpha `$-th O orbitals, and $`t_{i\alpha }`$ is the hopping amplitude between the $`i`$-th V and the $`\alpha `$-th O orbital. Defining the energy denominators $`\mathrm{\Delta }_\alpha =U+ϵ_Vϵ_\alpha `$ simplifies the expressions.

The initial state has each O orbital doubly filled and each V with one electron. The perturbation theory is given by three fourth-order terms and the direct second-order V-V term:

$`J`$ $`=`$ $`j_1+j_2+j_3+j_d`$ (3)   
$`=`$ $`{\displaystyle \frac{4}{U}}\left({\displaystyle \underset{\alpha }{}}{\displaystyle \frac{t_{1\alpha }t_{2\alpha }}{\mathrm{\Delta }_\alpha }}\right)^2+4{\displaystyle \underset{\alpha }{}}{\displaystyle \frac{(t_{1\alpha }t_{2\alpha })^2}{\mathrm{\Delta }_\alpha ^3}}`$ (4)   
$`+`$ $`4{\displaystyle \underset{\alpha <\beta }{}}{\displaystyle \frac{t_{1\alpha }t_{2\alpha }t_{1\beta }t_{2\beta }}{\mathrm{\Delta }_\alpha +\mathrm{\Delta }_\beta }}\left({\displaystyle \frac{1}{\mathrm{\Delta }_\alpha }}+{\displaystyle \frac{1}{\mathrm{\Delta }_\beta }}\right)^2+{\displaystyle \frac{4t_{12}^2}{U}}`$ (5)   
In the $`nnn`$ case, $`\alpha `$ and $`\beta `$ sum over the three orbitals in the single oxygen atom. In the $`nn`$ case, $`\alpha `$ and $`\beta `$ sum over the six orbitals in both oxygen atoms. The first three terms in (5) can be categorized by their configurations after the second hop of the four-hop process: 1) One vanadium empty; 2) One oxygen orbital empty; 3) Two oxygen orbitals half filled. The last term has an extra factor of two because it arises twice: the total spin singlet case is reduced in energy and the total spin triplet is increased by the same amount. The latter picks up a minus sign due to electron exchange.

This expression is evaluated with the SCAD model, which expresses the total density $`n(r)`$ as a sum over localized densities $`|\varphi _\alpha ^{(i)}(𝐫𝐑_i)|^2`$ centered at the atomic sites $`𝐑_i`$ . The orbitals $`\varphi _\alpha ^{(i)}`$ are solutions to atom-centered one-electron Hamiltonians $`H_i`$ for each site. The potentials in $`H_i`$ are determined self-consistently from the expression for the functional derivative of the total energy. It includes a local approximation for exchange and correlation energy and the Thomas-Fermi function for kinetic energy of overlapping densities.

Each V ion has the lowest of its five 3d levels occupied by a single electron, giving the V<sup>4+</sup>, O<sup>2-</sup> ionic description. $`U3.5`$ eV was computed by minimizing the SCAD energy subject to the constraint that one V ion has its charge increased by unity. The electron comes mainly from the other V ions with only a minor portion coming from the nearby O ions.

The matrix elements, $`t_{ij}=\psi _i|H|\psi _j`$ require the full Hamiltonian $`H`$ and orthogonalized orbitals $`\psi `$. The $`\psi `$’s are obtained from the SCAD orbitals using Löwdin’s method , and $`H`$ is determined from the site centered SCAD Hamiltonians by removing the kinetic energy overlap contributions from the latter. This gives expressions for $`H`$ that differ in the site selected for spherical harmonic expansion of the potential. We find the two possibilities, $`t_{ij}`$ and $`t_{ji}`$, may differ by $``$20%, which leads to a much larger uncertainty in the fourth-order $`J`$’s. Since the vanadium sites of a given pair of V ions are equivalent by symmetry, the direct interaction, $`j_d`$, has no such uncertainty. To be consistent with the direct interaction calculation, we use the vanadium-site-expanded potentials for evaluating matrix elements between oxygen-vanadium pairs. The net values obtained (labelled SCAD in Table II) agree rather well with those derived from LSDA energies for $`J_1`$, $`J_2`$, and $`J_2^{}`$. The close agreement may be fortuitous in view of the uncertainties mentioned above and the approximations inherent in the SCAD method. Nevertheless, we believe certain qualitative features of the SCAD results are real: 1) The values for $`J_1`$ and $`J_1^{}`$ result primarily from $`j_d`$, with relatively small contributions from fourth-order terms due to cancellation within $`j_1`$ and between $`j_2`$ and $`j_3`$. 2) The value for $`J_2^{}`$, the largest coupling, is dominated by a single term in $`j_1`$, resulting from V overlap with the middle O 2p level.

For each set of four coupling constants, we calculated the uniform susceptibility of the Hamiltonian (2) on periodic 20-spin clusters. The susceptibility is given by:

$$\chi (T)=\frac{n(g\mu _B)^2}{Nk_BT}\underset{ij}{}S_i^zS_j^z,$$
(6)   
where $`n`$ is the number of V atoms per gram and $`N`$ is the number of sites in the cluster. We take $`g=1.67`$ for all plots. This was determined from the fit to the experimental magnetic susceptibility data described below.

To evaluate (6), we calculate all eigenvalues of the Hamiltonian—eigenvectors are not required. We block-diagonalize the Hamiltonian with all possible symmetries: translations, rotations, S, and $`S_z`$ . The blocks are left with no degeneracies, so the eigenvalues are calculated very efficiently using the Lanczos algorithm with no reorthogonalization developed by Cullum and Willoughby. This allows $`\chi `$ to be calculated exactly at all temperatures using one Lanczos run for each symmetry sector. The Hamiltonian for the 20-spin cluster has blocks as large as 36950. Within each block at least the 400 lowest and highest eigenvalues are calculated, and an analytic density of states is assumed for the middle eigenvalues. This technique will be described elsewhere.

The susceptibility of the full Hamiltonian (2) calculated with each set of coupling constants in Table II is shown in Fig. 2. The experimental susceptibility of Taniguchi, et al.  is shown for comparison. All curves exhibit a spin gap, as evidenced by their low temperature behavior, $`\chi (T0)\mathrm{exp}(\mathrm{\Delta }/kT)`$, where $`\mathrm{\Delta }`$ is the gap. Both the LSDA and SCAD approaches overestimate the gap, indicating that the calculated coupling constants are too large. The coupling constants deduced from neutron scattering are also shown .

Also shown in Fig. 2 is a curve generated using the coupling constants obtained from a least-squares fit of the susceptibility to the experimental results. In the fitting procedure, we allow the $`g`$-value in eq. (6) and all four $`J`$’s to vary. At the best fit, we obtain the coupling constants listed as “Fit” in Table II and shown as the line thicknesses in Fig. 1. We find $`g=1.67`$, which is smaller than the $`g`$-value indicated by ESR measurements . Near the minimum, the fitting function is quadratic. The eigenvalues of the Hessian (scaled by an arbitrary constant) are 1, 0.046, 0.013, and 0.00039. The smallness of the last eigenvalue indicates that in the $`\delta \{J_1,J_1^{},J_2,J_2^{}\}=\{0.09,0.57,0.81,0.09\}`$ direction from the minimum, the least-squares fit is very soft.

The 20-spin cluster is sufficiently large compared with the correlation length to describe the infinite system accurately. The minimum triplet gap hardly varies between 20 and 32-spin clusters: $`\mathrm{\Delta }_{20}=9.92`$ meV while $`\mathrm{\Delta }_{32}=10.02`$ meV for the Fit Hamiltonian.

Fig. 3 shows the triplet dispersion $`\mathrm{\Omega }_Q`$ of the LSDA, SCAD, and susceptibility-fit coupling constants calculated with the expansion in Ref. . Since the LSDA and SCAD coupling constants overestimate the gap, we rescaled their $`J`$’s by 0.58 and 0.65, respectively. Both the Fit and rescaled LSDA $`\mathrm{\Omega }_Q`$ agree with the neutron scattering data reasonably well; in particular, they correctly have minima at $`Q=(0,0)`$.

To conclude, we have shown that the quantum-disordered phase in CaV<sub>4</sub>O<sub>9</sub> can be predicted in ab initio fashion. We calculated the coupling constants of the Heisenberg Hamiltonian for CaV<sub>4</sub>O<sub>9</sub> in two very different first-principles approaches. In both methods, the strongest coupling is found between next-nearest-neighbor V atoms on metaplaquettes—the weak coupling between nearest-neighbor V’s results from the cancellation among superexchange processes. The uniform magnetic susceptibility for each set of coupling constants is calculated using a novel finite-temperature exact diagonalization technique, which shows the Hamiltonians determined from both ab initio approaches have quantum-disordered phases. The Hamiltonian that best fits the experimental susceptibility is calculated, and the agreement is remarkable. Finally the triplet dispersion of the ab initio and best susceptibility-fit Hamiltonians are shown to agree well with the neutron scattering data.

We thank Z. Weihong for the code to calculate the curves in Fig. 3 and N.E. Bonesteel, J.L. Feldman, R.E. Rudd, M. Sato, R.R.P. Singh, and C.C. Wan for stimulating conversations. This work was supported by the Office of Naval Research. C.S.H was supported by the National Research Council, and W.E.P. by NSF Grant DMR-9802076. Computations were done at the Arctic Region Supercomputing Center and at the DoD Major Shared Resource Centers at NAVOCEANO and CEWES.