---
abstract: 'We report a study of exchange interactions in bulk CrO$_2$ calculated from first principles. We considered three near neighbor Cr-Cr exchange interactions: the interaction between corner and body center atoms mediated through a single oxygen atom; the interaction between a Cr and the Cr directly “above” it in the (001) direction, also mediated by a single O atom; and the interaction between a Cr and its neighbor in the (100) direction, mediated by two intervening oxygen atoms. The interactions were calculated by rotating the moments of one or more of the Cr ions while constraining the others to remain parallel. We then fit the resulting energy vs. angle data to the Heisenberg model and extracted exchange energy parameters with a least-squares method. We also calculated the exchange interactions using a “spin-spiral” technique, in which a relative angular displacement was imposed upon Cr moments in adjacent cells. Similar results were obtained with both approaches. The calculated $T=0$ K exchange interactions were subsequently used to determine the magnetization as a function of temperature via low-T spin-wave dispersion and a Monte-Carlo method. Reasonable agreement with experiment was obtained.'
author:
- 'H. Sims'
- 'S. J. Oset'
- 'W. H. Butler'
- 'James M. MacLaren'
- Martijn Marsman
title: 'Exchange Interactions in CrO$_2$'
---

Introduction
============

CrO$_2$ is one of only a few known ferromagnetic oxides, and is the only material which has been experimentally shown to be a ferromagnetic “half-metal”[@DFT-CrO2; @ji]. A half-metal is a material which is a metal for one spin channel and an insulator for the other. CrO$_2$ crystallizes in the rutile crystal structure (Figure \[rutile\]), as do TiO$_2$, VO$_2$, MnO$_2$, RuO$_2$, and SnO$_2$. The existence of isostructural oxides with a variety of different electronic and magnetic properties makes the rutile system interesting for theoretical investigations of spintronics because one can envisage the growth of layered devices with the same crystal structure throughout. Since CrO$_2$ offers such special opportunities for understanding oxide spintronics, it is important to establish how well our standard electronic structure tools work in dealing with the electronic and magnetic structure of this material. It is well known that they encounter difficulties in dealing with many transition metal oxides, including the very similar oxide VO$_2$, which DFT also predicts to be a half-metal at 0K, [@VO2paper] but is observed to be an insulator. An additional motiviation for understanding exchange interactions in CrO$_2$ is the fact that its Curie temperature ($T_c=386.5$ K)[@curie; @curie2] is sufficiently close to room temperature that its magnetic properties are significantly degraded at room temperature, hindering potential spintronics applications. A better understanding may point the way to improvement.

Electronic Structure of CO$_2$ Within Density Functional Theory
===============================================================

The electronic structure and density of states were calculated using density functional theory[@Kohn-Sham] (DFT). Our calculated density of states is similar to previous calculations.[@DFT-CrO2] The electronic structure of CrO$_2$ can be understood by comparing the DOS of CrO$_2$ with that of TiO$_2$. Rutile TiO$_2$ is a non-magnetic insulator. In the rutile structure, as in many transition metal oxides, each transition metal atom is surrounded by six oxygen nearest-neighbors. Each oxygen atom has three transition metal (TM) nearest neighbors. The TM ions form a body-centered tetragonal lattice, with a $\frac{c}{a}$ ratio of approximately $\frac{2}{3}$. The oxygen ions form rows collinear with the Cr ions. The rows containing the body-center Cr ions can be considered to run in the (110) direction, in which case the rows containing the corner Cr ions run in the (1$\overline{1}$0) direction. In addition to the lattice constants ($a=b$ and $c$), a single internal parameter is sufficient to describe the structure. This parameter may be taken to be the distance in the $a$-$b$ plane between a Cr ion and its neighboring oxygen ions. The oxygen atoms in the top and bottom faces of the tetragonal cell form a rectangle with a Cr ion at the center. The $\frac{c}{a}$ ratio and the internal parameter usually conspire such that the environment of each Cr ion is an approximate octahedron of oxygen ions. However, the octahedron is distorted because the square base of the octahedron with the Cr ion at its center is actualy a rectangle. In addition, the distance from a Cr ion to the oxygen ions at the apices of the octahedron may differ from the distance to the four equatorial oxygen ions that form the rectangular base. As shown in Figure \[rutile\], the edges of distorted octahedra form “ribbons” along the (001) direction.

It is straightforward to show that if we treat this system in a tight-binding approximation in which the TM atoms only interact directly with the oxygen atoms (i. e. hopping matrix elements only connect nearest neighbors), there will be an energy gap separating the oxygen $p$-states and the TM $d$-states. The gap extends from the O-$p$ onsite energy to the O-$d$ onsite energy. This gap is apparent in TiO$_2$, for which the oxygen $p$-states are filled and the Ti $d$-states are empty (Fig. \[DOSTiO2\]). When an energy gap occurs at the Fermi energy, it contributes significantly to reducing the energy of the structure, because all occupied states are pushed down in energy, while all unoccupied states are pushed up. In CrO$_2$, there are two additional electrons per TM atom compared to TiO$_2$, so some of the $d$-states above the gap must be occupied.

![Rutile structure. For CrO$_2$, we use $a=4.42$ and $\frac{c}{a}\approx0.670$ (experimental parameters).[]{data-label="rutile"}](structure.eps)

![Density of States for rutile TiO$_2$[]{data-label="DOSTiO2"}](TiO2.DoS.eps)

Let us consider three possibilities for the occupation of the $d$-states in CrO$_2$. The first possibility is that the system remains non-magnetic, similar to TiO$_2$ (Fig. \[DOSnmCrO2\]). In this case each Cr atom will have one spin-up and one spin-down $d$-electron. From Fig. \[DOSnmCrO2\] we can see that there remains a gap between the O-$p$ states and the Cr-$d$ states, but it is significantly smaller than for TiO$_2$. If the interactions are strictly nearest neighbor connecting only O ions with Ti or Cr ions, then the O-$p$ on-site energy in a tight-binding model would be at the top of the O-$p$ states and the transiton metal $d$ on-site energy would be at the bottom of the $d$-states. It should be understood that both states below the gap and the states above the gap are hybrids of O-$p$ and TM-$d$. However, the lower set of states are predominantly O-$p$, and there are 3 of them per spin channel per O ion, and the higher lying states are predominantly TM-$d$, and there are 5 of them per spin channel per TM ion.

![Density of States for Nonmagnetic CrO$_2$.[]{data-label="DOSnmCrO2"}](CrO2.NM.DoS.eps)

A second possibility is that the system becomes ferromagnetic, the experimentally-observed ground state. In this case there will be $1+x$ majority spin electrons and $1-x$ minority spin electrons (Fig. \[DOSferroCrO2\]). In this spin configuration, energy is gained from the intra-atomic “Hund’s first rule” energy, which in most versions of DFT is primarily a potential proportional to the local magnetization. Additional energy can be saved if $x=1$, because this restores the gap for the minority spin channel. It is clear from Fig. \[DOSferroCrO2\] that the system chooses $x=1$ for the ferromagnetic case. The system does this by lowering the on-site energy of the majority of Cr $d$-states and raising the on-site energy of the minority Cr $d$-states. If we were confident that there were only nearest neighbor interactions, the majority density of states would imply that the Cr majority $d$ on-site energy was essentially degenerate with the O-$p$ on-site energy. In reality, there are next-nearest-neighbor interactions that spread the O-$p$ states above the O-$p$ on-site energy and spread the Cr-$d$ states below the Cr-$d$ onsite energy.

![Density of States for Ferromagnetic CrO$_2$.[]{data-label="DOSferroCrO2"}](CrO2.FM.DoS.eps)

A third possibility is to maintain two electrons of the same spin on each Cr atom, but to make the system antiferromagnetic — for example, by having opposing moments on the corner and body-centered TM atoms of the 6-atom primitive cell (Fig. \[DOSAFCrO2\]). In this case, the spin-dependent DOS of the corner and body-center Cr ions are mirror images of each other. Assuming that the corner Cr ion moment is “up,” its $d$ on-site energy for spin up is pushed down while its $d$ on-site energy for spin down is pushed up. The opposite will happen for the body-center Cr ion. The result is that there is no gap in the [*total*]{} DOS for either spin channel, but if one looks only at the DOS for one of the Cr ions, one can see an approximate gap in one of the spin channels. Thus for this anti-ferromagnetic configuration, there would be a relatively large density of states at the Fermi energy for a given spin channel, but only on alernate Cr ions.

![Density of States for Antiferromagnetic CrO$_2$.[]{data-label="DOSAFCrO2"}](CrO2.AF.DoS.eps)

Comparing these three possibilities, it is not surprising that the ferromagnetic state has the lowest energy, the nonmagnetic state the highest (1.02 eV above ferromagnetic) with the anti-ferromagnetic intermediate between the two (0.30 eV above ferromagnetic). Thus, the tendency to form a moment in CrO$_2$ is very strong, and the energy associated with the ferromagnetic alignment of moments based on this initial test is moderately large within DFT. It should be recognized that other more complicated spin arrangements (e. g. , different antiferromagnetic states) may have lower energy than the simple one calculated here.

Exchange Interactions in CO$_2$
===============================

In order to investigate interatomic exchange interactions in CrO$_2$ in more detail, we have calculated the near-neighbor exchange interactions along the (100), (001), and (111) directions by rotating moments within specially-constructed supercells. We fit the resulting relationship between the energy of the system and the angle of rotation to the Heisenberg model $$H = \sum_{i,j}J_{ij}\boldsymbol{\mu}_i \cdot \boldsymbol{\mu}_j
\label{eq:heis}$$ where $\left|\boldsymbol{\mu}\right| = g\mu_B S = 2\mu_B$ is the spin moment, $g$ is the electron spin $g$-factor, $S = 1$ is the spin number, and $\mu_B$ is the Bohr magneton. To make contact with the standard Heisenberg model, we can pull the magnitude of the spin moment ($2\mu_B$) into the value of $J$ and treat the spins as unit vectors.

In addition to this supercell approach, we have taken advantage of a recently developed feature in the Vienna [*Ab-initio*]{} Simulation Package[@VASP] (VASP) to calculate a so-called helimagnetic state in which the moment in the $n^{\mathrm{th}}$ magnetic layer is canted by an angle $n\theta$ with respect to the $0^{\mathrm{th}}$ layer. In so doing, we are able to calculate several orders of $J_n$ of the form $$E = E_0 + \sum_nJ_n\cos{}n\theta$$ via Fourier analysis.

Near Neighbor Exchange Using Supercells
---------------------------------------

All of the calculations in this study were performed within DFT[@Kohn-Sham] in the generalized gradient approximation[@Perdew-Wang] (GGA) and in the local (spin) density approximation with onsite Coulomb interactions (LSDA+U)[@LDApU] using the Dudarev method,[@theDude] for which we use $U-J$ = 3 eV. We perform all calculations using the VASP software[@VASP] and pseudopotentials generated by Kresse [*e*t al]{}.[@Kresse]. To calculate the near-neighbor exchange interactions, we created a supercell containing two rutile unit cells (using both experimental and DFT-relaxed lattice parameters), stacked in either the (100) (Fig. \[100fig\]) or (001) (Fig. \[001fig\]) direction as appropriate. In all of the following calculations, we use an energy cut-off of 500 eV. For cells stacked along the (100) direction, we use a $5\times9\times15$ Monkhorst-Pack[@MP] grid of k-points, a $9\times9\times7$ grid for supercells stacked along (001), and a $9\times9\times15$ grid for the 6-atom cell used in the spin-spiral calculations. We also make use of the spin interpoloation method of Vosko-Wilk-Nusair.[@VoskoWN] Each of the 12-atom supercells has four Cr ions, whose magnetic moments we can individually constrain within the calculation. We chose three distinct magnetic configurations, designed to probe the exchange coefficients. In the first configuration, we rotated the moment of a corner Cr atom and held all other moments fixed using the constraining field method in VASP. In the second, we rotated the two Cr moments in the centers of their respective unit cells, and in the third we rotated a corner atom and its nearest center atom. A summary of the configurations can be found in Table \[exchtab\].

![The (100) supercell, with Cr ions numbered for comparison to Table \[exchtab\]. The bottom half of the structure in Figure \[rutile\] has been projected onto the (010) plane. []{data-label="100fig"}](cro2100s1.eps)

![The (001) supercell, with Cr ions numbered for reference. The left half of the structure in Figure \[rutile\] has been projected onto the (010) plane.[]{data-label="001fig"}](cro2001s1.eps)

                    [**Cr$_1$**]{}   [**Cr$_2$**]{}   [**Cr$_3$**]{}   [**Cr$_4$**]{}
  ---------------- ---------------- ---------------- ---------------- ----------------
  [**Case 1**]{}        fixed           rotated           fixed            fixed
  [**Case 2**]{}       rotated          rotated           fixed            fixed
  [**Case 3**]{}       rotated           fixed           rotated           fixed

  : Magnetic configurations used to calculated exchange coupling.[]{data-label="exchtab"}

To ensure that the systems would be sufficiently well-behaved, we rotated the moments through small angles (up to 60$^\circ$). We fit the energy vs. angle data to $A(1-\cos\theta)+B$, where $A$ is equivalent to the combined exchange parameter $J$ and $B$ is simply the angle-independent component of the energy.

For a given choice of supercell orientation, we have the following system of equations for the supercell method: $$A_{\textrm{Case 1}}=8J_{111}+2J_{100/001}
\label{eq:100fit}$$ $$A_{\textrm{Case 2}}=16J_{111}$$ $$A_{\textrm{Case 3}}=8J_{111} + 4J_{100/001}$$ Using a least-squares technique for overdetermined systems of equations,[@Williams] we can write $$AJ=b
\label{eq:lsquare}$$ $$A^TAJ=A^Tb$$ $$\overline{J}=(A^TA)^{-1}A^Tb$$ $$\sigma = \left|A\cdot\overline{J}-b\right|$$ where $\overline{J}$ is the calculated $J$ column vector, $\sigma$ is the error in the fit, and $$A=
\left(\begin{array}{cc}
8 & 2 \\
16 & 0 \\
8 & 4
\end{array}\right)\qquad
J=
\left(\begin{array}{c}
J_{111} \\
J_{100/001}
\end{array}\right)$$

We summarize the calculations performed within GGA and LSDA+U for experimental and relaxed lattice parameters using both the supercell method in Table \[tab:everything\] Throughout this work, the terms “experimental” and “relaxed” as they appear in this table denote structures with the experimental and the GGA- or LSDA+U-relaxed lattice parameters, respectively.

  ----------------- ----- --------------- ---------------- -------------- --------------
                                           ($U-J = 3$ eV)                 
                     $z$   Experimental       Relaxed       Experimental     Relaxed
  $J_{100}$ (meV)     4    $-11.8\pm2.5$   $-10.4\pm0.7$    $1.4\pm0.5$    $1.1\pm1.0$
  $J_{001}$ (meV)     2    $33.8\pm5.6$     $33.8\pm5.0$    $32.8\pm0.3$   $31.5\pm0.3$
  $J_{111}$ (meV)     8    $23.2\pm6.1$     $22.9\pm5.0$    $24.5\pm0.4$   $24.6\pm1.0$
  ----------------- ----- --------------- ---------------- -------------- --------------

  : Summary of all calculated exchange energies obtained using the supercell method. Uncertainties given arise from the error in the least-squares fit. In $J_{111}$, there is some (usually negligible) contribution to the error from the standard deviation of the values obtained through (100)- and (001)-stacked supercells. The quantity $z$ is the coordination number for the given interaction—the number of such interactions acting on a given Cr ion. Note that the (100) and (010) directions are equivalent and are referred to as (100) throughout this work. The columns “experimental” and “relaxed” refer to the experimental or DFT-relaxed structure.[]{data-label="tab:everything"}

The results of the calculations for the three cases are summarized as follows: in each case, we find a near-perfect fit to the cosine function, provided that we restrict the fit to small angles (less than or equal to $60^\circ$), as we did with the original calculations. We can see clearly in Table \[tab:everything\] the anisotropic nature of the exchange, which is to be expected given the shape of the cell. Most interestingly, we find that the interaction between Cr neighbors along the (100) or (010) directions (parallel to the $a$ or $b$ axes) is antiferromagnetic. However, the strength and multiplicity of the other interactions is sufficient to lead to a ferromagnetic ground state. Considering the dependence on lattice parameter, we notice that the (001) and (111) interactions seem to be almost unchanged with the small (0.6%) change in lattice constant. Somewhat surprisingly, however, the (100) interaction (calculated within the GGA) increases (becomes more positive) by more than an meV under this small expansion of the lattice. We also note that the LSDA+U calculations predict a positive, though nearly vanishing, $J_{100}$.

Helimagnetism
-------------

Helimagnetism is a noncollinear magnetic state in which the spins in adjacent layers along a certain direction are rotated with respect to one another by a fixed angle. Rutile MnO$_2$, for example, has been shown to exhibit helimagnetic ordering in the ground state.[@MnO2] We do not suspect that CrO$_2$ is a helimagnetic material, but by setting up a helimagnetic spin state, we can investigate the exchange using a different approach. The recently-added spin spiral capabilities of VASP[@spiral] allow us to calculate arbitrarily long-range exchange interactions within bulk CrO$_2$.

The spin spiral method modifies the periodic boundary conditions of the supercell approach, imposing helimagnetic order on the magnetic structure as determined by the propagation vector [**q**]{}. The vector [**q**]{} and the angle $\phi$ between any two spins are given by $$\phi=\mathbf{q}\cdot\mathbf{r}_{j}$$ $$\mathbf{q}=\frac{2\pi}{a_i}\xi\hat{\mathbf{e}}_i
\label{qeq}$$ where the azimuthal angle $\theta$ is restricted to $\frac{\pi}{2}$. We choose the unit vector $\hat{\mathbf{e}}_i$ to be either the (100) or (001) direction, and allow $\xi$ to vary between 0 and 1. Clearly, when $\xi = 0$, we recover the ferromagnetic state.

Because the unit cell contains two magnetic ions, varying the angles between neighboring layers of CrO$_2$ requires that one modify both $\xi$ and the initial orientation of the magnetic moments. For example, in this work, we use $\xi = 1$ in conjunction with an antiferromagnetic configuration of the two moments in the unit cell to create a periodic magnetic system wherein the moment of each atom’s nearest neighbor is antiparallel it. To obtain a system in which neighboring magnetic “layers” are oriented at an angle of $\frac{\pi}{4}$ from one another, we use $\xi=\frac{1}{4}$, so that each cell after the initial one is rotated by $\frac{\pi}{2}$. We then set up the moments in the initial cell such that the corner and body-centered Cr moments are oriented at the desired angle of $\frac{\pi}{4}$, leading to a smooth spin wave in the desired direction.

In this work, we choose a relatively short, commensurate spin wavelength in order to simplify the analysis, although the method allows for more general configurations as well. Using different values of [**q**]{}, and thus different values of $\theta$, we create a [**q**]{} spectrum. We then use Fourier analysis to extract the $J_n$. These $J_n$ differ in meaning from the $J$s calculated using the supercell method; for example, $J_1 = 8 J_{111}$ and $J_2 = 2 J_{100/001}$.

To calculate the helimagnetic state, we used a supercell composed of a single rutile unit cell. The angle of each subsequent Cr ion with respect to the first is given by . After acquiring five points (including the zero-frequency point q=0) of the $E(q)$ curve, we performed a discrete Fourier transform to obtain the first 4 $J_n$. We find good agreement between the $J_1$ calculated with (100) and (001) spin spirals, as expected. We also find a difference in sign between $J_2$ in the (100) and (001) cases, in agreement with the larger supercell calculations. Moreover, this method yields the additional parameters $J_3$ and $J_4$, corresponding to $J_{211/112}$ and $J_{200/002}$, respectively. These higher-order energies are smaller than the first- and second-order exchange energies, and will be neglected in further analysis. The results of the calculations are summarized in Table \[tab:heli\].

  ------------------------- -------------- --------- -------------- ---------
                                                                    
                             Experimental   Relaxed   Experimental   Relaxed
  [**$J_{100}$**]{} (meV)       -12.0        -12.2        -1.3        -1.8
  [**$J_{001}$**]{} (meV)        27.5        29.8         29.6        26.1
  [**$J_{111}$**]{} (meV)        20.8        20.7         25.6        25.0
  ------------------------- -------------- --------- -------------- ---------

  : Summary of calculated exchange interactions (in meV) using the spin spiral method (compare with Table \[tab:everything\]). The effect of the change in lattice parameter is decidedly smaller (perhaps negligible in most cases) in the spin spiral method. Note that in this method the LSDA+U calculation also predicts a negative $J_{100}$, albeit an order of magnitude smaller than in the GGA.[]{data-label="tab:heli"}

Comparison with Experiments
===========================

Spin Wave Stiffness
-------------------

In an effort to benchmark our calculations against known experimental results, we have calculated the spin wave stiffness constant for CrO$_2$ using the expression derived by Schlottmann:[@schlott] $$D_{100} = 2(J_{111} + J_{100}) S a^2 \label{eq:100D}$$ $$D_{001} = 2(J_{111} + J_{001}) S c^2\label{eq:001D}$$ where $a$ and $c$ are the lattice spacings in the appropriate directions and $S$ is the spin number (1 for CrO$_2$). These expressions can be easily understood as anisotropic extensions of results obtained for magnons in a one-dimensional chain (for which $D=2JSa^2$). In his work, Schlottmann considers the spins as quantum operators such that $\mathbf{S}_i\cdot\mathbf{S}_j=S(S+1)$, and he keeps this value separate from $J$. Additionally, he neglects $J_{100}$ in his expression for $D_{100}$. However, we use classical spins of magnitude $2 \mu_B$ (although the units are collapsed into the exchange constant $J$ as previously explained). Consequently, we must scale our $J$s by $1/\left|\boldsymbol{\mu}\right|^2 = 1/4$ in order to apply this expression. Further, our calculations indicate that $J_{100}$ it is roughly of the same order as $J_{111}$ and $J_{001}$, so we have included it in our analysis. Using this model, we calculate $D_{100}$ and $D_{001}$ for the various cells, exchange-correlation approximations, and methods considered throughout this work. Table \[tab:stiff\] reviews the values we obtained. Experimentally, we find several mostly-consistent values for the spin wave stiffness obtained through different methods. All of the experimental values assume an isotropic stiffness constant. Ji *et al.*[@stiff] fit the M(T) curve in order to obtain the coefficient on the $T^{3/2}$ term, from which they determine $D=1.8\times10^{-40}$ Jm$^2$. Zou *et al*[@zou] used magnetic force microscopy to determine the length and width of domain walls in CrO$_2$, from which they were able to calculate $D = 2.62\times10^{-40}$ Jm$^2$. Further, Rameev *et al.*[@fmr] used ferromagnetic resonance to measure the bulk magnon modes and obtained $D_B=3\times10^{-10}$ Oe cm$^2$, which is equivalent to $D=0.57\times10^{-40}$ Jm$^2$ via the relation $D_B = 2A / \mu_0 M_s$,[@units] which is a bit smaller than other reported values.

  --------------------- --------------------------------------------- -------------- --------- -------------- ---------
                                                                                                              
                                                                       Experimental   Relaxed   Experimental   Relaxed
  [**Supercell**]{}      $\bm{D_{100}}$ $(10\times10^{-40}$ J m$^2)$      $1.81$      $1.96$       $4.06$      $3.94$
                         $\bm{D_{001}}$ $(10\times10^{-40}$ J m$^2)$      $3.91$      $3.87$       $3.91$      $3.79$
                         $\bm{D_{avg}}$ $(10\times10^{-40}$ J m$^2)$      $2.34$      $2.46$       $4.00$      $3.89$
  [**Spin Spiral**]{}    $\bm{D_{100}}$ $(10\times10^{-40}$ J m$^2)$      $1.38$      $1.35$       $3.80$      $3.56$
                         $\bm{D_{001}}$ $(10\times10^{-40}$ J m$^2)$      $3.29$      $3.46$       $3.76$      $3.45$
                         $\bm{D_{avg}}$ $(10\times10^{-40}$ J m$^2)$      $1.84$      $1.84$       $3.79$      $3.52$
  --------------------- --------------------------------------------- -------------- --------- -------------- ---------

  : Comparison of the calculated spin stiffness constants $D$ for different methods of first-principles calculation. Here, $D_{avg} = \left(D_{100}\sqrt{D_{001}}\right)^{2/3}$.[]{data-label="tab:stiff"}

Using our calculated $D$s, we can predict the low-temperature spin-wave contribution to the magnetization as a function of temperature. A relatively straight-forward generalization of the argument found in Kittel[@Kittel] for a cubic system allows one to write the spin-wave dispersion relation for small excitations and long wavelengths as $$\omega(k,k_z) = D k^2 + D_z k_{z}^{2},\qquad k^2 = k_{x}^{2} + k_{y}^{2}, \qquad D = D_{100}, \qquad D_z = D_{001}
\label{eq:disp}$$ Integrating over a surface of constant $\omega$ in $k$-space, one obtains a density of states given by $$N(\omega) = \frac{1}{4\pi^2}\frac{1}{D\sqrt{D_z}}\sqrt{\omega}$$ Using this expression and the Planck distribution, we can calculate the coefficient $B$ in the $T^{3/2}$ model $$M(T) = M(0)(1-BT^{3/2})$$ $$B=\frac{0.0587}{SQ}\frac{1}{2S(J_{100}+J_{111})}\frac{1}{\sqrt{2S(J_{001}+J_{111})}}k_{B}^{3/2} = \frac{0.0587}{SQ}\frac{V}{D\sqrt{D_z}}k_{B}^{3/2}$$ where $Q$ is the number of magnetic ions per unit cell (2, in this case), $V$ is the volume of the cell, and $k_B$ is Boltzmann’s constant. Fitting the experimental[@stiff] M(T) curve yields $B=5\times10^{-5}$ K$^{-3/2}$. Using the spin-wave stiffnesses shown in Table \[tab:stiff\], we have, for supercells, $B_{GGA}^{expt}=2.40\times10^{-5} $K$^{-3/2}$, $B_{GGA}^{rel}=2.27\times10^{-5}$ K$^{-3/2}$, $B_{LSDA+U}^{expt}=1.07\times10^{-5}$ K$^{-3/2}$, and $B_{LSDA+U}^{rel}=1.09\times10^{-5}$ K$^{-3/2}$. For the spin spiral approach, $B_{GGA}^{expt}=3.43\times10^{-5}$ K$^{-3/2}$, $B_{GGA}^{rel}=3.48\times10^{-5}$ K$^{-3/2}$, $B_{LSDA+U}^{expt}=1.16\times10^{-5}$ K$^{-3/2}$, and $B_{LSDA+U}^{rel}=1.27\times10^{-5}$ K$^{-3/2}$. Thus, the coefficient obtained from GGA is within a factor of two ,while that derived from LSDA+U is off by about a factor of four. Assuming that DFT overestimates each exchange energy equally, this implies that our calculated values of $J$ may differ from experimental values by about 50% for GGA and a factor of about $2.5$ for LSDA+U (with $U-J=3$ eV). In each case, the spin spiral numbers are slightly closer to experiment. Figure \[fig:magnon\] shows the low-T M(T) curves from the calculated spin-wave dispersion compared to that from a fit to the experimental M(T) curve.

![The low-temperature M(T) curve. The GGA Spin Spiral and LSDA+U Supercell curves represent the extremes of the range of calculated M(T) curves. We compare against the actual experimental data[@curie] and a low-T fit to these data.[]{data-label="fig:magnon"}](magnon.eps)

Curie Temperature
-----------------

In light of the favorable agreement between calculated and experimental spin stiffness, we subsequently attempted to calculate the magnetic ordering temperature of CrO$_2$, comparing a mean field prediction to Monte Carlo simulations. A mean-field model using the calculated exchange parameters yields a Curie temperature several times larger than the measured value of 386.5 K.[@curie; @curie2] The mean-field expression is given by $$k_BT = \frac{3}{2}J_{tot}$$ where $J_{tot}$ is equivalent to half of the energy difference between a ferro- and an antiferromagnetic configuration in a 6-atom (2-Cr) cell. Using this expression, we obtain a mean-field Curie temperature for CrO$_2$ of 1160 K or 1240 K for the experimental and DFT-relaxed lattice parameters in the supercell method, respectively. This is somewhat surprising given the above analysis of our estimation of the exchange. However, it is not sufficient to consider only the low-temperature behavior, especially when the system in question has a magnetic ordering temperature above room temperature. In order to gain a simple yet ideally illuminating picture of the temperature dependence, we utilized a Monte Carlo simulation using the Metropolis-Hastings algorithm[@metropolis] with random numbers generated using the Mersenne Twister method[@mersenne]. For this simulation, we used a cubic grid of $10\times10\times10$ unit cells, where a unit cell consists of a corner and body-centered Cr ion. Only Cr ions are considered, and they are treated as simple constant-magnitude magnetic moments that initially point along the z axis. Our calculations indicate that the constant-magnitude approximation should be valid as long as the angle between adjacent moments is less than $100^\circ$.

We begin with a random spin configuration with the spin vectors chosen to be uniformly distributed on the unit sphere. In the Metropolis method, an iteration consists of a randomly chosen Cr ion being assigned a magnetic moment in a random direction. This will result in a change in energy $\Delta$E from the old configuration. If $\Delta$E is negative, meaning the new energy is lower, the new direction for that moment is kept. Otherwise, the new direction still has a probability of $e^{-\Delta{}E/k_BT}$ of being kept in its new orientation to simulate thermal agitation. If neither condition for keeping the moment’s new direction is met, then the change is undone, and the lattice of spins remains unmodified until the next iteration.

The calculation of $\Delta$E at each step considers all nearest neighbors along (100), (010), (001), and (111) directions, using a Heisenberg interaction between moments with the calculated exchange constants for GGA and LSDA+U with experimental and DFT-relaxed lattice parameters. Figure \[heiscurie\] shows the simulated results for the magnitude of the net magnetization versus temperature compared to reported values.[@curie] When interpreting these data, one must must be cognizant of the fact that the Monte Carlo simulations exhibit several shortcomings—namely, that it will necessarily not be able to predict the correct low-temperature T-dependence (as it uses a classical model), that there exists an unphysical tail on the curve arising from finite-size effects in the lattice, and that we assume that exchange remains constant with temperature, likely leading to an overestimation of the Curie temperature. The errors in the shape of the curve at low temperature should not have an impact in the accuracy of the result, as each value of $k_BT$ is run independently. Further, the high-temperature tail can be accounted for assuming we know where to look for it. The remaining discrepancy, that the exchange will reduce in strength as temperature rises, is a limitation of exploring this behavior with first-principles calculations. Considering these issues, we can attempt to scale the experimental curve to match the calculated one, as we do in Fig. \[heiscurie\]. By matching the portions of the curves at lower temperatures, we can determine with some confidence where the curve ought to go to zero and thus the true prediction of the model. By so doing, we find that the Monte Carlo simulation predicts a Curie temperature of between 450 K (for the GGA spin spiral method) and 650 K (for the LSDA+U supercell method), which, even in the latter case, is in much better agreement with experiment than the mean-field model.

![The calculated M(T) behavior predicts a Curie temperature of 450-550 K, depending on which exchange parameters one chooses from among those presented in this work. We compare with the experimental data[@curie].[]{data-label="heiscurie"}](MCs1.eps)

Conclusions
===========

We have calculated the near neighbor exchange interactions for bulk CrO$_2$ in the (100), (001), and (111) directions, finding them to be $-12.0$ meV, $24.5$ meV, and $23.6$ meV respectively. From our calculated spin stiffness parameters and the results of our classical Heisenberg Metropolis method, we obtain some confidence that DFT and VASP can determine the exchange interaction of CrO$_2$ with reasonable accuracy. Examining the calculated exchange parameters, we find that the sign of $J_{100}$, both in the supercell and the equivalent spin spiral calculations, indicates the possibility of non-collinear behavior in CrO$_2$ if the exchange parameters are modified. Thus, a mixed interface between CrO$_2$ and another material (such as RuO$_2$) might lead to non-collinear spins if the ratio between nearest and next-nearest neighbor interactions is pushed into a “favorable” zone. We investigate this possibility explicitly for CrO$_2$-RuO$_2$ interfaces in a companion paper by K. Chetry *et al*[@ruo2].

Acknowledgments
===============

This work was supported by the NSF-DMR under Grant No. 0706280 and MRSEC Grant No. 0213985. It was completed using computing resources from the University of Alabama’s High Performance Cluster.

[00]{}

K. Schwarz, J. Phys. F: Met. Phys. [**16**]{}, L211-215 (1986); I.I. Mazin, D.J. Singh, and C. Ambrosch-Draxl, Phys. Rev. B [**59**]{}, 411-418 (1999).

Y. Ji, G. S. Strijkers, F. Y. Yang, C. L. Chien, J. M. Byers, A. Anguelouch, G. Xiao, and A. Gupta, Phys. Rev. Lett. **86**, 5585-5588 (2001); J. S. Parker, S. M. Watts, P. G. Ivanov, and P. Xiong, Phys. Rev. Lett. **88**, 196601 (2002).

M. E. Williams, W. H. Butler, C. K. Mewes, H. Sims, M. Chshiev, and S. K. Sarker J. Appl. Phys. [**105**]{}, 07E510 (2009).

X. W. Li, A. Gupta, Giang Xiao, Appl. Phys. Lett. [**75**]{}, 713 (1999).

F. Y. Yang, C. L. Chien, X. W. Li, Gang Xiao, A. Gupta, Phys. Rev. B [**63**]{}, 092403 (2001).

G. Kresse and J. Hafner, Phys. Rev. B [**47**]{}, 558 (1993); G. Kresse and J. Hafner, Phys. Rev. B [**49**]{}, 14 251 (1994); G. Kresse and J. Furthmüller, Comput. Mat. Sci. [**6**]{}, 15 (1996); G. Kresse and J. Furthmüller, Phys. Rev. B [**54**]{}, 11169 (1996).

W. Kohn and L. J. Sham, Phys. Rev. [**140**]{}, A1133 (1965).

Y. Wang and J. P. Perdew, Phys. Rev. B [**43**]{}, 8911 (1991).

V. I. Anisimov, J. Zaanen, O. K. Andersen, Phys. Rev. B [**44**]{}, 943 (1991).

S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B [**57**]{}, 1505 (1998).

G. Kresse, and J. Joubert, Phys. Rev. B [**59**]{}, 1758 (1999).

H. J. Monkhorst and J. D. Pack, Phys. Rev. B [**13**]{}, 5188 (1976).

S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys. [**58**]{}, 1200 (1980).

M. Marsman and J. Hafner, Phys. Rev. B **66**, 224409 (2002).

Gareth Williams, [*The American Mathematical Monthly*]{}, Vol. 97, No. 6. (Jun. – Jul., 1990), pp. 511-513.

Stephen Blundell, [*Magnetism in Condensed Matter*]{}. Oxford University Press, New York. 2001.

A. Yoshimori, J. Phys. Soc. Jpn. [**14**]{}, 807 (1959).

P. Schlottmann, Phys. Rev. B. [**67**]{}, 174419 (2003).

C. Kittel, [*Introduction to Solid State Physics*]{}, 8th Ed. Hoboken, NJ: John Wiley and Sons, Inc., 2005.

Y. Ji, G. J. Strijkers, F. Y. Yang, C. L. Chien, J. M. Byers, A. Anguelouch, G. Xiao, A. Gupta, Phys. Rev. Lett. [**86**]{}, 5585 (2001).

Xiaozing Zou and Gang Xiao, Appl. Phys. Lett. [**91**]{}, 113512 (2007).

B. Z. Rameeva, R. Yilgina, B. Akta, A. Gupta and L. R. Tagirov, Microelec. Eng. [**69**]{}, 336 (2003).

J. Hamrle, O. Gaier, Seong-Gi Min, B. Hillebrands, Y. Sakuraba, and Y. Ando, J. Phys. D: Appl. Phys. [**42**]{}, 084005 (2009).

Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, and Augusta H. Teller, J. Chem. Phys. [**21**]{}, 1087 (1953); W.K. Hastings, Biometrika [**57**]{}, 97 (1970).

Makoto Matsumoto and Takuji Nishimura, ACM Transactions on Modeling and Computer Simulation. [**8**]{}, 3 (1998).

K. Chetry, H. Sims, Zhihong Lu, W. H. Butler, and A. Gupta, unpublished.
