Magnetic Fluctuations and the Spin-Orbit Interaction in Mott Insulating CoO

Motivated by the presence of an unquenched orbital angular momentum in CoO, a team at Chalk River, including a recently hired research officer Roger Cowley, performed the first inelastic neutron scattering experiments on the classic Mott insulator [Sakurai $\textit{et al.}$ 1968 Phys. Rev. $\mathbf{167}$ 510]. Despite identifying magnon modes at the zone boundary, the team was unable to parameterise the low energy magnetic excitation spectrum below $T\rm{_{N}}$ using conventional pseudo-bosonic approaches. It would not be for another 40 years that Roger, now at Oxford and motivated by the discovery of the high-$T_{c}$ cuprate superconductors [Bednorz&Muller 1986 Z. Phys. B $\mathbf{64}$ 189], would make another attempt at the parameterisation of the magnetic excitation spectrum that had previously alluded him. Upon his return to CoO, Roger found a system embroiled in controversy, with some of its most fundamental parameters still remaining undetermined. Faced with such a formidable task, Roger performed a series of inelastic neutron scattering experiments in the early 2010's on both CoO and a magnetically dilute structural analogue MgO. These experiments would prove instrumental in the determination of both single-ion [Cowley $\textit{et al.}$ 2013 Phys. Rev. B $\mathbf{88}$ 205117] and cooperative magnetic parameters [Sarte $\textit{et al.}$ 2018 Phys. Rev. B $\mathbf{98}$ 024415] for CoO. Both these sets of parameters would eventually be used in a spin-orbit exciton model [Sarte $\textit{et al.}$ 2019 Phys. Rev. B $\mathbf{100}$ 075143], developed by his longtime friend and collaborator Bill Buyers, to successfully parameterise the complex spectrum that both measured at Chalk River almost 50 years prior. The story of CoO is of one that has come full circle, one filled with both spectacular failures and intermittent, yet profound, little victories.


Introduction
In the years following Brockhouse's invention of the neutron triple axis spectrometer [1][2][3], Chalk River was considered a premier neutron scattering facility, attracting some of the most brilliant scientific minds of the 20 th century [4][5][6][7]. Many of whom would go on to have long and distinguished careers, whose profound influences are still felt throughout modern physics today [8]. One such scientist was Roger A. Cowley, whose 6 year tenure (1964)(1965)(1966)(1967)(1968)(1969)(1970) at Chalk River would prove to be particularly influential on the newly minted PhD. A period in his career that would see both high productivity and the establishment of lifelong collaborations. ∆, difficulties arose due to limitations inherent to optical techniques [45,59,[83][84][85][86][87]. As will be described later, significant entanglement of the spin-orbit split manifolds due to the presence of multiple far-reaching large exchange interactions with values comparable to that of spinorbit coupling prevented the extraction of both the spin-orbit coupling constant λ and exchange constants J using conventional approaches [16,17,22,[31][32][33]88]. In the face of a seemingly Herculean task, little progress was seen for almost another 20 years [27].
In a sense, Roger's renewed interest in CoO was ultimately driven by the same motivation. He believed that CoO, with its unquenched orbital angular momentum, could be used as a tool in the establishment of a complete and comprehensive theoretical framework to describe the electronic structure of the cuprates and ultimately their low temperature properties [108][109][110][111]. Troubled by the ever increasing complexity of the Hamiltonians employed in conventional linear spin wave theory approaches to describe such systems [112][113][114][115][116][117][118], Roger desired a framework based on a minimalist Hamiltonian given bŷ describing explicitly both single-ion and cooperative magnetism via Heisenberg exchange. Despite recent inelastic neutron scattering experiments reported by Tomiyasu & Itoh [31] and Yamani et al. [32,33] revealing a magnetic excitation spectrum that was much more complex than what was first measured by the team at Chalk River ( Fig. 1(b)) -to be later confirmed by his own measurements [119] (Figs. 1(c)-(f)) -Roger still believed it was possible that the low energy magnetic fluctuations of a system that had alluded him decades prior could be fully parameterised utilising such a minimalist approach. From Roger's perspective, the use of such a Hamiltonian would provide a much clearer understanding of the complex low temperature magnetism in systems such as CoO by permitting a direct input of spin-orbit, anisotropy, and magnetic exchange parameters, instead of relying on indirect contributions based on expressions derived from higher order perturbation theory [120][121][122][123].

Single-Ion Physics of CoO
Despite a myriad of both experimental and theoretical studies spanning decades and a rejuvenated interest in CoO spurred on by the cuprates, Roger was still faced with the fact that both terms comprisingĤ in Eq. 1 were still not well-understood at the turn of the new millennium [27,31,37,58,83,101,124]. It was at this time that he realised that in order for his vision to come to fruition, the single-ion physics of the 3d 7 Co 2+ in CoO would need to be first thoroughly understood and ultimately parameterised, before addressing the additional complications arising from cooperative magnetism [27,[31][32][33]120,[125][126][127]. Such a reinvestigation of the single-ion physics would be accomplished experimentally with the technique of neutron scattering [111]. In contrast to optical spectroscopy corresponding to the technique of choice at the time, neutron scattering is not limited by the selection rules of ∆l = ±1 [84,128], reducing the difficulty in measuring the often subtle dd transitions (∆l = 0) [85,86].
In this section, we describe how Cowley et al. [111] employed neutron scattering to parameterise the two most dominant contributions toĤ CF , thereby laying the foundation for the significant progress that would be seen throughout the next decade in the understanding the low energy magnetic excitations in CoO.
In the form presented in Eq. 1,Ĥ is valid only in the paramagnetic regime [129,130]. In the case of T < T N , it can be shown [131,132] that the introduction of a contributionĤ M F from the mean molecular field stemming from long range magnetic order that is given bŷ where H M F (i) is the molecular field constant for site i that is defined as to Eq. 1 allowsĤ to be rewritten as a sum of a single-ion (Ĥ 1 ) and an inter-ion (Ĥ 2 ) term given byĤ respectively, where Ŝ z (j) denotes a thermal average given by Ŝ η = n f n n|Ŝ η |n , weighted by the Boltzmann thermal population factor f n . The importance of this particular separation ofĤ into two terms will later become apparent when discussing the spin-orbit exciton model [119,131]. The expression for the single-ion Hamiltonian of CoO can then be rewritten aŝ consisting of the individual contributions from the crystalline electric field, spin-orbit coupling, structural distortion, and mean molecular field, respectively [22,31,88,111,119,129,130].

Crystalline electric field, Tanabe-Sugano & determination of 10Dq
At the time of Roger's neutron scattering experiments, the crystal field scheme and the formalism to addressĤ CEF in CoO was still controversial, stemming from the limitations of optical measurements and from the question concerning which theoretical framework was appropriate, respectively [45,59,83,[133][134][135][136][137][138]. At one extreme was weak crystal field theory, where the crystal field splitting of the free-ion states was much smaller than the energy difference between these free-ion states, while at the other extreme was strong crystal field theory, where the opposite was true [15]. For the former case,Ĥ CEF is treated as a perturbation to the free-ion states, with a complete set of commuting observables of:L 2 ,L z ,Ŝ 2 andŜ z with corresponding good quantum numbers L, m L , s and m s . For the case of octahedral coordination [139,140], H CEF can be expressed in terms of the Stevens operatorsÔ 0 4 andÔ 4 4 , and the numerical coefficient In contrast, the basis in strong crystal field theory corresponds to the |t 2g and |e g orbitals.
In the case of CoO, whereas in weak crystal field theory, where Hund's rules are applied to the degenerate 3d orbitals to determine the 4 F (L = 3, s = 3 2 ) ground state free-ion electronic configuration, in the case of strong crystal field theory, Hund's rules are applied to the split 3d orbitals to determine the ground state orbital electronic configuration 2 E (s = 1 2 ). Armed with high quality single crystals and a high flux of epithermal neutrons present on the direct geometry chopper spectrometer MAPS at the ISIS spallation source [141,142], three distinct dd transitions at 0.870(9) eV, 1.84(3) eV, and 2.30 (15) eV were identified at 300 K [111]. Their identity being confirmed by both eliminating a non-electronic background component, and by noting the Q and temperature dependence of the excitations of interest.
Having experimentally determined the dd transitions, the parameterisation ofĤ CEF was accomplished by a comparison of these transitions to the energy levels of Co 2+ calculated from electrostatic matrices and character tables based on the formalism previously established by Tanabe and Sugano [143,144]. These electrostatic matrices are defined by both the crystal field splitting 10 Dq and electron-electron interaction parameters given by B, C or J(dd), C(dd) in Racah or Hubbard formalism, respectively [15,145,146].
As illustrated in Fig. 3, the experimental data was successfully reproduced [111] by the Tanabe-Sugano diagram with values for J(dd), C(dd), and 10Dq of 1.3(2) eV, 0.49(10) eV, and 0.94(10) eV, respectively, with both the values of J(dd) and its ratio to C(dd) showing excellent agreement with previously calculated values [105,146,147]. The ratio 10Dq J(dd) = 0.72 (18) confirmed that Co 2+ in CoO lies significantly below the spin-crossover value of 10Dq J(dd) ∼2.5, assuming a high spin (s = 3 2 ) orbital triplet 4 T 1 ground state, despite a value of ∼1 eV for 10Dq, placing CoO in the so-called weak field-intermediate crystal field regime [15]. In this particular regime, the complete set of commuting observables is identical to the set found in weak crystal field theory. In the case of magnetic properties, the values of L and s are fixed to those in the 4 F Co 2+ free-ion ground state (L = 3, s = 3 2 ). Diagonalisation ofĤ CEF for the 4 F free-ion ground state in an octahedral coordination (Eq. 7) yields an orbital triplet ground state 4 T 1 , an excited orbital triplet 4 T 2 , and an orbital singlet 4 A 2 , where ∆( 4 T 1 → 4 T 2 ) = 480B 4 and ∆( 4 T 2 → 4 A 2 ) = 600B 4 , where the Stevens factor B 4 is related to the crystal field splitting by 10Dq = 400B 4 [15,139,140,148,149].
The identification of the 1.84(3) eV excitation as the 4 T 1 → 4 A 2 transition was particularly noteworthy because there were two possibilities for the second excited state, corresponding to either 4 A 2 or 2 E [45,83]. Both transitions are dipolar forbidden with both exhibiting transition matrix elements | |0|M ±,z |f | = 0, whereM ±,z =L ±,z + 2Ŝ ±,z [1,150], and thus could not distinguished based on intensities alone in the dipolar approximation [151]. Instead, Cowley et al. [111] noted that the 4 T 1 → 4 A 2 transition corresponded to a large quadrupolar matrix element which had been previously predicted, but at the time had not been resolved experimentally [152,153].

Spin-orbit coupling
The low value of Z for Co 2+ places CoO in the Russell-Saunders L-S coupling scheme where spin-orbit coupling is treated as a pertubration of the formĤ SO = λL ·Ŝ to the |L = 3, m L , s = 3 2 , m s basis defined byĤ CEF in the low energy limit [16,17,31,129,130,149,154]. As was the case forĤ CEF , where the 4 F free-ion state was exclusively considered, a similar approach may be taken forĤ SO [137,138]. Since 10Dq ∼1 eV [111], the magnetic properties must solely be determined by the 4 T 1 ground state. Such a restriction of the analysis to the 4 T 1 ground state manifold requires a projection from the original |L = 3, m L basis onto a smaller basis |l = 1, m l , yielding a new projected spin-orbit Hamiltonian [130,155,156] given bŷ consisting of new orbital angular momentum operators that act on the projected |l = 1, m l , s = 3 2 , m s basis, accompanied by a projection factor α. Although it was the convention at the time to determine the value of α using representation theory [15], Sarte et al. [148,157], inspired by work done on 4d and 5d transition metal oxides by Stamokostas and Fiete [158], later developed an alternate method based on the matrix representation of angular momentum operators for obtaining these projection factors α.
Having projectedL onto a fictitious operatorl with l = 1, the basis forĤ SO (Eq. 8) is now the 12 |l = 1, m l ; s = 3 2 , m s states. Based on both the Landé interval rule and the addition theorem [151],Ĥ SO yields three unique effective total angular momentumĵ =l +Ŝ manifolds corresponding to j eff = 1 2 , 3 2 , and 5 2 , with their energy eigenvalues given by [35,111,149,159,160] Since the energy differences between j eff manifolds are fixed, λ (in theory) can be determined experimentally by measuring the energy transfer associated with spin-orbit transitions. To accomplish such task, Cowley et al. [111] performed neutron scattering experiments at low temperatures. The choice of neutron scattering limited the measurement of spin-orbit transitions to those transitions associated with ∆j eff = 0, ±1 [1]; while the use of low temperatures limited the measurement of spin-orbit transitions exclusively out of the ground state.

Tetragonal distortion & mean molecular field
As described so far, the experimental determination of λ appears deceptively simple. Such simplicity is based on two assumptions. The first is that the degeneracy for each j eff manifold is maintained, thus restricting ∆E = 0 transitions to those with ∆j eff = ±1. The second assumption is that the j eff manifolds are not shifted significantly away from the eigenvalues of H SO , such that ∆E remain proportional to λ. In the case of pure CoO, both assumptions are invalid [16,17,22,31].
As summarised by Fig. 4(b), the invalidity stems from the contributions from the two remaining terms inĤ 1 that have been conveniently absent in our discussion so far. The first is the distortion HamiltonianĤ dis , which corresponds to the structural deformation of the unit cell that accompanies long-range antiferromagnetic order [11][12][13][51][52][53]. While the exact symmetry of the low temperature phase has proven to be particularly contentious [37], we have chosen to consider the case of a tetragonal distortion which we will show captures the essential low energy physics. We note that our inelastic data is not able to uniquely identify the nuclear and magnetic structures. Utilising symmetry arguments, the influences of a tetragonal distortion is given by [140,149,[161][162][163] where Γ corresponds to a distortion parameter, and whose presence inĤ 1 results in the removal of the degeneracy of the spin-orbit manifolds [164][165][166]. The second additional term is the molecular field HamiltonianĤ M F . Behaving as a Zeemanlike term,Ĥ M F breaks time reversal symmetry, thus removing the degeneracy of the individual j eff manifolds [149]. As a first approximation, by considering exclusively a single dominant next-nearest neighbour 180 • Co 2+ -O 2− -Co 2+ superexchange pathway with a magnetic exchange constant J 2 ,Ĥ M F (Eq. 2) can be simplified to [22,129,130] where z 2 denotes the number of next-nearest neighbours. As illustrated in Fig. 4(b), a strong value of this exchange interaction results in significant mixing of spin-orbit split levels [16,17]. In the case of CoO, such entanglement is further exacerbated by the presence of various long range spin-spin superexchange interactions with exchange constants comparable to λ [27,[31][32][33]. With the additional complications of both complex magnetic ordering and structural distortions [12,13,37,167], it quickly becomes apparent how such a deceptively "simple" 3d monoxide yields such a complex magnetic excitation spectrum. The limitations imposed on approaches based on conventional linear spin wave theory [23][24][25][26], combined with the necessity of a multiparameter spin-orbital Hamiltonian incorporating both exchange and spin-orbit coupling to describe such a spectrum, made both the modelling and understanding of the magnetic excitations in CoO particularly elusive for decades [28,31,119,[168][169][170]. In order to circumvent the problematicĤ M F contribution, it is important to note that the significant degree of entanglement present in CoO was not inherent to all Co 2+ magnets, with some possessing a small H M F due to |J| being significantly smaller than |λ| [129, 130, 148-150, 161, 171]. In contrast to CoO, these Co 2+ -magnets exhibited a much weaker degree of mixing, such that the spin-orbit split levels remained well-separated in energy. Although limited by the large value of J 2 in CoO, Cowley et al. [111] effectively eliminated H M F by the chemical dilution of Co 2+ by non-magnetic Mg 2+ . The resulting dilute monoxide Mg 0.97 Co 0.03 O, in which the majority of Co 2+ are isolated from one another, does not assume long range magnetic order down to 5 K. Furthermore, as its unit cell remains cubic upon cooling (Ĥ dis = 0),Ĥ 1 for the case of Mg 0.97 Co 0.03 O can be reduced toĤ CEF +Ĥ SO (Fig. 4(c)).
Neutron scattering experiments on a polycrystalline sample of Mg 0.97 Co 0.03 O at 5 K on the direct geometry chopper spectrometer MARI [172,173] at ISIS identified one distinct excitation at 37.1(5) meV [111]. This particular excitation exhibited a Q dependence consistent with the Co 2+ magnetic form factor, thus confirming its origin as electronic, while its FWHM of 6(1) meV being significantly broader than the instrumental resolution ∼ 2% ∆E E i ∼ 2 meV [174] is consistent with a distribution of different local environments that would be expected from such a dilute monoxide.
Having established that the excitation was indeed the j eff = 1 2 → j eff = 3 2 spin-orbit excitation and given the gap ∆E = −9λ 4 [88,129,130,148,161], an energy transfer of 37.1(5) meV yields a value for λ of −16(3) meV. Despite not requiring any detailed knowledge of the Q dependence of the spin-excitations, the inclusion of multiple exchange constants, and a structural distortion contribution, the value of λ determined through chemical dilution was in fact consistent with values obtained from nominally pure CoO through indirect and model dependent approaches [15,31,78,167,169,[175][176][177][178][179].

Extraction of Exchange Constants J
Having parameterised the two largest contributions toĤ 1 , Roger's attention now shifted towards addressing the cooperative magnetism in CoO. Faced with both the aforementioned limitations imposed on conventional linear spin wave theory approaches and a body of often contradictory literature spanning decades [16, 17, 22, 27, 31, 34, 76-78, 82, 169], Roger knew that the problem of extracting J would need to be accomplished through an alternative approach.
Motivated by this alternative, and seemingly more "direct" method for extracting J [180,196], additional inelastic neutron scattering measurements were performed on polycrystalline Mg 0.97 Co 0.03 O [88]. As illustrated in Fig. 4(a), measurements at 5 K displayed a hierarchy of seven dispersionless magnetic excitations up to an energy transfer ∆E ∼15 meV. These excitations exhibited clear modulation in Q that is characteristic of magnetic clusters [196][197][198][199], thus distinguishing them from single-ion dispersionless crystal-field excitations that were previously observed at higher energy transfers [111]. It was at this time that Roger was faced with two important questions. The first was what were these magnetic clusters? It was clear that these excitations were not single-ion in origin, but it was not clear if these excitations were from interactions between pairs (N = 2), triples (N = 3), etc. The second was how would the J constants be extracted? In contrast to d 5 Mn 2+ , the magnetism for d 7 Co 2+ is arguably made much more complicated by the orbital degree of freedom in the t 2g channel, requiring multiple projections of angular momentum operators [15-17, 130, 148, 155, 156, 159].
In this section, we describe how Roger's collaborators Sarte et al. [88] exploited the paramagnetic nature Mg 0.97 Co 0.03 O to disentangle the spin exchange and spin-orbit interactions. It would take almost three additional years for this disentanglement, combined with the interpretation of the energy hierarchy of excitations as Co 2+ pairwise interactions, to result in the extraction of 7 exchange constants out to the fourth coordination shell in the rocksalt lattice.

Sample characterisation and dominance of Co 2+ pairs
Although the value of 3% Co 2+ concentration at first appears nonsensical, especially in the context of an often flux-limited technique such as inelastic neutron scattering [1], this particularly low concentration actually significantly simplifies the magnetism under consideration by restricting the cooperative magnetism in Mg 0.97 Co 0.03 O to be dominated by pairwise interactions between isolated Co 2+ pairs [180]. Such dominance can be rationalised from probabilistic arguments, where it can be shown that the ratio of numbers of N to N + 1 spin clusters in Mg 1−x Co x O is given by 1−x x ; thus, for small values of x, the number and hence scattering from Co 2+ pair excitations far outweigh those from larger clusters, as has been observed with dilute Mn 2+ -based magnets for x < 0.05 [186,188].
Although such arguments provide an excellent rationale for the dominance of pairwise interactions in Mg 0.97 Co 0.03 O, these probabilistic arguments are ultimately based on the assumptions that the Co 2+ concentration is actually x = 0.03 and its distribution is homogenous. In order to confirm the value of x, its value was experimentally determined using both x-ray diffraction and DC magnetometry.
The question of chemical homogeneity and the possibility of Co 2+ clustering was addressed by first noting the absence of a ZFC/FC split in the DC susceptibility. Secondly, the powderaveraged dynamic structure factor S(Q, E) for Mg 0.97 Co 0.03 O was shown to be exhibit completely diffferent features compared to those present in either CoO and MgO. A third and final argument is that the value of λ that was extracted from Mg 0.97 Co 0.03 O agreed very well with the reported values in the literature [15,31,111,169,175]. An agreement only made possible by the disentanglement of the individual spin-orbit manifolds, a process that would hindered by any type of clustering [180].
Direct experimental confirmation of the absence of any Co 2+ clustering was accomplished by two vastly different methods. The first involved performing the same scattering experiments at identical temperatures with slightly different incident energies on a second polycrystalline sample of Mg 0.97 Co 0.03 O that was synthesised by a solution-based (sol-gel) technique [200] in lieu of the traditional ceramic approach [201,202]. The second method consisting of EDX-SEM measurements [203] confirmed the uniform distribution of Co 2+ , Mg 2+ , and O 2− , with no evidence for significant clustering. The resulting elemental analysis confirmed that the Co 2+ concentration was within experimental error to the values previously deduced by both DC magnetometry and x-ray diffraction measurements.

Effective pair spin-orbital Hamiltonian
Having established that the cooperative magnetism and thus the low energy response in Mg 0.97 Co 0.03 O is dominated by Co 2+ pairs and their pairwise interactions, the discussion now shifts to the neutron scattering response of an isolated pair of magnetic ions and how it can be used to extract both the interaction distance R and exchange constant J.
As summarised by Fig. 4(c), the low energy fluctuations in Mg 0.97 Co 0.03 O correspond exclusively to excitations within the ground state doublet, and thus requiring the projection of the spin operatorŜ onto the j eff = 1 2 manfiold [16,17,31,129,130,148,150,189]. By defining the projected angular momentum operatorŜ = βĵ with j = 1 2 , the interaction energy between a Co 2+ spin pair given byĤ ex = 2JŜ 1 ·Ŝ 2 can be approximated bŷ where the effective projection factorα = 2β 2 is a proportionality constant between the magnetic exchange constant J and the gap between the triplet (Γ eff = 1) and singlet (Γ eff = 0) states being given by ∆E =α|J| [180,196]. This value ofα can be calculated numerically by diagonalising an effective pair Hamiltonian describing the superexchange interaction of the spin-orbit split manifolds of a Co 2+ pair in the non-distorted (Ĥ dis = 0) and magnetically dilute paramagnetic (Ĥ M F = 0) Mg 0.97 Co 0.03 O lattice [111]. As predicted by the projection theorem of angular momentum [15,151,204], the energy splitting within the j eff = 1 2 manifold calculated fromĤ pair (Fig. 4(d)) is linear when |λ| |J| with α = 50 9 [22,130]. Therefore, pair excitations in dilute Mg 0.97 Co 0.03 O provides a direct experimental route for determining |J|. The limitation to a magnitude is based on principle that ∆E is simply an energy difference between the two eigenvalues of the ground state doublet, corresponding to either a transition from Γ eff = 1 to Γ eff = 0 state or vice versa, and thus its value is independent of the sign of J [205]. As will be discussed later, the determination of the sign of J is more subtle, ultimately relying on statistical mechanical arguments.
While the excitation energy provides |J|, the Q dependence can be used to extract the corresponding intra-pair distance R m , where m denotes the coordination shell. By applying the powder-averaged first moment sum rule and the single mode approximation, the Q dependence of the energy-integrated dynamic structure factor for an isolated Co 2+ pair is proportional to |F (Q)| 2 1 − sin(QRm) QRm [149,[206][207][208][209][210][211][212]. Since the modulation is solely dependent on R m , an excitation can be assigned to a particular pair and corresponding coordination shell [196,213].

Experimental determination of J
Having established the theoretical framework to describe the excitations from the Co 2+ pairwise interactions in Mg 0.97 Co 0.03 O, the extraction of the values for J from the low energy hierarchy that was measured [88] on MARI and IRIS will now be addressed. As illustrated in Fig. 4(d), exchange constants were assigned to each of seven excitations based on the value of excitations' measured energy transfers via the diagonalisation ofĤ pair (Eq. 13) using the extracted value of λ = −16(3) meV [111]. By fitting the energy-integrated intensity of the excitations to the first moment sum rule [206] (Fig. 5(a)), each of the 7 exchange constants was assigned to one of the relative coordination shells ranging from m = 1 to m = 4.
Inspired by neutron scattering studies on the dimer compound Ba 3 Mn 2 O 8 [214][215][216], an experimental route based on the temperature dependence of the integrated intensity was used to determine the sign of J. Since the integrated intensity scales as the thermal population difference between the ground and excited states [205,210,213,214], the temperature dependence of the integrated intensity for antiferromagnetically coupled (J > 0) pairs and its ferromagnetic counterpart are expected to exhibit distinctly different dependence on temperature. As illustrated in Fig. 5(b), it was shown that all integrated intensities fall onto either one of two universal curves describing antiferromagnetism or ferromagnetism. The extracted values of J based on the energy, momentum, and temperature dependence are summarised in Tab 1.

The "dual" hierarchy
Based on the literature available at the time [217,218], the observation that all coordination shells, with the exception of m = 2, display two closely spaced excitations with differing signs for J was initially confounding. Although unexpected, the presence of "dual" ferroand antiferromagnetic interactions is in fact consistent with the GKA rules [219][220][221][222][223], since each of these exchange pathways consists of at least one 90 • Co 2+ -Co 2+ exchange interaction involving the overlap of half and filled t 2g orbitals. As illustrated in Fig. 4(d,inset), the GKA rules indeed predict that the combination of the orbital degree of freedom for each Co 2+ and a lack of orbital ordering (or anisotropy) in Mg 0.97 Co 0.03 O would manifest itself as either a direct antiferromagnetic t 1 2g -t 1 2g or a weaker ferromagnetic t 1 2g -t 2 2g direct exchange interaction. As is summarised in Tab. 1, the experimental values of J are not only consistent with the GKA rules, as the antiferromagnetic interaction is stronger than its ferromagnetic counterpart for all the m = 2 excitations, whilst the 180 • Co 2+ -O 2− -Co 2+ m = 2 coupling leads only to a strong antiferromagnetic interaction, but futhermore,the values of J show excellent agreement with general trends reported by both experiment [27,31,34,[76][77][78][79][80] and recent GGA+U DFT calculations [47].
Finally, the mean field definition of the Curie-Weiss temperature given by [148,224,225] enables a direct comparison of the extracted values of J to thermodynamic data for both Mg 0.97 Co 0.03 O and CoO. In the case of the former, the concurrent presence of AFM and FM implies that the value of J 2 will ultimately determine θ CW . Therefore, by inserting the values of s = 3 2 , J 2 = 35.9(6) K (3.09(5) meV), setting the value of z 2 to 1 since the chemical dilution reduces the cooperative magnetism to individual pairwise interactions, and dividing by a correction factor of ∼ 1.9 that was deduced by Kanamori [16] to account for spin-orbit coupling, Eq. 14 yields a value of −44.9(7) K, in agreement with the experimentally determined value of −41(6) K [88].
In the case of CoO, such a comparison is slightly more involved, and involves the key assumption that CoO assumes a type-II collinear antiferromagnetic structure, corresponding to (111) ferromagnetic sheets that are stacked antiferromagnetically along [111] [11-13, 47, 55, 58]. As illustrated in Fig. 6, the determination of z i was accomplished by assigning F and AF labels to Co 2+ cations that are located on an even integer (or same) number and an odd number of (111) planes away from a reference Co 2+ , respectively. By inserting the values of s = 3 2 , J i as determined from the measured energy transfers, z i from the procedure outlined above, and Kanamori's second perturbation theory correction factor into Eq. 14, a Curie-Weiss temperature of −295(5) K is obtained, in good agreement with both the experimental value of θ CW = −330 K [82,[226][227][228], and in particular, T N = 291 K [14,16,17,22,37,[229][230][231].
As will be discussed in the next section, the determination of J represents the parameterisation of the last key contribution to Roger's minimalist Hamiltonian (Eq. 1), laying the foundation for the model that would ultimately be used to address the low energy magnetic fluctuations of CoO.

Spin-Orbit Excitons in CoO
In this section, we describe how Roger's collaborators Sarte et al. [119] extended previous theoretical work on PrTl 3 [131] to construct a mean-field and multilevel spin-orbit exciton model based on Green's functions to address the complex magnetic excitation spectrum of CoO. By employing the spin-orbit and spin exchange coupling parameters that were previously determined experimentally [88,111], the model, based on a tetragonally distorted type-II antiferromagnetic unit cell, successfully captured both the sharp low-energy excitations at the magnetic zone centre, and the energy broadened peaks at the zone boundary. Despite the model's success at low energy transfers, the failure of the model to describe the higher energy excitations leaves an important avenue open for future investigations.
for the response function G αβ of generic spin operatorsÂ andB that is defined as and proportional to the magnetic neutron cross section by the fluctuation-dissipation theorem [1,213,[232][233][234][235]. It is the presence of [Â,Ĥ] in Eq. 15 that necessitates a thorough understanding ofĤ, its commutator withÂ, and why its parameterisation has dominated the conversation so far. As first outlined by Buyers et al. [131] and described later in more detail by Sarte et al. [119], the derivation of an expression for the neutron response function begins by first rotating all components of the spin operators comprising the inter-ion HamiltonianĤ 2 onto a basis consisting of the eigenstates ofĤ 1 . It can be shown that such a rotation reduces Eq. 15 to four commutators, and when combined with the random phase decoupling method [236][237][238][239] , allows the equation-of-motion to be written as describing the coupling of the single-site response function g αβ (ω) by the Fourier transform of the exchange interaction J(Q) [240]. With the only nonzero single-site response functions in such a highly symmetric environment being: g +− , g −+ , and g zz [130,132,241], the restriction of αβ combinations to +−, −+, or zz allows Eq. 17 to be rewritten as where the prefactor Φ = 1 when α = +,− or 2 when α = z. By approximating CoO as a type-II collinear antiferromagnet [11-13, 47, 55, 58]., the site indices i and j are restricted to values 1 or 2, yielding four coupled equations for each αβ combination. These equations when combined together yield where G −+ (Q, E) has the same form as G +− (Q, E) with indices + ←→ −. Here, J s and J d denotes J(Q) on the same (i = j) and different (i = j) sublattices (Fig. 6), respectively, while ω has been relabelled as E =hω, sinceh is conventionally set to 1. Finally, by the fluctuation-dissipation theorem, the imaginary part of the total response function G(Q, E) = G +− (Q, E) + G −+ (Q, E) + G zz (Q, E) is proportional to the dynamical structure factor [242][243][244][245], thus providing a direct comparison between the calculated model and experiment.

Model
As presented in Eq. 19, G(Q, E) is a function of both g αβ (E) and J(Q). Since g αβ is itself a function ofĤ 1 (Eq. 4), the single-site response in the limit of T → 0 K is defined by three parameters: λ, Γ, and H M F . An initial estimate for λ was taken to be its reported value of −16 meV [111]. An initial estimate for H M F was determined from the experimentally determined θ CW [82,[226][227][228] via its mean field definition (Eq. 14), yielding an initial estimate of 64.8 meV. The initial estimate of the distortion parameter Γ=−8.76 meV was determined by scaling the reported value of Γ for KCoF 3 [130] by the ratio of their respective tetragonal distortions.
In the case of J(Q) = i =j J ij e iQ·d ij , it is both a function of J and the magnetic structure under consideration [240,246]. As a first approximation, CoO was treated as a tetragonally distorted type-II collinear antiferromagnet with isotropic exchange constants equal to those reported for Mg 0.97 Co 0.03 O [88]. In order to account the possibility for both the "dual hierarchy" and the tetragonal distortion, 16 different orbital configurations of the form xAxxγ were considered. Whereas x could either be antiferromagnetic (A) or ferromagnetic (F) to account for the possibility for either anti-/ferromagnetic coupling in coordination shells m =1, 3, and 4, the index γ could assume labels 1 or 2, distinguishing the presence or absence of distorted bonding configurations. Finally, in order to account the large experimental beam present on MERLIN [247], the model considered the mean contribution from all 16 equally weighted xxAxγ orbital configurations, each with the same value of λ, Γ, and J n that are subject to different values of H M F , corresponding physically to a unique type of "domain" in the bulk CoO single crystal.

Comparison between model & experiment: success and failures
As illustrated in Fig. 7 and summarised in Tabs. 2 and 3, by allowing the value of H M F to refine independently for each of the equally weighted 16 xAxxγ domains, with each possessing identical refined values of λ, J, and Γ, the spin-orbit exciton model [119,131] successfully reproduced both the fine structure at the magnetic zone centre and the broad excitations at the (1.5, 1.5, −1) zone boundary, whilst capturing the steeply dispersive columns of scattering observed at higher energy transfers.
The need for all 16 xAxxγ domains stems from the low energy fine structure at the zone centre [31][32][33]. Although an individual xAxxγ domain does reproduce the correct bandwidth, the neutron response is dominated by a single highly dispersive G −+ mode, while both G −+ and G zz modes are significantly weaker in intensity and weakly dispersive at higher energy transfers. By considering the fine structure to consist of multiple overlapping G −+ modes exclusively, other additional features of the data are captured by the model including: the flat band at ∼40 meV, the weak dispersionless mass of scattering at ∼80-90 meV, and the steeply dispersive columns of high energy scattering at the (1.5, 1.5, −1) zone boundary.
In addition to the necessity for 16 xAxxγ domains, the discrepancy between the calculated dispersion and that measured experimentally indicated the need for the optimisation of the model's parameters. Corresponding to the main parameter forĤ pair , λ plays a key role in the determination of |J|, in particular |J 2 |, from the experimentally determined ∆E. In the case of a fixed value of H M F , an increase in the value of |λ| corresponds to a decrease in |J| for a given ∆E. The result is a weaker dispersion that shifts the dominant G −+ component to higher energy transfers at the zone centre, while minimally affecting the two other weakly dispersive modes. Such behaviour contrasts the effect of H M F which is not on the dispersion itself, but rather a change in its value results in a systematic shift in the value of energy transfers for a fixed Q, with the shift being significantly larger compared for the same relative change for |λ|. Finally, in contrast with λ (J) and H M F , the influence of the Γ is most pronounced G +− , providing a mechanism to shift the G +− mode without inducing a comparable shift for the prominent G −+ mode when all other parameters are fixed.
With such a large number of domains under consideration, constraints on the parameter space for J via λ, Γ, and H M F were required to ensure convergence for a least squares optimisation. The parameter spaces for both λ and Γ were defined by the values previously reported in the literature [15,31,78,111,130,162,167,169,[175][176][177][178][179]. Initial modelling attempts fixed the values of J to those obtained from the diagonalisation ofĤ pair for a fixed λ [88]. The failure of these initial attempts to reproduce both the fine structure at the zone centre and the broad excitations at the (1.5, 1.5, −1) zone boundary simultaneously was ultimately remedied by allowing the values of J to vary ±20% from its value obtained fromĤ pair . In contrast to all the other parameters under consideration that were set to be equal for all 16 domains, no such restriction was applied to H M F which was allowed to vary from 0 to an arbitrarily large upper limit.
Despite the success of the spin-orbit exciton model [119,131] to account for some of the experimental data, its success appeared to be limited to low energy transfers along (1.5, 1.5, L). Along (2, 2, L), both the spin-orbit exciton model and linear spin wave theory predicted nearly zero intensity for magnetic fluctuations for an antiferromagnetic structure, in stark contrast with the prominent magnetic scattering that was observed experimentally (Fig. 7) [22,126]. On the other hand, although the model did successfully reproduce the steep columns of scattering at higher energy transfers at the zone boundaries along (1.5, 1.5, L), a closer inspection of the Q dependence of these high energy modes (Fig. 8(a) revealed that their intensity decayed more rapidly than what is predicted by the magnetic form factor f (Q), indicating the presence of delocalised magnetism [248,249]. While the similarity between the dispersion of phonons at high Q (Fig 8(b)) and the magnetic scattering absent in the spin-orbit exciton model suggests that a magneto-vibrational contribution [250][251][252] may account for the model's failures along (2, 2, L), there is no clear physical mechanism to account for the failure at high energy transfers. Promising possibilities include spatial extended magnetism due to strong covalent bonding and significant hybridisation of the 3d orbitals that was found to be the case for Sr 2 CuO 3 [253], or possibly multi-magnon decay processes, supported by distinct similarities of the high energy response in CoO to those for itinerant magnets [109,110,[254][255][256][257][258][259][260][261][262].

Concluding Remarks
The story of CoO is of one that has come full circle. Here we have described how Roger, through his extraction of the single-ion parameters of Co 2+ via chemical dilution [88,111], and Bill, through the establishment of the spin-orbit exciton model [131], laid for the foundation for the future success of Sarte et al. [119] in parameterising the low energy magnetic fluctuations in CoO that were first measured with neutrons by these two at Chalk River more than 50 years ago [22]. While the failure of the model at the zone boundaries can be attributed to a magnetovibrational contribution to the neutron cross section [250][251][252], the physical origin underlying the rapid decay of the column of magnetic fluctuations at high energy transfers remains an open question. Despite the strong insulating nature of CoO [45], it can be speculated that the model's failures at high energy transfers corresponds to the breakdown of spin-orbit excitons which may be accompanied by a crossover from localised to spatially-extended magnetism, reminisicent of an itinerant-like response [249] or strong covalency [253].
The success of the spin-orbit exciton model [119,131] with such a simple Hamiltonian in addressing the low energy magnetic fluctuations for a system with an orbital degree of freedom that is far displaced from the λ J limit such as CoO, suggests that such a model would be of great interest for a community who has recently paid particular interest on Mott insulators with strong spin-orbit coupling in the search for unconventional, and often novel magnetic ground states [97,98,[263][264][265][266][267]. The parameters extracted from this analysis will be of importance in future work in understanding the response of exchange bias in thin films and also in resolving the low temperature nuclear and magnetic structures in CoO. Representing a powerful tool that allows for the direct parameterisation of a simple Hamiltonian to model the excitation spectrum under consideration, it would be of interest to see if the spin-orbit exciton model [131] would be able to achieve the same success in addressing low energy magnetic fluctuations in high Z magnets (e.g. Ru 4+ , Ir 4+ ) as calculations currently reported in the literature that are based on linear spin wave theory, often employing very complex Hamiltonians [112][113][114][115][116][117][118][120][121][122][123][268][269][270]. Figure 1. Comparison of Q-integrated cuts of inelastic neutron scattering data on single crystal CoO reported in the literature [22,[31][32][33] and recent scattering measurements performed on MERLIN [247] at 5 K by Sarte et al. [119] with an E i = 75 meV and 45 meV at the magnetic zone (a) boundary and (b) centre, respectively. Solid lines in (a) indicate the location of excitations previously determined by IR spectroscopy [176]. Horizontal bars indicate instrumental resolution [174]. (Q, E) slices of CoO measured on MERLIN at 5 K with an E i of (c) 110 meV, (d) 75 meV, and (e) 45 meV reported by Sarte et al. [119] illustrates an excitation spectrum that is highly structured in both energy and momentum. All (Q, E) slices have been folded along [001]. Adapted from Reference [119] with permission.  [143,144] for d 7 Co 2+ in octahedral coordination calculated by Cowley et al. [111]. Shaded rectangles correspond to excitations for cubic CoO experimentally measured on MAPS [141] at room temperature, with the height and width corresponding to experimental error in energy and the statistical error of the refined value for 10Dq/J(dd), respectively. The red dashed line at 10Dq/J(dd) ∼ 2.5 denotes the spin crossover from (left) a high-spin S = 3 2 , 4 T 1 orbital configuration to (right) the low-spin S = 1 2 , 2 E. Adapted from References [111,119] with permission. with an E f of 1.84 meV revealing seven low-energy dispersionless magnetic excitations [88]. (b) Calculated energy eigenvalues as a function of the tetragonal distortion (Ĥ dis ) and the magnetic molecular field (Ĥ M F ) perturbations to the j eff manifolds of the 4 T 1 ground state for d 7 Co 2+ in octahedral coordination. Both the energy eigenvalues and individual parameters are presented to scale. (c) Relevant energy scales for the effective pair HamiltonianĤ pair (Eq. 13) [130]. (d) Energy splitting of the j eff = 1 2 manifold calculated from the diagonalisation ofĤ pair (black line). The non-linearity constrasts the behaviour predicted by the projection theorem [15,151,204] (gray line). (inset) The mechanism for antiferromagnetism (top) and weaker ferromagnetism (bottom) is a result of a combination of the 90 o Co 2+ −O 2− −Co 2+ exchange pathway and the orbital degree of freedom in the t 2g channel on each Co 2+ , in agreement with the predictions of the Goodenough-Kanamori-Anderson rules [219][220][221][222][223]. Yellow arrows denote local t 2g spin configurations and teal arrows denote total spin configurations on each Co 2+ . Adapted from References [88,119] with permission. Figure 5. (a) Form-factor-corrected neutron scattering intensity Q dependence for all seven magnetic excitations previously identified in Fig. 4. The Q dependence for each magnetic excitation has been rescaled by the intradimer distance R obtained from the first moment sum rule [149,[206][207][208][209][210][211][212]. The black curve corresponds to 1 − sin(QR) QR . (b) Normalised emperature dependence of the integrated intensity for the low energy magnetic excitations. All seven excitations fall onto one of two universal curves calculated for antiferromagnetic and ferromagnetic coupled pairs. All quantities were normalised by the analytical expression for the temperature dependence of ferromagnetically coupled pairs, I F (T ) [205,210,213,214]. Adapted from Reference [88] with permission.   . Arrows indicate fluctuations that may exhibit itinerant-like behaviour, as (b) a Q-integrated (ξ,ξ,L) cut of (a) falls off more rapdily than the Co 2+ magnetic form factor, in stark contrast with the behaviour of (c) an equivalent Q-integrated (ξ,ξ,L) cut for an E i = 45 meV. Comparison between the phonon scattering at large Q centred about (d,f) the nuclear zone boundaries and (e) modes located near the magnetic zone boundary along the (2, 2, 0) direction that were unaccounted for by the spin-orbit exciton model [119]. Dashed lines indicate the overlap of the energy transfers range between both gapped optical phonon modes centered near the nuclear zone boundaries and the low-Q scattering near the (2, 2, 0) magnetic zone boundary. All three (Q, E) slices have been folded along [001] and have been renormalised to share a common relative intensity scale. Adapted from Reference [119] with permission.