Spin-Orbit Excitons in CoO

TL;DR

This study explores the complex magnetic excitations in CoO, revealing how strong spin-orbit coupling and exchange interactions lead to entangled spin-orbit excitons, with a model capturing low-energy features but failing at high energies.

Contribution

It introduces a multi-level spin exciton model for CoO that incorporates spin-orbit coupling and exchange interactions, extending previous theoretical frameworks.

Findings

The model reproduces low-energy magnetic excitations.

High-energy excitations show faster decay than the model predicts.

Failure at high energies suggests breakdown of localized spin-orbit excitons.

Abstract

CoO has an odd number of electrons in its unit cell, and therefore is expected to be metallic. Yet, CoO is strongly insulating owing to significant electronic correlations, thus classifying it as a Mott insulator. We investigate the magnetic fluctuations in CoO using neutron spectroscopy. The strong and spatially far-reaching exchange constants reported in [Sarte et al. Phys. Rev. B 98 024415 (2018)], combined with the single-ion spin-orbit coupling of similar magnitude [Cowley et al. Phys. Rev. B 88, 205117 (2013)] results in significant mixing between $j_{eff}$ spin-orbit levels in the low temperature magnetically ordered phase. The high degree of entanglement, combined with the structural domains originating from the Jahn-Teller structural distortion at $\sim$ 300 K, make the magnetic excitation spectrum highly structured in both energy and momentum. We extend previous theoretical work on PrTl$_{3}$ [Buyers et al. Phys. Rev. B 11, 266 (1975)] to construct a mean-field and multi-level spin exciton model employing the aforementioned spin exchange and spin-orbit coupling parameters for coupled Co$^{2+}$ ions on a rocksalt lattice. This parameterization, based on a tetragonally distorted type-II antiferromagnetic unit cell, captures both the sharp low energy excitations at the magnetic zone center, and the energy broadened peaks at the zone boundary. However, the model fails to describe the momentum dependence of the excitations at high energy transfers, where the neutron response decays faster with momentum than the Co$^{2+}$ form factor. We discuss such a failure in terms of a possible breakdown of localized spin-orbit excitons at high energy transfers.

Tables

Table 1. Table 1: Initial values (in meV) for the parameters of the spin-orbit exciton model.

Parameter	Initial Value	Reference
$\lambda$	$-16$	Cowley et al.,2013
$\Gamma$	$-8.76$	Buyers et al., 1971,Gladney,1966
$J_{1F}$	$-0.918$	Sarte et al.,2018a
$J_{1AF}$	1.000
$J_2$	3.09
$J_{3F}$	$-0.182$
$J_{3AF}$	0.262
$J_{4F}$	$-0.0504$
$J_{4AF}$	0.0759
$H_{MF}$	64.8	Kanamori, 1957,Singer, 1956,Nagamiya et al.,1955

Table 2. Table 2: Refined values (in meV) of the mean molecular field parameter H M ​ F subscript 𝐻 𝑀 𝐹 H_{MF} for all 16 x ​ A ​ x ​ x ​ γ 𝑥 𝐴 𝑥 𝑥 𝛾 xAxx\gamma orbital configurations considered in a mean-field multi-level spin-orbit exciton model, each with refined values of λ 𝜆 \lambda , Γ Γ \Gamma , and J m , ξ subscript 𝐽 𝑚 𝜉 J_{m,\xi} listed in Tab.  3 . Numbers in parentheses indicate statistical errors.

Orbital Configuration	Refined Value
AAAA1	62.4(2)
AAAA2	46.2(1)
AAAF1	55.3(2)
AAAF2	53.9(1)
AAFA1	46.2(1)
AAFA2	49.9(1)
AAFF1	56.2(2)
AAFF2	56.0(2)
FAAA1	47.3(1)
FAAA2	47.9(1)
FAAF1	58.9(3)
FAAF2	58.8(2)
FAFA1	61.6(3)
FAFA2	60.9(3)
FAFF1	48.9(1)
FAFF2	59.5(2)
Average	54.4(4)

Table 3. Table 3: Summary of the initial values, parameter spaces, and refined values for the parameters of the mean-field multi-level spin-orbit exciton model. All values are reported in meV and numbers in parentheses indicate statistical errors.

Parameter	Initial Value	Range	Refined Value
$\lambda$	$-16$	[$-19$,$-13$]	$-19.00\left(1\right)$
$\Gamma$	$-8.76$	[$-8.76$,$-6.16$]	$-6.16\left(1\right)$
$J_{1F}$	$-0.918$	[$-1.134$,$-0.730$]	$-0.780\left(1\right)$
$J_{1AF}$	1.000	[0.798,1.24]	0.848(1)
$J_2$	3.09	[2.29,4.55]	2.43(1)
$J_{3F}$	$-0.182$	[$-0.220$,$-0.145$]	$-0.154\left(1\right)$
$J_{3AF}$	0.262	[0.209,0.316]	0.223(1)
$J_{4F}$	$-0.0504$	[$-0.0581$,$-0.0402$]	$-0.0428\left(1\right)$
$J_{4AF}$	0.0759	[0.0606,0.0874]	0.0645(1)
$H_{MF}$	64.8	[0,100]	Tab. 2

Equations

$$S(\mathbf{Q},\omega )=g_L^2{f}^2(\mathbf{Q})\sum _{\alpha \beta }({\delta }_{\alpha \beta }-{\hat{Q}}_{\alpha }{\hat{Q}}_{\beta })S^{\alpha \beta }(\mathbf{Q},\omega ),$$

$$S^{\alpha \beta }(\mathbf{Q},\omega )=\frac{1}{2\pi }\int dt{e}^{i\omega t}⟨{\hat{S}}^{\alpha }(\mathbf{Q},t){\hat{S}}^{\beta }(-\mathbf{Q},0)⟩,$$

$$S^{\alpha \beta }(\mathbf{Q},\omega )=-\frac{1}{\pi }\frac{1}{1-\exp (\omega /k{}_{\mathrm{B}}T)}\Im G^{\alpha \beta }(\mathbf{Q},\omega ).$$

$$\begin{array}{rl}{G}^{\alpha \beta }(ij,t)& =G({\hat{S}}^{\alpha }(i,t),{\hat{S}}^{\beta }(j,0))\\ & =-i\Theta (t)⟨[{\hat{S}}^{\alpha }(i,t),{\hat{S}}^{\beta }(j,0)]⟩,\end{array}$$

$$\omega G(\hat{A},\hat{B},\omega )=⟨[\hat{A},\hat{B}]⟩+G([\hat{A},\hat{\mathcal{H}}],\hat{B},\omega ),$$

$$\begin{array}{rl}\hat{\mathcal{H}}& ={\hat{\mathcal{H}}}_{CF}+\sum _{ij}J(ij)\hat{\mathbf{S}}(i)\cdot \hat{\mathbf{S}}(j).\end{array}$$

$${\hat{\mathcal{H}}}_{MF}(i)=\sum _i{H}_{MF}(i){\hat{S}}_z(i),$$

$$H_{MF}(i)=2\sum _{i>j}J(ij)⟨{\hat{S}}_z(j)⟩,$$

$${\hat{\mathcal{H}}}_1=\sum _i{\hat{\mathcal{H}}}_{CF}(i)+\sum _i{\hat{S}}_z(i)\left(2\sum _{i>j}J(ij)⟨{\hat{S}}_z(j)⟩\right),$$

$$\begin{array}{rl}{\hat{\mathcal{H}}}_2& =\sum _{ij}J(ij){\hat{S}}_z(i)[{\hat{S}}_z(j)-2⟨{\hat{S}}_z(j)⟩]\\ & +\frac{1}{2}\sum _{ij}J(ij)[{\hat{S}}_{+}(i){\hat{S}}_{-}(j)+{\hat{S}}_{-}(i){\hat{S}}_{+}(j)],\end{array}$$

$$⟨{\hat{S}}_{\alpha }⟩=\sum _n{f}_n⟨n|{\hat{S}}_{\alpha }|n⟩\equiv \sum _n{\hat{S}}_{\alpha nn}{f}_n,$$

$${\hat{\mathcal{H}}}_1=\sum _n\sum _i{\omega }_n{C}_n^{†}(i)C_n(i),$$

$$\begin{array}{rl}{\hat{\mathcal{H}}}_1& ={\hat{\mathcal{H}}}_{CF}+{\hat{\mathcal{H}}}_{MF}=\\ & ({\hat{\mathcal{H}}}_{CEF}+{\hat{\mathcal{H}}}_{SO}+{\hat{\mathcal{H}}}_{dis})+{\hat{\mathcal{H}}}_{MF},\end{array}$$

$${\hat{\mathcal{H}}}_{CEF}=B_4\left({\hat{\mathcal{O}}}_4^0+5{\hat{\mathcal{O}}}_4^4\right).$$

$${\hat{\mathcal{H}}}_{SO}=\lambda \hat{\mathbf{L}}\cdot \hat{\mathbf{S}},$$

$$\mathcal{C}=\left[{\begin{array}[]{ccccccc}0&0&-0.79&0.61&0&0&0\\ 0&0&0&0&-0.71&0&-0.71\\ 0.61&0&0&0&0&-0.79&0\\ 0&1.00&0&0&0&0&0\\ 0&0&-0.61&-0.79&0&0&0\\ 0&0&0&0&-0.71&0&0.71\\ 0.79&0&0&0&0&0.61&0\\ \end{array}}\right].$$

$$\mathcal{C}^{-1}\hat{L}_{z}\mathcal{C}=\left[{\begin{array}[]{ccc|ccc|c}\mathbf{1.50}&\mathbf{0}&\mathbf{0}&0&0&-1.94&0\\ \mathbf{0}&\mathbf{0}&\mathbf{0}&0&0&0&0\\ \mathbf{0}&\mathbf{0}&\mathbf{-1.50}&-1.94&0&0&0\\ \hline\cr 0&0&-1.94&\mathbf{-0.50}&\mathbf{0}&\mathbf{0}&0\\ 0&0&0&\mathbf{0}&\mathbf{0}&\mathbf{0}&2.00\\ -1.94&0&0&\mathbf{0}&\mathbf{0}&\mathbf{0.50}&0\\ \hline\cr 0&0&0&0&2.00&0&\mathbf{0}\\ \end{array}}\right],$$

$$\hat{L}_{z}=\left[{\begin{array}[]{ccc}-1&0&0\\ 0&0&0\\ 0&0&1\\ \end{array}}\right]$$

$${\hat{\mathcal{H}}}_{SO}=\alpha \lambda \hat{\mathbf{l}}\cdot \hat{\mathbf{S}},$$

$$diag\left(\hat{\mathcal{H}}_{SO}\right)=\left[\begin{array}[]{cccccccccccc}-60&0&0&0&0&0&0&0&0&0&0&0\\ 0&-60&0&0&0&0&0&0&0&0&0&0\\ 0&0&-24&0&0&0&0&0&0&0&0&0\\ 0&0&0&-24&0&0&0&0&0&0&0&0\\ 0&0&0&0&-24&0&0&0&0&0&0&0\\ 0&0&0&0&0&-24&0&0&0&0&0&0\\ 0&0&0&0&0&0&36&0&0&0&0&0\\ 0&0&0&0&0&0&0&36&0&0&0&0\\ 0&0&0&0&0&0&0&0&36&0&0&0\\ 0&0&0&0&0&0&0&0&0&36&0&0\\ 0&0&0&0&0&0&0&0&0&0&36&0\\ 0&0&0&0&0&0&0&0&0&0&0&36\\ \end{array}\right],$$

$$\mathcal{C}^{-1}\hat{j_{z}}\mathcal{C}=\left[\begin{array}[]{cc|cccc|cccccc}\mathbf{-\frac{1}{2}}&\mathbf{0}&0&0&0&0&0&0&0&0&0&0\\ \mathbf{0}&\mathbf{\frac{1}{2}}&0&0&0&0&0&0&0&0&0&0\\ \hline\cr 0&0&\mathbf{-\frac{3}{2}}&\mathbf{0}&\mathbf{0}&\mathbf{0}&0&0&0&0&0&0\\ 0&0&\mathbf{0}&\mathbf{-\frac{1}{2}}&\mathbf{0}&\mathbf{0}&0&0&0&0&0&0\\ 0&0&\mathbf{0}&\mathbf{0}&\mathbf{\frac{1}{2}}&\mathbf{0}&0&0&0&0&0&0\\ 0&0&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{\frac{3}{2}}&0&0&0&0&0&0\\ \hline\cr 0&0&0&0&0&0&\mathbf{-\frac{5}{2}}&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\ 0&0&0&0&0&0&\mathbf{0}&\mathbf{-\frac{3}{2}}&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\ 0&0&0&0&0&0&\mathbf{0}&\mathbf{0}&\mathbf{-\frac{1}{2}}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\ 0&0&0&0&0&0&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{\frac{1}{2}}&\mathbf{0}&\mathbf{0}\\ 0&0&0&0&0&0&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{\frac{3}{2}}&\mathbf{0}\\ 0&0&0&0&0&0&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{\frac{5}{2}}\\ \end{array}\right].$$

$${\hat{\mathcal{H}}}_{dis}=B_2{\hat{\mathcal{O}}}_2^0=\Gamma \left({\hat{l}}_z^2-\frac{2}{3}\right),$$

$${\hat{\mathcal{H}}}_{MF}=\sum _i{H}_{MF}(i){\hat{S}}_z=2z_2{J}_2⟨\widehat{S_z}⟩{\hat{S}}_z,$$

$${\hat{\mathcal{H}}}_{pair}=\alpha \lambda {\hat{\mathbf{l}}}_1\cdot {\hat{\mathbf{S}}}_1+\alpha \lambda {\hat{\mathbf{l}}}_2\cdot {\hat{\mathbf{S}}}_2+2J_{1,2}{\hat{\mathbf{S}}}_1\cdot {\hat{\mathbf{S}}}_2,$$

$${\hat{S}}_{(\pm ,z)}=\sum _{mn}{\hat{S}}_{(\pm ,z)mn}{C}_m^{†}{C}_n.$$

$$G^{\alpha \beta }(i,j,\omega )=\sum _{mn}{\hat{S}}_{\alpha mn}{\hat{G}}^{\beta }(m,n,i,j,\omega ),$$

$$g^{\alpha \beta }(E)=\sum _n\left\{\frac{S_{\alpha 0n}{S}_{\beta n0}}{E+i\Delta -E_{n0}}-\frac{S_{\alpha n0}{S}_{\beta 0n}}{E+i\Delta +E_{n0}}\right\},$$

$$J(\mathbf{Q})=\sum _{i\ne j}{J}_{ij}{e}^{i\mathbf{Q}\cdot {\mathbf{d}}_{ij}},$$

$$G^{\alpha \beta }(\mathbf{Q},E)=g^{\alpha \beta }(E)+g^{\alpha +}(E)J(\mathbf{Q})G^{-\beta }(\mathbf{Q},E)+g^{\alpha -}(E)J(\mathbf{Q})G^{+\beta }(\mathbf{Q},\omega )+2g^{\alpha z}(E)J(\mathbf{Q})G^{z\beta }(\mathbf{Q},E),$$

$$G^{+-}(\mathbf{Q},E)$$

Full text

I Introduction

Mott insulators are materials where conventional band theory fails, predicting metallic behavior owing to half-filled bands, with the origin of the insulating response indicative of strong electronic correlations Anisimov et al. (1991); Mott (1949); Phillips (2006); Brandow (1976, 1977). Mott insulators are parent materials for high-temperature cuprate superconductivity Lee et al. (2006); Kastner et al. (1998); Zhao et al. (2011). Moreover, there have been some suggestions that these insulators may even be implicated as being parent to some iron-based superconductors Si and Abrahams (2008); Zhu et al. (2010); McCabe et al. (2014); Stock and McCabe (2016); Wilson et al. (2010). These Mott insulators display well-defined spin excitations, however rapidly breakdown Plumb et al. (2014); Headings et al. (2010); Stock et al. (2010, 2007) on charge doping towards superconductivity. Ho et al. (2001) More recently, Mott insulators with strong spin-orbit coupling have been of particular interest in the search for unconventional topological states Moore (2010); Schaffer et al. (2016); Witczak-Kremp et al. (2014). These studies have focussed on 4$d$ and 5$d$ transition metals with strong spin-orbit coupling resulting in $j\rm{{}_{eff}}=\frac{1}{2}$ ground states, and new Kitaev bond directional phases Jackeli and Khaliullin (2009); Winter et al. (2017); Takagi et al. (2019). However, much of the single-ion physics that results in these $j\rm{{}_{eff}}=\frac{1}{2}$ ground states is present in Co2+-based compounds that also display strong spin-orbit coupling Sano et al. (2018); Liu and Khaliullin (2018).

In this context, it is timely to investigate the classic Mott insulator CoO, where significant spin-orbit coupling is present and comparable to the magnetic exchange. In this study, we investigate the mixed spin-orbit transitions in CoO through their parameterization with a multi-level spin-orbit exciton model extending previous theoretical work on PrTl3 Buyers et al. (1975). While this model reproduces the experimental data at low-energy transfers, we show that its failure at high energy transfers is accompanied by a possible breakdown of these excitations.

For the past several decades, CoO has been one of the most extensively studied Mott insulators. The 3$d$ metal monoxide was among the first orbitally-ordered materials to be investigated with neutron diffraction Shull et al. (1951). Its primitive unit cell consists of one 3$d^{7}$ Co2+ and one 2$p^{6}$ O2- with 2 electrons in the n=2 shell, which corresponds to 15 valence electrons. With an odd number of valence electrons, conventional band theory Norman (1989); Terakura et al. (1984) would predict CoO to be metallic. However, CoO is a very strong insulator Engel and Schmid (2009) with a room temperature resistivity of 108 $\Omega\cdot\rm{cm}$, and an optical band gap of 2.5 eV van Elp et al. (1991); Gillen and Robertson (2013); Gvishi and Tanhauser (1972), with evidence for metallic behavior being found only under extremely high pressures on the order of 100 GPa Cohen et al. (1997).

Possessing a cubic $Fm\bar{3}m$ structure Jauch et al. (2001); Sasaki et al. (1979); Tombs and Rooksby (1950); Ok and Mullen (1968); Liu et al. (2006) at room temperature (Fig. 1(a)), CoO assumes long-range antiferromagnetic order at $T\rm{{}_{N}}\sim 290$ K Greenwald (1953), in contrast to the long-range ferromagnetism predicted by general band coupling models that assume a dominant direct exchange Deng et al. (2010); Walsh et al. (2008). Despite being the subject of many neutron diffraction studies, its magnetic structure has proven to be particularly contentious, with both collinear Roth (1958); Ressouche et al. (2006) and non-collinear van Laar (1965); Tomiyasu et al. (2004) models describing diffraction patterns equally well Fiebig et al. (2003); van Laar et al. (1966); Kernavanois et al. (2003); Shishidou and Jo (1998); Timm et al. (2016).

As illustrated in Fig. 1(b), the ${}^{4}T_{1}$ crystal field ground state for CoO corresponds to the $d^{7}$ Co2+ assuming a high spin ($S=\frac{3}{2}$) configuration, yielding an orbital triplet with one hole in the $t_{2g}$ orbital manifold. The resulting orbital degeneracy, coupled with both a Jahn-Teller driven Jauch et al. (2001) unit cell distortion and various far-reaching large exchange interactions, yields a complex magnetic excitation spectrum that results from the strong entanglement of multiple spin-orbit levels (Fig. 1(c)). The resulting multi-parameter spin-orbital Hamiltonian incorporating both exchange and spin-orbit coupling of similar magnitude, further complicated by the complex magnetic ordering and structural distortions, has made the understanding of the magnetic excitations in this material particularly difficult Sakurai et al. (1968).

By employing both the spin-orbit coupling constant $\lambda$ and the magnetic exchange constants $J$ that were experimentally determined in our previous work Cowley et al. (2013); Sarte et al. (2018a) on the magnetically diluted monoxide Mg0.97Co0.03O, we will show that the low energy magnetic excitation spectrum of CoO measured in the Néel regime by inelastic neutron spectroscopy is reproduced by a mean-field multi-level spin-orbit exciton model based on Green’s functions Buyers et al. (1975). Our parameterization successfully captures the fine structure of well-defined low energy spin excitations present at the magnetic zone center and also the broadening in both momentum and energy at the zone boundaries. In contrast, the model fails to reproduce the high energy response consisting of high velocity excitations that decay with momentum faster than the Co2+ form factor. We suggest that this failure of the model provides evidence for a breakdown of localized spin-orbit excitations at high energy transfers, possibly replaced by delocalized or itinerant-like fluctuations, despite CoO being a strong Mott insulator.

This paper is divided into three general sections. In the first section, we describe the theoretical framework that we apply to parameterize the neutron scattering response in CoO. We first outline the single-ion response defining the crystal field Hamiltonian, and then discuss the coupled equations-of-motion that were used to numerically derive the neutron response. In the second section, we first present the experimental data as measured with neutron spectroscopy, followed by a direct comparison to our multi-level spin-orbit exciton model. To conclude, we discuss the high energy excitations and the poor agreement with the Co2+ form factor and speculate as to their origin.

II Theory

We first discuss the theoretical framework used to describe the localized magnetic response in CoO. The neutron magnetic cross section is proportional to the magnetic dynamic structure factor $S(\bf{Q},\omega)$ defined by

[TABLE]

where $g_{L}$ is the Landé $g$-factor, $f({\bf{Q}})$ is the magnetic form factor, Kernavanois et al. (2003, 2004) and $S^{\alpha\beta}(\bf{Q},\omega)$ corresponds to the dynamic spin structure factor. Since the orbital angular momentum is quenched ($i.e.$ $\left\langle\hat{\mathbf{L}}\right\rangle=0$) for $3d^{7}$ Co2+, the orbital contribution to the scattering cross section is assumed to be weak, and therefore the spin operators provide the dominant contribution to the neutron scattering cross section Zaliznyak and Lee (2005). This assumption allows $S^{\alpha\beta}(\bf{Q},\omega)$ to be defined in terms of expectation values of spin operators $\hat{S}^{\nu}(i,t)$ of index $\nu=$ $+$, $-$, or $z$, acting on a site $i$ at a time $t$. Such a definition of the dynamic structure factor is given by

[TABLE]

whose relationship to the response function $G^{\alpha\beta}(\bf{Q},\omega)$ is given by the fluctuation-dissipation theorem as

[TABLE]

Motivated by previous work on PrTl3 Buyers et al. (1975); Fulde and Peschel (1971), the theoretical portion of this paper begins by first writing down the equations-of-motion for the response function in terms of commutators involving the magnetic Hamiltonian $\hat{\mathcal{H}}$. We then investigate this magnetic Hamiltonian in CoO to define both the single-ion states and how these states are coupled from site to site on the rocksalt lattice. Finally, we apply mean-field theory to decouple the equations-of-motion, thereby reducing the formula for the response functions to a set of coupled linear equations that can be computed numerically and directly compared with experiment. In this approach, we use creation and annihilation operators of the single-ion states rather than the Holstein-Primakoff transformation for a single spin operator. This approach allows for both the incorporation of spin-orbit level mixing, and the explicit inclusion of the single-ion terms in the Hamiltonian, such as spin-orbit coupling, rather than employing anisotropy terms that incorporate the orbital physics through perturbation theory Yosida (1991).

II.1 The Equation-of-Motion for the Response Function

According to linear response theory, the response function measured with neutrons is proportional to the Fourier transform of the retarded Green’s function that is given by

[TABLE]

where $\Theta(t)$ is the Heaviside function. As shown in Section I of the Supplementary Information sup , by taking the first time derivative of $G^{\alpha\beta}(ij,t)$, applying the Heisenberg equation-of-motion, and Fourier transforming from time to energy, one arrives at the following equation-of-motion

[TABLE]

where $\hat{A}$ and $\hat{B}$ denote generic spin operators. Eq. 4 indicates that deriving a model for the neutron scattering response functions relies both on the understanding of the Hamiltonian $\hat{\mathcal{H}}$ and its commutator with the spin operators. We now investigate the individual contributions to magnetic Hamiltonian in CoO.

II.2 The Total Magnetic Hamiltonian $\hat{\mathcal{H}}$

The total magnetic Hamiltonian consisting of crystal field ($CF$) contributions and coupling between Co2+ on sites $i$ and $j$ can be written as

[TABLE]

By defining a molecular field Hamiltonian

[TABLE]

where

[TABLE]

$\hat{\mathcal{H}}$ can be written as a sum of a single-ion ($\hat{\mathcal{H}}_{1}$) and an inter-ion ($\hat{\mathcal{H}}_{2}$) term given by

[TABLE]

and

[TABLE]

where $\langle\hat{S}_{z}(j)\rangle$ denotes a thermal average given by

[TABLE]

with $\hat{S}_{\alpha nn}=\langle n|\hat{S}_{\alpha}|n\rangle$, and $f_{n}$ is the Boltzmann thermal population factor. The inclusion of a factor of 2 in Eq. 6 follows the convention that was established in Ref. Sakurai et al., 1968 to explicitly account for the double counting in the sum over sites.

The procedure we follow to derive the neutron response consists of two parts. First, we diagonalize the single-ion component $\hat{\mathcal{H}}_{1}$ for a given molecular field such that

[TABLE]

where $C(i)$ and $C^{\dagger}(i)$ are ladder operators satisfying the commutation relations $[C_{n}(i),C^{\dagger}_{m}(j)]=\delta_{ij}\delta_{nm}$, and $\omega_{n}$ are the energy eigenvalues. The second step consists of using these states to apply mean-field theory on the inter-ion $\hat{\mathcal{H}}_{2}$ term to then compute the neutron response using the equation-of-motion given in Eq. 4. We will now discuss the eigenstates of the single-ion component of the Hamiltonian $\hat{\mathcal{H}}$, followed by the inter-ion component. This section concludes with the application of mean-field theory on the inter-ion component, allowing for the derivation of an expression for the Green’s function that can be calculated numerically and compared directly to experiment.

II.2.1 Single-Ion Hamiltonian $\hat{\mathcal{H}}_{1}$

As schematically illustrated in Fig. 1(c), the single-ion component of the Hamiltonian consists of four components

[TABLE]

corresponding to the contributions from the crystalline electric field $\hat{\mathcal{H}}_{CEF}$, spin-orbit $\hat{\mathcal{H}}_{SO}$, structural distortion $\hat{\mathcal{H}}_{dis}$, and mean molecular field $\hat{\mathcal{H}}_{MF}$. Hyperfine nuclear transitions are neglected since previous measurements Chatterji and Schneider (2009) have indicated that these are on the order of $\sim\mu$eV, and thus beyond the experimental resolution of our study. We now discuss each term of the Hamiltonian $\hat{\mathcal{H}}_{1}$.

The Crystalline Electric Field, $\hat{\mathcal{H}}_{CEF}$: As illustrated in the Tanabe-Sugano Tanabe and Sugano (1954a, b); Liu and Khaliullin (1963) diagram in Fig. 1(b), for the case of CoO, the $3d^{7}$ Co2+ is octahedrally coordinated by the weak field O2- ligand resulting in a crystal field splitting $10Dq$ that is weaker than the energy differences between the free-ion terms Pratt and Coelho (1959). Consequently, the crystalline electric field contribution $\hat{\mathcal{H}}_{CEF}$ can be treated as a perturbation to the free-ion basis states that are defined by Hund’s rules incorporating the effects of electron-electron Coulomb repulsion and the Pauli Exclusion Principle. A combination of Hund’s first rule of maximum multiplicity and second rule requiring the total orbital angular momentum $L$ be maximized yields a total spin of $S=\frac{3}{2}$ and orbital angular momentum $L=3$, for the $3d^{7}$ Co2+, corresponding to an orbital ground state term symbol of ${}^{4}F$.

The crystalline electric field $\hat{\mathcal{H}}_{CEF}$ contribution corresponding to the octahedral coordination of the ${}^{4}F$ free-ion ground state by O2- ligands can be expressed in terms of the Stevens operators $\hat{\mathcal{O}}^{0}_{4}$ and $\hat{\mathcal{O}}^{4}_{4}$, and the numerical coefficient $B_{4}<0$ Abragam and Bleaney (1986); Hutchings (1964) as

[TABLE]

Since the spin-orbit coupling is expected to be considerably weaker than the crystal field contribution for the $3d$ Co2+, the complete set of commuting observables are $\hat{L}^{2}$, $\hat{L}_{z}$, $\hat{S}^{2}$ and $\hat{S}_{z}$ with corresponding good quantum numbers $L$, $m_{L}$, $s$ and $m_{s}$ in the Russell-Saunders $L$-$S$ coupling scheme; thus, by the Wigner-Eckart theorem, both Stevens operators in $\hat{\mathcal{H}}_{CEF}$ can be defined in the $|L,m_{L}\rangle$ basis, as summarized in Section II of the Supplementary Information sup .

The diagonalization of $\hat{\mathcal{H}}_{CEF}$ results in an orbital triplet ground state (${}^{4}T_{1}$), an excited orbital triplet (${}^{4}T_{2}$), and an orbital singlet ${}^{4}A_{2}$, where $\Delta(^{4}T_{1}\rightarrow^{4}T_{2})=480B_{4}$ and $\Delta(^{4}T_{2}\rightarrow^{4}A_{2})=600B_{4}$. The Stevens factor $B_{4}$ is related to the crystal field splitting by $10Dq=400B_{4}$ (Fig. 1(b)), where $10Dq$ was previously measured to be $\sim$1 eV Cowley et al. (2013); van Schooneveld et al. (2012); Magnuson et al. (2002); Chiuzbaian et al. (2008); Sakellaris and Mavridis (2012); Gorschlüter and Merz (1994). Since the ${}^{4}T_{1}$ crystal field ground state and ${}^{4}T_{2}$ first excited state are separated by $\sim$ 1 eV, it is a valid approximation that the ${}^{4}T_{1}$ ground state will exclusively determine the magnetic properties of CoO Larson et al. (2007); Haverkort et al. (2007).

Spin-Orbit Coupling, $\hat{\mathcal{H}}_{SO}$: The second perturbation to the ${}^{4}F$ free-ion ground state is spin-orbit coupling given by

[TABLE]

where $\lambda$ is the spin-orbit coupling constant. A common approach is to exclusively consider the ${}^{4}T_{1}$ ground state, requiring a projection from the original $|L=3,m_{L}\rangle$ basis onto a smaller basis $|l=1,m_{l}\rangle$ that defines the subspace that is spanned by the crystal field ground state. As discussed by Abragam & Bleaney Abragam and Bleaney (1986), this particular projection can be performed using representation theory. Here, we outline an alternate method based on the matrix representation of angular momentum operators Sarte et al. (2018b) that was inspired by the work on $4d$ and $5d$ transition metal oxides by Stamokostas and Fiete Stamokostas and Fiete (2018).

The matrix approach begins by first determining the set of eigenvectors $|\phi_{CEF}\rangle$ of the crystalline electric field Hamiltonian $\hat{\mathcal{H}}_{CEF}$ in the $|L=3,m_{L}\rangle$ basis. Since $|\phi_{CEF}\rangle$ is also a basis, a transformation matrix $\mathcal{C}$ can be constructed that rotates from the $|L=3,m_{L}\rangle$ to the $|\phi_{CEF}\rangle$ basis. The matrix $\mathcal{C}$ consists of columns corresponding to eigenvectors of $\hat{\mathcal{H}}_{CEF}$ in the $|L=3,m_{L}\rangle$ basis arranged in order of increasing energy eigenvalues. In the case of degenerate eigenvalues, a small perturbative Zeeman term of the form $\epsilon\hat{S}_{z}$, with $\epsilon$ being a small constant, was applied to remove the degeneracy and uniquely define the column order. For Co2+ in octahedral coordination with $B_{4}$ set to $-1$, $\mathcal{C}$ is given by

[TABLE]

Having obtained the transformation matrix, the rotation from the $|L=3,m_{L}\rangle$ to the $|\phi_{CEF}\rangle$ basis can then be accomplished by $\hat{\mathcal{O}}_{|\phi_{CEF}\rangle}=\mathcal{C}^{-1}\hat{\mathcal{O}}_{|L,m_{L}\rangle}\mathcal{C}$.

For the $\hat{L}_{z}$ operator, this transformation yields

[TABLE]

illustrating the ground state orbital triplet ${}^{4}T_{1}$, and the excited orbital triplet ${}^{4}T_{2}$ and singlet ${}^{4}A_{2}$ states. We note that the opposite sequence exists in the case of a tetrahedral environment with an orbital singlet ground state Decaroli et al. (2015). A comparison of the top 3 $\times$ 3 block matrix to the $\hat{L}_{z}$ operator in the $|l=1,m_{l}\rangle$ basis given by

[TABLE]

confirms that the block matrix is equivalent to the $\hat{L}_{z}$ operator in the $|l=1,m_{l}\rangle$ basis, with a projection factor $\alpha=-\frac{3}{2}$, in agreement with previous approaches based on representation theory. Therefore, in the low temperature/energy limit, the spin-orbit Hamiltonian (Eq. 10) can be rewritten as

[TABLE]

corresponding to a new Hamiltonian consisting of new orbital angular momentum operators that act on the projected $|l=1,m_{l}\rangle$ basis.

By assigning an effective angular momentum operator to the subspace spanned by the $|l=1,m_{l}\rangle$ basis, it is implied these new operators must follow the same commutation relations for general angular momentum operators. To check this fundamental requirement is satisfied, we have transformed the three $\hat{L}_{x,y,z}$ operators, each of which is a $7\times 7$ matrix, to the $|\phi_{CEF}\rangle$ basis with the procedure outline above. We then extracted the top $3\times 3$ block matrices of the projected matrices $\mathcal{C}^{-1}\hat{L}_{x,y,z}\mathcal{C}$ to define $\hat{l}_{x,y,z}$, and confirmed that these matrices do follow the commutation relations of angular momentum given by $\hat{\bf{l}}\times\hat{\bf{l}}=i\hat{\bf{l}}$. We note that the presence of a thermally isolated low energy triplet does not guarantee that these commutation relations are followed. An example of such a failure has been recently discussed in the context of low energy doublets in the heavy fermion CeRhSi3 Pásztorová et al. (2019).

Having projected $\hat{\mathbf{L}}$ onto a fictitious operator $\hat{\mathbf{l}}$ to reflect the triplet orbital degeneracy of the ${}^{4}T_{1}$ ground state of $\hat{\mathcal{H}}_{CEF}$, we now derive the eigenstates of the perturbative $\hat{\mathcal{H}}_{SO}$ term. The basis is now the 12 $|l=1,m_{l};s={3\over 2},m_{s}\rangle$ states, and based on both the Landé interval rule and the addition theorem of angular momentum, we expect this Hamiltonian to yield three levels defined by $j\rm{{}_{eff}}$=${1\over 2}$, ${3\over 2}$, and ${5\over 2}$. Using the projection factor $\alpha=-\frac{3}{2}$, and the experimentally determined Cowley et al. (2013) spin-orbit coupling constant $\lambda=-16$ meV, the diagonalization of the spin-orbit Hamiltonian $\hat{\mathcal{H}}_{SO}$ matrix yields

[TABLE]

confirming a doublet ground state with quartet and sextet excited manifolds (Fig. 1(c)) with $\Delta$(doublet$\rightarrow$quartet) = $-36$ meV $\equiv\frac{3}{2}\alpha\lambda$ and $\Delta$(quartet$\rightarrow$sextet) = $-60$ meV $\equiv\frac{5}{2}\alpha\lambda$. Finally, we may confirm that each manifold corresponds to $j\rm{{}_{eff}}=\frac{1}{2}$, $\frac{3}{2}$, and $\frac{5}{2}$, respectively, by projecting the components of the effective total angular momentum operator $\hat{\mathbf{j}}=\hat{\mathbf{l}}+\hat{\mathbf{S}}$ onto the subspaces spanned by the spin-orbit manifolds of $\hat{\mathcal{H}}_{SO}$. As was the case for the projection onto the subspace spanned by the ${}^{4}T_{1}$ crystal field ground state, such a projection is accomplished by first defining a transformation $\mathcal{C}$ which rotates operators from the $|l,m_{l},s,m_{s}\rangle$ basis to the $|\phi_{SO}\rangle$ basis with the columns being the eigenvectors of $\hat{\mathcal{H}}_{SO}$ arranged in increasing energy. Rotating the $\hat{j}_{z}$ operator onto the $|\phi_{SO}\rangle$ basis yields

[TABLE]

A comparison of the block matrices of the projected $\hat{j}_{z}$ operator as given above to the $\hat{J}_{z}$ operator in the $|j{\rm{{}_{eff}}}=\frac{1}{2},m_{j}\rangle$, $|j{\rm{{}_{eff}}}=\frac{3}{2},m_{j}\rangle$, and $|j{\rm{{}_{eff}}}=\frac{5}{2},m_{j}\rangle$ bases confirms that the top 2 $\times$ 2, middle $4\times 4$ and bottom $6\times 6$ block matrices corresponds to $j\rm{{}_{eff}}=\frac{1}{2},\frac{3}{2}$, and $\frac{5}{2}$ manifolds, respectively. By projecting the $\hat{j}_{x}$ and $\hat{j}_{y}$ operators, it can be shown that these block matrices also satisfy the canonical commutation relations of angular momentum.

The Distortion Hamiltonian, $\hat{\mathcal{H}}_{dis}$: The next perturbation to the single-ion Hamiltonian corresponds to the structural deformation of the CoO unit cell that accompanies long-range antiferromagnetic order, resulting in the distortion of the crystalline electric field from ideal octahedral coordination Schrön et al. (2012). While the exact symmetry of the low temperature phase has proven to be particularly contentious, we will consider a simple tetragonal distortion, corresponding to a uniaxial distortion along the $z$ axis. Utilizing symmetry arguments, Walter (1960) the influences of such a distortion on the crystalline electric field is given by

[TABLE]

with a distortion parameter $\Gamma$.

The Molecular Field Hamiltonian, $\hat{\mathcal{H}}_{MF}$: The final term in the single-ion Hamiltonian corresponds to the effect of the molecular field that results from the magnetic order of Co2+ moments in the rocksalt lattice. $\hat{\mathcal{H}}_{MF}$ behaves as a Zeeman-like term, resulting in a splitting of the nearly degenerate $j\rm{{}_{eff}}$ levels. By considering a single dominant next-nearest neighbor $180^{\circ}$ Co2+-O2--Co2+ superexchange pathway with a magnetic exchange constant $J_{2}$, the corresponding mean molecular field is given by

[TABLE]

where $H_{MF}(i)$ is defined by Eq. 5 above, and $z_{2}$ denotes the number of next-nearest neighbors. As summarized by Fig. 1(c), the result of such a strong value of this exchange interaction is the significant entanglement between individual $j\rm{{}_{eff}}$ levels, in contrast with other Co2+-based magnets such as CoV2O6, Wallington et al. (2015) CoV3O8, Sarte et al. (2018b) and CoNb2O6 Cabrera et al. (2014), where the degree of mixing is much weaker. In these particular magnets, the spin-orbit split levels remain well-separated in energy, and therefore can be considered as $j{\rm{}_{eff}}=\frac{1}{2}$ magnets. The strong intertwining in CoO represents a limitation imposed on approaches based on conventional linear spin wave theory, and the reason why a multi-level spin-orbit exciton model needs to be considered.

II.2.2 Inter-Ion Hamiltonian $\hat{\mathcal{H}}_{2}$

As summarized by Eq. 7, the inter-ion Hamiltonian is defined by the exchange parameters $J(ij)$ between sites $i$ and $j$. In contrast to the parameters for the single-ion Hamiltonian $\hat{\mathcal{H}}_{1}$: $10Dq$, $\lambda$, and $\hat{\mathcal{H}}_{MF}$ discussed above, there does not exist a widely accepted set of experimentally determined exchange constants for CoO Feygenson et al. (2011). Given the complexity of the mixed $j\rm{{}_{eff}}$ levels (Fig. 1(c)), we have previously investigated the pair response in the dilute monoxide Mg0.97Co0.03O, Sarte et al. (2018a) where chemical dilution removes both the magnetic order-induced molecular field and the accompanying structural distortion that are originally present in CoO Buyers et al. (1984); Furrer et al. (2011); Furrer and Waldmann (2013). A summary of the experimental results is presented in Fig. 2, taken from Ref. Sarte et al., 2018a. Utilizing probabilistic arguments, it was shown that the series of well-defined low energy magnetic excitations (Fig. 2(a)) present in Mg0.97Co0.03O correspond to excitations of Co2+ pairs. These pairs are described by the effective pair Hamiltonian given by

[TABLE]

corresponding to a 144 $\times$ 144 matrix in terms of the two-particle basis of $|l_{1}=1,m_{l_{1}},s_{1}=\frac{3}{2},m_{s_{1}}\rangle\otimes|l_{2}=1,m_{l_{2}},s_{2}=\frac{3}{2},m_{s_{2}}\rangle$, where $l_{i}$, $m_{l_{i}}$, $s_{i}$, and $m_{s_{i}}$ denote the eigenvalues of the $\hat{\mathbf{l}}$, $\hat{l}_{z}$, $\hat{\mathbf{S}}$, and $\hat{S}_{z}$ operators, respectively, for the $i^{th}$ particle. As summarized schematically in Fig. 2(b), the pair Hamiltonian $\hat{\mathcal{H}}_{pair}$ describes individual $j\rm{{}_{eff}}=\frac{1}{2}$ pair excitations as transitions between triplet ($\tilde{\Gamma}=1$) and singlet ($\tilde{\Gamma}=0$) levels separated by an energy $\Delta E=\tilde{\alpha}J$, where $\tilde{\alpha}$ behaves as an effective conversion factor between the energy transfer measured experimentally and the corresponding desired magnetic exchange constants. The solution to $\hat{\mathcal{H}}_{pair}$ as a function of exchange energy $J_{1,2}$ is shown in Fig. 2(c), with the solid black line and colored points corresponding to the exact solution to the above Hamiltonian and the measured energy positions, respectively. For comparison, the behavior predicted by the projection theorem of angular momentum is also presented. The deviation of the exact solution from the linear behavior predicted by the projection theorem is a consequence of the coupling of the ground state $j\rm{{}_{eff}}=\frac{1}{2}$ and higher energy $j\rm{{}_{eff}}={3\over 2}$ manifolds. Since the degree of coupling increases as $|J|\rightarrow|\lambda|$, the predicted values for exchange constants with larger magnitudes, specifically $J_{2}$, are particularly sensitive to the value of the spin-orbit coupling constant $\lambda$. This point will be addressed below in the context of the analysis of the single crystal data.

While the energy dependence affords estimates of the exchange constants, the relative distance $\mathbf{R}$ between the two Co2+ spins that participate in the exchange interaction, and thus the relative coordination shell (Fig. 2(a)) can be determined for each magnetic excitation from their momentum dependence via the first moment sum rule Hohenberg and Brinkman (1974); Sarte et al. (2018b); Haraldsen et al. (2005), as is summarized in Fig. 2. Given the ground state for antiferromagnetically/ferromagnetically coupled Co2+ ions is a triplet/singlet, the temperature dependence was used to establish the sign of the exchange constant Xu et al. (2000). A final summary of the estimates of the exchange constants for the first four coordination shells of CoO is presented in Table 1 below.

II.3 Mean-Field Theory for Multi-Level Spin-Orbit Excitons

As discussed above, the modeling of the neutron response requires an understanding of both the Hamiltonian and its commutator with the spin operators. In the previous section, we diagonalized the single-ion Hamiltonian $\hat{\mathcal{H}}_{1}$ such that $\hat{\mathcal{H}}_{1}|n\rangle=\omega_{n}|n\rangle$. Since all terms of the inter-ion Hamiltonian $\hat{\mathcal{H}}_{2}$ are based on the components of the spin operator $\hat{\mathbf{S}}$, these operators can be rotated onto the basis states of the single-ion Hamiltonian by use of the ladder operators that were previously defined in Eq. 9 with such a coordinate rotation being given by

[TABLE]

By writing the full Hamiltonian $\hat{\mathcal{H}}=\hat{\mathcal{H}}_{1}+\hat{\mathcal{H}}_{2}$ in terms of the ladder operators as defined in Eq. 12, and using the definition of the inter-level susceptibility $\hat{G}$ defined by

[TABLE]

where the indices $\alpha,\beta$ are either $+$, $-$, or $z$, the second term on the RHS of the equation-of-motion of the Green’s function (Eq. 4) reduces to three sets of commutators, termed diagonal, transverse, and longitudinal, with each involving spin operators rewritten in terms of ladder operators, as discussed in Section III of Supplementary Information sup . Buyers et al., Buyers et al. (1975) demonstrated that by combining the random phase decoupling method Wolff (1960); Cooke (1973); Yamada and Shimizu (1967, 1966) (Section IV of the Supplementary Information sup ) with the definitions of both the single-site response function

[TABLE]

and the Fourier transform of the exchange interaction $J(\mathbf{Q})$

[TABLE]

the full neutron response Fourier transformed into momentum $\mathbf{Q}$-space can be written as a set of coupled linear equations given by

[TABLE]

where $\omega$ has been relabeled as $E=\hbar\omega$. Here, we have employed the $T\rightarrow 0$ K single-site response function since the energy transfers under consideration in the current investigation ($\geq~{}20$ meV) are much larger than the sample temperature ($\sim$ 0.5 meV). The denominator of $g^{\alpha\beta}$ consists of $E_{no}=\hbar\omega_{n}-\hbar\omega_{0}$ corresponding to the energy associated with a transition from the ground state $|0\rangle$ to the $|n\rangle$ eigenstate of the single-ion Hamiltonian $\hat{\mathcal{H}}_{1}$, while the presence of the positive infinitesimal $\Delta$ is to ensure analyticality, and was set to 50% of the experimental resolution width (HWHM) on MERLIN Bewley et al. (2006) that was calculated by PyChop PyC . Coupling between the single-site response functions, and thus the dispersion of the total response functions $G^{\alpha\beta}$ is defined by $J(\mathbf{Q})$ which is parameterized by both $J_{ij}$ and $\mathbf{d}_{ij}$ denoting the exchange constant and displacement vector, respectively, between sites $i$ and $j$. As a first approximation, our calculations have considered the simplest case where the exchange interaction is spatially isotropic. We note that in general this is not case, owing to the anisotropy of the orbital configuration of Co2+.

By considering all possible combinations of indices $\alpha,\beta$ in Eq. 16 and noting that the non-zero single-site response functions for Co2+ in such a highly symmetric environment are: $g^{+-}$, $g^{-+}$, and $g^{zz}$, only three non-zero Green’s functions are obtained:

[TABLE]

with both $G^{++}$ and G*–* being both zero, as required by definition of the retarded Green’s function.

The simplest model for the long-range antiferromagnetic order in CoO is a type-II collinear antiferromagnetic magnetic structure Shull et al. (1951). Corresponding to (111) ferromagnetic sheets stacked antiferromagnetically along the [111] direction, this type of magnetic structure has been observed in CoO under pressure, despite the suppression of the structural distortion Ding et al. (2006). Such a model implies that CoO can be reduced to two unique magnetic sublattices; thus, the site indices $i$ and $j$ assume labels of either 1 or 2, and Eq. 17 becomes four coupled linear equations

[TABLE]

and

[TABLE]

with $J_{s}$ and $J_{d}$ denoting $J(\mathbf{Q})$ on the same ($i=j$) and different ($i\neq j$) sublattices, respectively. Solving these four coupled equations yields

[TABLE]

where $G^{-+}(\mathbf{Q},E)$ has the same form as $G^{+-}(\mathbf{Q},E)$ with indices $+$ $\longleftrightarrow$ $-$. The equations above are a function of the single-site response function $g^{\alpha\beta}$ and the Fourier transform of the exchange interaction $J(\mathbf{Q})$.

In contrast to the single-site response function, $J(\mathbf{Q})$ does explicitly depend on the particular magnetic sublattice under consideration, stemming from the presence of the position indices $i$ and $j$ in its definition given by Eq. 15. The Co2+ sites in the $d$- and $s$-sublattices for a particular coordination shell $m$ was determined by first selecting a reference Co2+ cation, thus defining a reference (111) plane as illustrated in Fig. 3. By definition of the type-II antiferromagnetic structure, Co2+ cations located on odd integer number of $(111)$ planes away from the reference plane are defined as belonging to the $d$ sublattice, while Co2+ cations located in the same or an even integer number of $(111)$ planes away are on the $s$ sublattice.

The program VESTA Momma and Izumi (2011) was used to determine the displacement vectors to calculate $J_{s,d}$ from Eq. 15 giving:

[TABLE]

and

[TABLE]

where the contributions from each coordination shell $m$ have been labeled explicitly.

By employing the definitions of the single-site response function $g^{\alpha\beta}$ (Eq. 14), and the Fourier transform of the exchange interaction $J_{s,d}(\mathbf{Q})$ (Eqs. 20 and 21), the total response function $G(\mathbf{Q},E)$ given by

[TABLE]

can be calculated numerically. In the T $\rightarrow$ 0 K limit, the imaginary part of $G(\mathbf{Q},E)$ is proportional to the dynamical structure factor (Eq. 2), and thus Eq. 1 may be reduced to

[TABLE]

demonstrating that the imaginary component of the total response function given by Eq. 22, with the inclusion of the magnetic form factor which here has been approximated by the isotropic magnetic form factor $f(Q)$, is directly proportional to the neutron magnetic cross section.

II.4 Parameters: Initial Values & Orbital Configurations

Having presented our model, we now discuss the parameterization of the spin-orbit excitations in CoO. Since our model approximates CoO as a tetragonally distorted type-II antiferromagnet, the single-site response $g^{\alpha\beta}$, itself being a function of the single-ion Hamiltonian $\hat{\mathcal{H}}_{1}$ (Eq. 6), is defined by three parameters: $\lambda$, $\Gamma$, and $H_{MF}$. The spin-orbit coupling parameter $\lambda$ was taken to be $-16$ meV, corresponding to its value reported by Cowley et al. Cowley et al. (2013). An initial estimate for the mean molecular field $H_{MF}$ was determined by first extracting the value for $\sum\limits_{i>j}z_{ij}J_{ij}$ from the reported Kanamori (1957); Singer (1956); Nagamiya et al. (1955) Curie-Weiss temperature $\theta{\rm{{}_{CW}}}=-330$ K ($-28.4$ meV) via its mean field definition

[TABLE]

where $\zeta$ is a scale factor of 1.9 calculated by Kanamori Kanamori (1957) accounting for mixing between the ${}^{4}F$ and ${}^{4}P$ free-ion states. Inserting the value of $\sum\limits_{i>j}z_{ij}J_{ij}$ into the definition of $H_{MF}$ given by Eq. 5, yields an initial estimate of 64.8 meV. An initial estimate of the tetragonal distortion parameter $\Gamma$=$-8.76$ meV was determined by scaling the value of $-1.49$ meV that was reported for KCoF3 Buyers et al. (1971) by an empirical factor of 0.0116/0.00197=5.89 corresponding to the ratio of their respective tetragonal distortions $\delta a/a$.

To define $J(\mathbf{Q})$, we have taken the values for the exchange constants for the dilute monoxide Mg0.97Co0.03O Sarte et al. (2018a) as estimates for pure undiluted CoO since these exchange constants correspond to a Curie-Weiss temperature (Eq. 23) in close agreement with the value reported for CoO. However, our investigation on Mg0.97Co0.03O also revealed that each coordination shell possessed the possibility for both antiferromagnetic and ferromagnetic coupling, with the exception of the second nearest neighbor which is fixed to be antiferromagnetic by the 180∘ Co2+-O2--Co2+ superexchange pathway. As illustrated in Fig. 2(c,inset), such dual behavior is a direct consequence of the $t_{2g}$ degeneracy of the high spin $d^{7}$ configuration of Co2+, and thus a particular choice of $J$, be it antiferromagnetic or ferromagnetic, corresponds to a specific local orbital configuration. By incorporating this dual behavior for coordination shells 1, 3, and 4, we must consider all $2^{3}$=8 sets of exchange constants of the form $xAxx$, where $x$ can be either antiferromagnetic (A) or ferromagnetic (F). Furthermore, since our model incorporates the effects of a tetragonal (uniaxial) distortion, we must also distinguish the involvement of a distorted or undistorted bonding configuration for each of the 8 $xAxx$ orbital configurations. Thus, with these 2 additional degrees of freedom, our model must consider 16 different orbital configurations of the form $xAxx\gamma$, where the index $\gamma$ =1 or 2, distinguishing the presence or absence of distorted bonding configurations.

Each of these 16 $xAxx\gamma$ combinations may be interpreted as a unique orbital configuration, physically corresponding to a unique type of “domain” in the bulk CoO single crystal, each of which is subject to a different mean molecular field $H_{MF}$. In contrast, since all cations under consideration are assumed to be Co2+ in octahedral coordination subject to a cooperative Jahn-Teller tetragonal distortion, $\lambda$, $\Gamma$, and the individual $J_{n}$ values (where $n$ denotes a particular type of coupling in a coordination shell $m$) are fixed to be equal for each of the 16 $xAxx\gamma$ orbital configurations. By noting that our neutron spectroscopic measurements were performed with a large experimental beam that irradiated a macroscopic number of domains in the single crystal of CoO, our model considers the mean contribution from all 16 equally weighted $xxAx\gamma$ orbital configurations. The initial parameters considered in the model and their initial values are summarized in Tab. 1.

III Experiment

Having discussed both the underlying theory of Co2+ magnetism and the corresponding physical parameters that constitute our spin-orbit exciton model, we will now address the experimental results from neutron scattering experiments on a single crystal of CoO. This section begins with a description of the experimental techniques, followed by a summary of the neutron spectroscopic data. We conclude this section with a description of how our model was used to interpret the low-energy fluctuations of CoO deep within the Néel regime.

III.1 Experimental Details

Sample Preparation: Polycrystalline samples of CoO were synthesized by annealing high purity Co3O4 ($>$ 99.99%) under flowing Ar at 1200oC for 36 h with intermittent grinding until laboratory x-ray diffraction confirmed the absence of the Co3O4 precursor. The phase pure CoO powder was compressed into cylindrical rods using a hydraulic press and subsequently annealed under flowing Ar at 1275oC for 24 h in a horizontal annealing furnace. Crystal growth was performed using the floating zone technique with a four-mirror optical floating zone halogen furnace (CSI system Inc.), yielding a 10 g single crystal of CoO ($l$ = 50 mm, $\phi$ = 8 mm). The feed and seed rods were counter-rotated at 35 rpm with a vertical translation of 2 to 4 mm hr*-1* in a pure Ar atmosphere. The initial polycrystalline seed rod was replaced for subsequent runs by single crystal seeds from earlier growths. Previous Cowley et al. (2013) optical and scanning electron microscopy, x-ray diffraction and DC susceptibility measurements on the single crystal confirmed the presence of a single growth domain with a mosaic spread of approximately 0.1o and the absence of both multiple magnetic or strain domains on the crystal surface and Co3O4 impurities, respectively.

Neutron Inelastic Scattering: Neutron inelastic scattering measurements were performed on the MERLIN direct geometry chopper spectrometer Bewley et al. (2006) at the ISIS neutron spallation source (Didcot, U.K.). The $t_{0}$ chopper was spun at 50 Hz in parallel with the “sloppy” Fermi chopper package to fix the incident energy with the energy transfer defined as $E=E_{i}-E_{f}$. To access a large dynamic range, three fixed incident energies $E_{i}$ of 110 meV, 75 meV, and 45 meV were selected with Fermi chopper frequencies of 350, 300, and 250 Hz, providing a resolution at the elastic line (E=0) of 7.3, 4.8, and 2.7 meV, respectively. A 5 g portion of the CoO single crystal was mounted in a top loading closed cycle refrigerator such that the [110] and [001] crystallographic axes lay within the horizontal plane. A tomographic reconstruction in momentum-space was accomplished by rotating the crystal about the [0$\overline{1}$0] axis over 120o in 0.5o steps.

The four dimensional $(\mathbf{Q},E)$ experimental data at each angle $\Psi$ and $E_{i}$ was collected at 5 K for an accumulated charge of 30 $\rm{\mu A\cdot hr}$ on the spallation target. The raw experimental data was normalized by accumulated proton charge, corrected for detector-efficiency using a vanadium reference sample, and reduced by the Mantid data analysis software.Arnold et al. (2014); Taylor et al. (2012) Visualization and manipulation of reduced experimental data including rebinning and projections were performed using the HORACE software package distributed by ISIS Ewings et al. (2016).

III.2 Experimental Data

We begin by first presenting a summary of the experimental data from single crystal neutron spectroscopic measurements allowing for a direct comparison with previous work on CoO to establish consistency. The experimental data from the MERLIN chopper spectrometer at 5 K is presented in Fig. 4 in the form of $(\mathbf{Q},E)$ slices along both $(1.5\pm 0.1,1.5\pm 0.1,L)$ (Figs. 4(a,c)) and $(2.0\pm 0.1,2.0\pm 0.1,L)$ (Fig. 4(b)) capturing both the magnetic zone center and boundary, where both are compared to previously published work in the form of $\mathbf{Q}$-integrated cuts presented in Figs. 4(d) and (e), respectively. As illustrated in Fig. 4(a), a $(\mathbf{Q},E)$ slice with an incident energy $E_{i}$=110 meV exhibits a band of excitations extending from $\sim$ 20 meV up to $\sim$ 60 meV energy transfer, corresponding to a similar range in energy transfer reported by previous THz Satoh et al. (2017) and Raman Hayes and Perry (1974) spectroscopic measurements.  These excitations decrease in intensity with increasing $L$, as is expected for the Co2+ magnetic form factor, thus indicating these excitations are possibly magnetic. A higher resolution slice with an $E_{i}$= 45 meV (Fig. 4(c)) reveals that the band of excitations corresponds to a fine structure consisting of a series of modes that are unevenly spaced in energy, in agreement with previous triple-axis, Yamani et al. (2010); Sakurai et al. (1968) time-of-flight data, Tomiyasu and Itoh (2006) and Raman spectroscopy, Struzhkin et al. (1993) with the exception of a broader peak reported for triple-axis measurements at $\sim 40$ meV. However, it is important to note that the previously reported triple-axis measurements employed final energies $E_{f}$=14.6 meV and 30.5 meV, both of which potentially produce spurious signals near 40 meV as a result of weak elastic scattering corresponding to $E_{i}$ $\rightarrow$ 4 $E_{f}$ ($\lambda_{f}/2$) and 4 $E_{i}$ ($\lambda_{i}/2$) $\rightarrow$ 9 $E_{f}$ ($\lambda_{f}/3)$ Shirane et al. (2004). This weak elastic process may have contributed to the extra scattering intensity that was observed in the previously reported triple-axis data, yet is absent in our current time-of-flight data. We also note the presence of nearly dispersionless optical phonon branches near this energy may also contribute to the overall neutron cross section Wdowik and Parlinski (2007); Goyal and Goyal (1977). We will later discuss how this fine structure near the magnetic zone center can be understood in terms of the spin-orbit excitations.

Slices and $\mathbf{Q}$-integrated cuts through the magnetic zone boundary presented in Figs. 4(b) and (e) shows additional complexity present in the neutron response. At low energy transfers below $\sim$20 meV, strong acoustic phonons can be seen in Fig. 4(b) to disperse from the even integer positions in momentum transfers. The phonon nature of these excitations is confirmed by the fact that the cross section grows with increasing momentum transfer $L$, in contrast to the magnetic cross section that is subject to the Co2+ form factor. In addition to the acoustic phonons, a flat band is observed at $E\sim$ 20 meV. This band also exhibits higher intensity at large momentum transfers, indicating lattice fluctuations as the origin. However, it should be noted that a peak in the aluminum phonon density of states exists near this energy transfer, thus suggesting that this particular band likely corresponds to scattering from the sample can and walls of the cryostat. At higher energy transfers, two distinct bands of excitations are present over a large range of energy transfers spanning from $\sim$30 to $\sim$70 meV. The intensity of these bands weaken with momentum transfer $L$; however, as illustrated in Fig. 4(e), these excitations overlap with lattice vibrational modes previously identified by infrared spectroscopy Austin and Garbett (1970). We will later discuss the origin of these excitations in terms of magneto-vibrational scattering Egelstaff (1965) by comparing the dispersion of these two high energy bands at the magnetic zone boundary to those of phonons measured at large momentum transfers.

III.3 Comparison between Experimental Data and the Spin-Orbit Exciton Calculation

Having summarized our experimental data, we now present a direct comparison to our calculated parameterization based on the Green’s function approach that has been discussed above. Here, we first present the final model used to describe the experimental data. This is followed by a discussion concerning how such a model and its refined parameters were obtained.

As illustrated in Fig. 5, by allowing the value of mean molecular field $H_{MF}$ to refine independently for each of the equally weighted 16 $xAxx\gamma$ domains, each possessing identical refined values of $\lambda$, $J$, and $\Gamma$, our mean-field multi-level spin-orbit exciton model successfully reproduces both the fine structure at the magnetic zone center and the broad excitations at the $(1.5,1.5,-1)$ zone boundary, while capturing the steeply dispersive columns of scattering observed at higher energy transfers (Fig. 6). The refined values of the mean molecular field $H_{MF}$ are listed in Tab. 2 for the 16 $xAxx\gamma$ orbital configurations, each with refined values of $\lambda$, $\Gamma$, and $J_{m,\xi}$ summarized in Tab. 3.

The need for all 16 $xAxx\gamma$ domains is illustrated in Figs. 7(e) and (f), where despite the success of the spin-orbit exciton model in reproducing the bandwidth of the excitations between $\sim$20 meV up to $\sim$60 meV with the initial parameters listed in Tab. 1, the use of a single $xAxx\gamma$ orbital configuration fails to reproduce the low energy fine structure at the zone center and the steeply dispersive columns of high energy scattering at the zone boundary, instead predicting the presence of a single dominant highly dispersive $G^{-+}$ mode (Figs. 7(b) and (g)). While both $G^{++}$ and $G^{--}$ modes exhibit negligible intensity (Fig. 7(d)), a direct consequence of the definition of $G^{\alpha\beta}$ given by Eq. 16 above, the spin-orbit exciton model does predict two additional gapped modes corresponding to the longitudinal $G^{zz}$ mode and the transverse $G^{+-}$ mode. As summarized in Figs. 7(a,c), it is clear that both the weakly dispersive $G^{+-}$ and Gzz modes cannot account for the missing spectral weight in the fine structure, as both modes are significantly weaker in intensity relative to the dominant $G^{+-}$ mode, while contributing solely at higher energy transfers ($\gtrsim$40 meV). As illustrated in Figs. 6 and 8, while the longitudinal $G^{zz}$ corresponds to a band of scattering that is broad in both momentum and energy while being located above the dominant $G^{-+}$ component, the transverse $G^{+-}$ corresponds to a flat band of scattering centered about $\sim$40 meV.

In addition to the necessity for 16 $xAxx\gamma$ domains, the discrepancy between experiment and the calculated dispersion for each of the $G^{\alpha\beta}$ components that were presented in Fig. 7 confirmed the need for the optimization of the spin-orbit exciton model’s parameters. A summary of the influences of each of the model’s parameters for a fixed orbital configuration: $J$ via $\lambda$, $H_{MF}$, and $\Gamma$ on the calculated spectra is presented in Fig. 9. Calculations are shown by false colormaps and compared against constant momentum cuts at the $\bf{Q}$=$(1.5,1.5,-1.0)$ magnetic zone boundary and $\bf{Q}$=$(1.5,1.5,-0.5)$ zone center. For illustrative purposes, two different orbital configurations $AAAA$ and $AAAF$ are shown being denoted by solid and dashed lines, respectively.

As was previously noted in the discussion of the inter-ion Hamiltonian $\hat{\mathcal{H}}_{2}$, the conversion factor $\tilde{\alpha}$ (Fig. 2) between the energy transfers measured by neutron spectroscopy and the magnitude of the corresponding magnetic exchange constants were determined by diagonalizing the pair Hamiltonian $\hat{\mathcal{H}}_{pair}$ (Eq. 11). Since deviations away from the linear dependence predicted by the projection theorem occurs when the value of $|J|\rightarrow|\lambda|$, the value of larger magnetic exchange constants, such as the strong 180*∘* antiferromagnetic superexchange $J_{2}$, is particularly sensitive to $|\lambda|$. As illustrated in Fig. 9(a), the strong sensitivity of $|J_{2}|$ on the value of $|\lambda|$ is of particular concern for CoO since the experimentally determined value of $\lambda=-16(3)$ meV possesses a significant relative error of almost 20% Cowley et al. (2013), corresponding to a large range of possible exchange values.

The influence of the large relative error for $\lambda$ is summarized in Fig. 9(b), illustrating that an increase of $|\lambda|$ by 20% to $-19$ meV from its initial value of $-16$ meV for a fixed value of a given molecular field $H_{MF}$ results in a significant shift to higher energy transfers for the dominant $G^{-+}$ component at both the magnetic zone center and boundary. This increase in energy transfer is much more prominent for the $G^{-+}$ component compared to the less intense $G^{+-}$, with the greatest increase for both components occuring at the magnetic zone center. Similar behavior in the calculated response is observed for a change in the mean molecular field $H_{MF}$, as illustrated in Fig. 9(c), with shifts in energy transfers significantly larger compared to the same relative change in the value of $|\lambda|$. In contrast with the three parameters presented so far, the influence of the tetragonal distortion parameter $\Gamma$ is most pronounced on the less intense weakly dispersive $G^{+-}$ mode. As illustrated in Fig. 9(d), $\Gamma$ provides a mechanism to shift the energy transfers of the $G^{+-}$ mode without inducing a comparable shift for the dominant $G^{-+}$ mode for a fixed set of values for $H_{MF}$ and $J$ via a fixed value for $\lambda$.

With such a large number of domains under consideration, constraints on the parameter space for $J$ via $\lambda$, $\Gamma$, and $H_{MF}$ were required to ensure convergence for a least squares optimization. As summarized in Tab. 3, the model presented in Fig. 5 allowed $\lambda$ to vary from $-13$ to $-19$ meV, corresponding to the experimental error associated with its reported value Cowley et al. (2013), while the value of $\Gamma$ varied from $-8.76$ to $-6.16$ meV, corresponding to the range of its values reported in the literature Cowley et al. (2013). The exchange constants $J_{m,\xi}$ for each coordination shell $m$ = 1 …4, and type $\xi$=$A$ and $F$ were allowed to vary $\pm$20% from their starting values obtained from the excitation energies derived from Mg0.97Co0.03O Sarte et al. (2018a) with a Hamiltonian $\hat{\mathcal{H}}_{pair}$ employing the value of $\lambda$ under consideration. Such a deviation of the exchange constants was rationalized by noting that the estimates of the exchange constants listed in Tab. 1 are based on the superexchange pathways present in Mg0.97Co0.03O. Due to the prevalence of the non-magnetic Mg2+ in the dilute monoxide, these pathways and their superexchange constants most likely differ compared to those present in pure CoO. As a first approximation, the relative deviation of all magnetic exchange constants $J$ from their respective values in Mg0.97Co0.03O were set to be equal for all coordination shells $m$ and type $\xi$. As previously discussed in Section II.4, the values of $\lambda$, $J$, and $\Gamma$ were constrained to be equal for each of the 16 domains, while no such restriction was applied to the mean molecular field $H_{MF}$ which was itself allowed to vary independently for each domain from a value of 0 to an arbitrarily large upper limit. In our model, this limit was set to 100 meV, corresponding to $\sim$1.5 times the value of the initial value of 64.8 meV.

Among all the spin-orbit exciton model’s parameters, the expansion of the parameter space for the exchange constants $J_{m,\xi}$ proved to be of particular importance for the model’s success in reproducing both the fine structure at the magnetic zone center and the broad excitations at the $(1.5,1.5,-1)$ zone boundary simultaneously. As summarized in Section VI of the Supplementary Information sup , the restriction of the exchange constants $J_{m,\xi}$ to be equal to their respective values originally reported for Mg0.97Co0.03O Sarte et al. (2018a) for a given value of $\lambda$, resulted in failure of the spin-orbit exciton model to reproduce both broad excitations at the $(1.5,1.5,-1)$ zone boundary (Fig. S2). As summarized by Tab. SIII, such a restriction placed on the values of exchange constants resulted in the value of $\lambda$ being refined to its most negative permissible value. Such behavior is a reflection of the model’s attempts to minimize the average value for $H_{MF}$ in order to capture the intensity at lower energy transfers at the $(1.5,1.5,L)$ zone boundary. As illustrated in Fig. S3, the clear failure of the model, even by expanding the parameter space of $\lambda$ to include all values down to $-23$ meV, confirmed the inability of the model to reproduce the data at the zone boundary while employing the exact exchange constants measured in Mg0.97Co0.03O.

III.4 Magneto-Vibrational Scattering and the Magnetic Zone Boundary

Despite the success of the spin-orbit exciton parameterization of the experimental data along $(1.5,1.5,L)$, the model still fails at the $(2,2,L)$ zone boundary, suggesting the presence of additional physics that is beyond our localized model. The presence of such prominent magnetic scattering near $(2,2,L)$ is particularly unusual due to the predictions from both our spin-orbit exciton model and linear spin-wave theory for nearly zero intensity for magnetic fluctuations at these points in reciprocal space for an antiferromagnetic structure.

By comparing the scattering along $(2,2,L)$ to phonon branches near different nuclear zone boundaries, we note that there is a distinct overlap in dispersion. As is summarized in Fig. 10(a), a strong and steeply dispersing optical phonon mode near ${\bf{Q}}=(5.5,3.5,4.5)$ exhibits an identical dispersion to that of the mode centered about $(2,2,0)$, a mode that is not accounted for by our localized model. The phonon nature of the gapped mode in Fig. 10(a) is confirmed by the fact that the intensity increases with ${\bf{Q}}$.  A similar observation is illustrated in Fig. 10(c) for energy transfers greater than 60 meV, where a gapped phonon mode, identified by both its ${\bf{Q}}$ dependence and previously reported first principles calculations Wdowik and Parlinski (2007), exhibits a dispersion that is identical to the magnetic scattering around $(2,2,0)$.

The apparent similarity between the dispersion of phonons at high $\mathbf{Q}$ and the magnetic scattering that was originally unaccounted for in our model may suggest that our data at the zone boundaries includes a magneto-vibrational contribution to the neutron cross section. Corresponding to an indirect energy exchange between the neutron and a vibrating nucleus via the electro-magnetic interaction between the neutron and the unpaired electrons of the nucleus, the magneto-vibrational neutron cross section Egelstaff (1965) is identical to the more commonly employed direct one-phonon cross section employing the nuclear force, with the exception that the nuclear scattering length $b$ is replaced by the magnetic scattering length given by

[TABLE]

where $\alpha$ is the angle between the momentum transfer ${\bf{Q}}$ and the ordered magnetic moment direction ${\bf{\mu}}$. The presence of the isotropic magnetic form factor $f(Q)$ guarantees that the cross section will ultimately decay with momentum transfer, regardless of the underlying phonon origin of the scattering. This particular cross section requires an ordered magnetic structure and has been used previously to characterize the dynamic magnetic form factor Brown et al. (1988). Similar cross sections have been reported in the rare earth magnetic pyrochlores Fennell et al. (2014) and doped manganites Hennion et al. (2002). We note that in the case of CoO, the magnetic ions are strongly correlated and coupled in all three dimensions, implying that that the correlations are maintained at large energy transfers in the range of 30-60 meV where prominent phonon modes exist. Based on the comparison presented above, we suggest that both the excitations along $(2,2,L)$, and other excitations that were unaccounted for by our model such as the weak magnetic scattering at $\sim$40 meV at the $(1.5,1.5,L)$ zone boundary all originate from phonon modes, rather than the underlying magnetic Hamiltonian $\hat{\mathcal{H}}$.

While we have cast this discussion in terms of the magneto-vibrational cross section, which itself does not provide any new information on the underlying Hamiltonian, these excitations may be indicative of a coupling between the lattice and magnetic degrees of freedom. Such magneto-elastic coupling in the underlying Hamiltonian has been investigated recently in CeAuAl3 Čermák et al. (2019). It is particularly compelling to consider such a claim given that these modes are located near the expected energy scale for the single-ion spin-orbit transitions from $|j{\rm{{}_{eff}}}=\frac{1}{2}\rangle\rightarrow|j{\rm{{}_{eff}}}=\frac{3}{2}\rangle$. We will not investigate this coupling further given the need for first principles calculations to reconcile past reported phonon data and the modes identified in our current study.

III.5 The High-Energy Response and the Co2+ Form Factor

We have demonstrated that the experimental data at low energy transfers can be succesfully parameterized in terms of spin-orbit excitons with crystal field and exchange parameters based on our previous work reported on MgO substituted with Co2+. In Section III.3, such a mean-field multi-level spin-orbit exciton model was shown to successfully reproduce the data in pure CoO near the zone center for energy transfers below 40 meV, while the failure of the model near the magnetic zone boundaries up to energy transfers of $\sim$60 meV was attributed in Section III.4 to either magneto-vibrational scattering or a coupling to underlying phonon excitations that were identified at large momentum transfers. In this final subsection, we will address the magnetic excitations at large energy transfers above $\sim$ 65 meV, where no phonon scattering is expected, exceeding the dynamic range predicted by first principle calculations Wdowik and Parlinski (2007) and measured by both our time-of-flight experiment and previous experiments Kant et al. (2008); Sakurai et al. (1968); Goyal and Goyal (1977).

In Fig. 11 we present a comparison between the scattering for high energy transfers ([70, 105] meV) with the lower energy magnetic fluctuations that were discussed above in terms of our spin-orbit exciton parameterization. As summarized in Fig. 11(a), excitations located at low momentum transfers extend up to high energies. Previously identified in Fig. 6 as being magnetic in origin, it appeared that these excitations were successfully reproduced by the spin-orbit exciton model. These “columns” of scattering were interpreted as a result of the overlap at the magnetic zone boundary of multiple $G^{-+}$ components from different local orbital arrangements that were required to capture the fine structure of the strongly dispersing fluctuations (Fig. 11(b)) centered about the zone centers and found at lower energy transfers. While the analysis presented so far suggest that the $\mathbf{Q}$-dependence for all magnetic excitations appear to follow the Co2+ form factor, the energy-integrated slice in Fig. 11(a) demonstrates this is not the case for high energy transfers, where the intensity decays more rapidly than $f(Q)$. Such an observation can be confirmed by comparing momentum cuts of the magnetic fluctuations at small momentum transfers with the magnetic form factor for Co2+. In contrast with the lowest energy transfers ([20-40] meV, Fig. 11(e)) where the magnetic fluctuations follow the form factor, deviations from such localized magnetic behavior begin to appear at $(1.5,1.5,\pm 0.5)$ for intermediate energy transfers ([40,70] meV, Fig. 11(d)), with the deviations being particularly prominent around $(1.5,1.5,0)$ at the highest energy transfers ([70,105] meV, Fig. 11(e)). The rapid decay of intensity with momentum transfer maybe indicative of a real space object that is extended spatially, and not due to localized magnetism.

We now speculate as to the possible origin for such delocalized magnetism. In the case of the highest energy transfers in CoO, where the relative deviation from the form factor is greatest, the scattering is steeply “dispersive” in energy, indicative of a large underlying energy scale. These fluctuations are highly reminiscent of the magnetic response found in itinerant magnets such as CeRhIn5 Stock et al. (2015), the cuprates Stock et al. (2007, 2010), and also iron-based systems Stock et al. (2014); Plumb et al. (2018). In the case of CeRhIn5, the high-energy steeply dispersive excitations were found to be longitudinally polarized and occupied a region in $(\mathbf{Q},E)$ phase space where two magnon processes were kinematically allowed. Termed the “1+2” model, multi-magnon decay processes were used to provide a heuristic description of the data in the case of CeRhIn5. Motivated by the qualitative similarities in the excitation spectrum between CeRhIn5 and CoO, we investigated the possibility that the steeply dispersing excitations observed in CoO overlap a similarly allowed region of $(\mathbf{Q},E)$ phase space. The phase space permitted for such a decay of the low energy magnetic fluctuations was calculated using a simple model based on energy and momentum conservation that is given by Huberman et al. (2005); Stock et al. (2018); Songvilay et al. (2018)

[TABLE]

where $E_{\mathbf{Q}_{1,2}}$ are the energies of transverse excitations at a given momentum transfer. As shown in Fig. 12, the kinematically allowed region overlaps in both momentum and energy with the steeply dispersing excitations at the highest energy transfers; a region where the spin-orbit exciton model predicts the presence of a longitudinally polarized excitations ($i.e.$ $G^{zz}$ mode), as was previously illustrated in Fig. 6(b). Employing the observed overlap that is summarized in Fig. 12, we speculate that the spectral weight for these steeply dispersive high energy excitations that are localized in momentum draws from this longitudinally polarized mode, analogous to what was observed in previously investigated itinerant magnets. Magnetism spatially extending beyond the Co2+ site has been suggested theoretically with some moment expected to be presented on the oxygen atom van Elp et al. (1991). Given the extended nature of the magnetism in such a scenario, the magnetism would be expected to decay faster in momentum transfer than the isotropic Co2+ form factor. An intermediate example has been reported in the case of Sr2CuO3 Walters et al. (2009) which found excellent agreement between the magnetic form factor and a model including strong covalent bonding and hybridization of the $3d$ orbitals. It could be that the higher energy excitations are more sensitive to such a situation in CoO.

IV Concluding Remarks

We have presented a neutron spectroscopic study of the magnetic fluctuations in the Mott insulator CoO. We have parameterized the low energy magnetic excitations near the magnetic zone center in terms of a mean-field multi-level spin-orbit exciton model incorporating multiple structural domains owing to multiple local orbital arrangements, as well as a prominent tetragonal Jahn-Teller distortion. Dispersive excitations at the zone boundaries mimicking magnetic scattering at low $\mathbf{Q}$ that were originally unaccounted for by the spin-orbit exciton model were found to exhibit similar dispersions as phonons measured at larger momentum transfers, suggesting that the model’s failures at the zone boundaries may be attributed to magneto-vibrational scattering or possibly coupling to lattice degrees of freedom. Finally, we report a discrepancy between the $\mathbf{Q}$-dependence of excitations at high energy transfers and the behavior predicted by the Co2+ form factor. Despite the strong insulating nature of CoO, we speculate that such a discrepancy corresponds to a breakdown of spin-orbit excitons that may be accompanied by a crossover from localized to spatially-extended magnetism, reminisicent of an itinerant-like response or strong covalency.

Acknowledgements.

We would like to convey our sincere gratitude to both W. J. L. Buyers and the late R. A. Cowley for their assistance throughout the preliminary stages of this investigation. We acknowledge fruitful conversations with C. R. Wiebe, G. M. McNally, J. A. M. Paddison, G. Perversi, S. E. Maytham, and K. J. Camacho. We are grateful to the Royal Society, the STFC, the ERC, the EPSRC, and the Carnegie Trust for the Universities of Scotland for financial support. P. M. S. acknowledges additional financial support from the University of California, Santa Barbara through the Elings Fellowship, the CCSF, RSC, ERC, and the University of Edinburgh through the GRS and PCDS. A portion of this work was supported by the DOE, Office of Science, Basic Energy Sciences under Award DE-SC0017752.