# Single-ion properties of the transverse-field Ising model material CoNb2O6

J. A. Ringler  Department of Physics, Colorado State University, Fort Collins,
CO 80523, USA  A. I. Kolesnikov  Neutron Scattering Division, Oak Ridge
National Laboratory, Oak Ridge, Tennessee 37831, USA  K. A. Ross  Department
of Physics, Colorado State University, Fort Collins, CO 80523, USA  Quantum
Materials Program, CIFAR, MaRS Centre, West Tower 661 University Ave., Suite
505, Toronto, ON, M5G 1M1, Canada

###### Abstract

CoNb2O6 is one of the few materials that is known to approximate the one-
dimensional transverse-field Ising model (1D-TFIM) near its quantum critical
point. It has been inferred that Co2+ acts as a pseudo-spin 1/2 with
anisotropic exchange interactions that are largely Ising-like, enabling the
realization of the TFIM. However, the behavior of CoNb2O6 is known to diverge
from the ideal TFIM under transverse magnetic fields that are far from the
quantum critical point, requiring the consideration of additional anisotropic,
bond-dependent (Kitaev-like) terms in the microscopic pseudo-spin 1/2
Hamiltonian. These terms are expected to be controlled in part by single-ion
physics, namely the wavefunction for the pseudo-spin 1/2 angular momentum
doublet. Here, we present the results of both inelastic neutron scattering
measurements and electron paramagnetic resonance spectroscopy on CoNb2O6,
which elucidate the single-ion physics of Co2+ in CoNb2O6 for the first time.
We find that the system is well-described by an intermediate spin-orbit
coupled Hamiltonian, and the ground state is a well-isolated Kramers doublet
with an anisotropic $g$-tensor. We provide the approximate wavefunctions for
this doublet, which we expect will be useful in theoretical investigations of
the anisotropic exchange interactions.

###### pacs:

xx.xx.mm

††preprint: APS/123-QED

##  I Introduction

The transverse-field Ising model (TFIM) is arguably the simplest and most
tractable model that exhibits a quantum phase transition, yet it leads to a
wealth of intriguing results [1]. Of particular current interest, this model
has been used to predict non-equilibrium phenomena such as Kibble-Zurek
scaling, quantum annealing in glassy spin systems, and potential violations of
the eigenstate thermalization hypothesis and the emergence of the subsequent
many-body localized phase [2, 3, 4]. Materials that realize the TFIM are
highly sought-after for both testing of such theoretical ideas and materials
engineering. However, the requirement of having a low enough (e.g. laboratory-
accessible) transverse magnetic field strength to induce the quantum critical
point (QCP) strongly limits the number of TFIM materials. There are only three
known solid-state magnetic systems which approximate the TFIM, including
dipolar-coupled 3D LiHoYF4 [5], as well as quasi-1D systems (Ba/Sr)Co2V2O8 [6]
and CoNb2O6 [7], the latter being the system of interest in this study.

The quasi-1D TFIM magnet CoNb2O6 is known to exhibit signatures of a 1D QCP
111It should be noted that the 3D transition actually occurs at
$\mu_{0}H_{c}=$5.5\text{\,}\mathrm{T}$$, as observed by inelastic neutron
scattering measurements [7]. However, the 1D QCP at zero temperature can be
inferred from properties of the quantum critical fan [12]. under an applied
transverse field strength of $\mu_{0}H_{c}=$5.25\text{\,}\mathrm{T}$$, as
observed via various thermodynamic measurements such as specific heat and
magnetization, as well as nuclear magnetic resonance [9, 10, 11, 12]. CoNb2O6
belongs to a family of $3d$ transition metal niobates $M$Nb2O6 which form in
the columbite structure (space group Pbcn (60), $M$ = Cd, Co, Cu, Fe, Mn, Ni,
Zn) [13, 14], where zig-zag chains of edge-sharing $M$-O^2- octahedra run
along the crystallographic c-axis (Fig.1(a)). For Co^2+ in this structure, the
crystal electric field (CEF) of the local environment and spin-orbit coupling
(SOC) conspire to form effective spin-$\frac{1}{2}$ degrees of freedom with
magnetic moments that tend to lie in the ac-plane (in the absence of a
transverse field). These pseudo-spins form strongly-coupled ferromagnetic (FM)
chains along c [15], the details of which are discussed in Sec. III. The
chains are linked together antiferromagnetically (AFM) in the ab-plane (Fig.
1(b)), and the b-axis is perpendicular to the net ordered moments. The
pseudospins have been assumed to be coupled via Ising interactions, a model
that correctly reproduces many of the features near the critical point with a
field along the transverse axis.

Figure 1: Visualization of the crystal and magnetic structures of CoNb2O6. (a)
For CoNb2O6 in the columbite structure, the Co^2+ is octahedrally coordinated
with O^2- ions, which run in a zig-zag pattern along the c-axis, as described
in the main text. (b) The effective spin-$\frac{1}{2}$ moments orient in the
ac-plane and form an isosceles triangular network with AFM exchange in the ab-
plane.

However, away from the critical point there have been some notable deviations
from ideal TFIM behavior, which until recently had been explained only
phenomenologically [16, 4, 17, 7, 18, 19]. Recently there has been a renewed
interest in the microscopic details of the Hamiltonian in the regime far from
the QCP. A recent study by Fava et al. [19] demonstrated that the symmetry
being broken near the QCP is not that of a global Ising symmetry, but that of
a glide symmetry which is a consequence of its 3D space group. This symmetry
analysis showed that there are off-diagonal terms in the Hamiltonian beyond
the usual Ising exchange that are needed to account for all discrepancies
between the model and experimental data. Shortly thereafter, Morris et al.
[20] modeled their time-domain terahertz (THz) spectroscopy data with a
Hamiltonian that explained the overall ordered moment direction as the average
of two different observed Ising exchange axes in the $ac$ plane, which
alternate along the c-axis. These bond-dependent Ising interactions inspired
the authors to form a new microscopic Hamiltonian the authors dubbed the
“twisted” Kitaev chain, which is another model of fundamental importance in
non-equilibrium physics and contains off diagonal interactions equivalent to
those proposed for CoNb2O6in Ref. 19. Additional THz measurements supported
their findings and ultimately led to the possible first direct measurement of
Kramers-Wannier duality in nature [21]. In order to match their field-
dependent terahertz data, Morris et al. proposed an orientation of the Ising
exchange axes; the best fit had the same ordered moment direction in the
$ac$-plane, but additionally included a component that alternated $\pm$ 17
degrees along the $b$-axis. However, their fit did not take into account the
$g$-tensor, which should play an important role in this determination, since
it combines with the Ising exchange axes to determine the final ordered moment
direction.

In both of these recent works proposing terms beyond the simple Ising
ferromagnet, the models map back to the TFIM and produce quantum critical
behavior within the same Ising universality class as the 1D TFIM. However,
questions remain about these new Kitaev-like terms that become increasingly
relevant farther from the QCP.

Preceding these efforts in CoNb2O6, there has been a growing interest in
understanding bond-dependent anisotropic interactions (such as Kitaev
interactions) in $3d$ transition metal compounds, particularly Co^2+-based
magnetic systems [22, 23, 24]. Liu et al. recently showed that $3d$ cobaltates
should display bond-dependent interactions, and in some situations the Kitaev
term would dominate the microscopic Hamiltonian [22]. These calculations were
applied to “honeycomb lattices”, formed by layers of edge-sharing, trigonally-
distorted CoO6 octahedra. The authors parameterized the Co^2+ CEF ground state
wavefunctions by the trigonal distortion strength, which determines the
strength of the non-Kitaev terms and leads to a rich phase diagram. Although
CoNb2O6 is not a honeycomb lattice material, and the octahedra are not simply
trigonally distorted, the work of Liu et al. inspires the notion that the
details of the Kitaev-like interactions in CoNb2O6 can be understood by
studying the single-ion wavefunctions of Co^2+ in the CoNb2O6 lattice. This is
part of the motivation for this work.

In this article, we present the results of inelastic neutron scattering (INS)
measurements on pure CoNb2O6 which used a wide incident neutron energy range
$E=60-$2500\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$, and use an intermediate
spin-orbit coupled model to determine its single-ion energy levels (Sec. IV.1)
as well as the wavefunctions of the Co^2+ CEF ground state doublet (App. A).
We then present the EPR spectrum measured from a diluted variant (described in
Sec. II). These data conspire to definitively show that the Co^2+ single-ion
ground state is a well-isolated Kramers doublet with strongly anisotropic
moments described by a rhombic $g$-tensor. Although this scenario is expected
based on prior work, this study provides the details that were previously
missing, and which can be used to develop a theoretical understanding of the
anisotropic exchange in CoNb2O6; namely a direct determination of the
transverse field $g$-value, $g_{b}\approx 3.3$, and the details of ground
state CEF wavefunctions.

The contents of the paper are structured as follows: In Sec. II, we describe
the experimental details of all measurements. In Sec. III, we explain the
single-ion models used to fit the INS and EPR data. In Sec. IV, we present the
results and discuss them in the context of the single-ion model, from which
the conclusions are drawn in Sec. V.

##  II Experimental Methods

In order to investigate the large dynamic range expected for the Co^2+ single-
ion energy levels, we used the SEQUOIA time-of-flight chopper spectrometer in
the Spallation Neutron Source at Oak Ridge National Laboratory for the INS
experiment [25]. This instrument is optimized for magnetic neutron scattering
studies, with ample low-$|\mathbf{Q}|$ detector coverage, and allows a wide
range of incident neutron energies. The experiment was performed on an 8.62 g
powder sample of pure CoNb2O6, which we synthesized via a standard solid state
synthesis as reported elsewhere [26]. The powder was loaded into a cylindrical
aluminum sample can under helium exchange gas. Data were taken with a wide
range of incident energies ($E_{i}=60,150,250,700$ and
$2500\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$) at two temperatures where
appropriate ($T=$5\text{\,}\mathrm{K}$$ and $T=$200\text{\,}\mathrm{K}$$).
Instrumental settings and configurations, such as the $T_{0}$ and Fermi
chopper (FC1 and FC2) and elastic energy resolutions for each of the incident
neutron energies, can be found in the appendix.

We acquired the low-temperature ($T=$5\text{\,}\mathrm{K}$$) EPR spectra using
a Bruker EMX Elexsys spectrometer equipped with a ColdEdge closed-cycle helium
cryostat. We used a signal with microwave frequency
$f=$9.6453\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$ and a modulating field
amplitude of $A_{mod}=$0.4\text{\,}\mathrm{m}\mathrm{T}$$, corresponding to a
microwave power of $P=$6.325\text{\,}\mathrm{m}\mathrm{W}$$. For EPR
measurements, we substituted 1% Co into the isomorphic non-magnetic analog
MgNb2O6 system (Mg_0.99Co_0.01Nb2O6). This created magnetically isolated Co2+
ions in the the same crystal field environment as in CoNb2O6, ensuring that
spin-spin interactions did not interfere with the determination of the single-
ion $g$-values. The spectra were analyzed using the EasySpin toolbox for
MATLAB [27], where we simulated a spectrum that was calculated by defining a
system consisting of a pseudospin $(S_{\text{eff}}=1/2)$ (described in Sec.
III.3) and a single ^59Co nucleus, which has a natural abundance of 1.

##  III  Single-Ion Modeling

In this section we present the methods used in modeling the single-ion data
throughout this study, and were used to fit the data presented in Sec. IV. We
first introduce the intermediate spin-orbit coupled Hamiltonian used in
calculating the Co^2+ single-ion energy levels, as well as the resulting INS
response. This method, described in Sec. III.1, has been demonstrably
effective in $3d$ transition metal compounds, particularly in octahedrally-
coordinated Co^2+ complexes [28, 29, 30, 31] and is analogous to those
typically used in rare-earth magnets [32, 33]; the difference being in that
for the rare-earth case the single-ion Hamiltonian is formulated firmly in the
$J=||L+S||$ basis, rather than treating the spin-orbit coupling term on equal-
footing with the CEF term. We then describe the EPR Hamiltonian. The resulting
phenomenology is that the single-ion Co^2+ levels in a distorted octahedron
are Kramers doublets with magnetic dipole moments described by an anisotropic
g-tensor.

Figure 2: Energy level scheme depicting the splitting due to perturbations on
the ${}^{4}F$ free-ion ground state. Note that the left three panels in this
diagram corresponds to a purely cubic CEF (i.e., an octahedral field with no
distortions). The cubic CEF splits the free-ion ground state into three
orbital multiplets; the ground-state ${}^{4}T_{1}$ triplet is well-separated
in energy from the excited states [34]. The ground state degeneracy is reduced
further into three spin-orbit-split multiplets, which is further split into 14
doublets. The energy eigenvalues of the lowest doublets resultant from the INS
fits are displayed in the right panel (see Table 1 for a complete and specific
list of all energy eigenvalues).

###  III.1  Single-ion Hamiltonian

We begin by noting that for a $3d^{7}$ transition metal ion Hund’s rules
dictate that the Co^2+ free ion ground state is ${}^{4}F$ $(S=3/2,L=3)$,
leading to a $28\times 28$ Hilbert space. The odd number of unpaired electrons
($d^{7}$) implies all eigenstates are at least doubly degenerate, per Kramers
degeneracy theorem. As is common for $3d$ transition metals [32], we
approximate the single-ion Hamiltonian ($\hat{\mathcal{H}}_{\text{SI}}$) to
consist of a contribution from the crystal electric field and a contribution
from spin-orbit coupling (SOC) as follows,

$\hat{\mathcal{H}}_{\text{SI}}(\mathbf{L},\mathbf{S})=\hat{\mathcal{H}}_{\text{CEF}}(\mathbf{L})+\hat{\mathcal{H}}_{\text{SOC}}(\mathbf{L},\mathbf{S})$
(1)

where the CEF Hamiltonian $\hat{\mathcal{H}}_{\text{CEF}}(\mathbf{L})$ can be
written in Stevens form [32],

$\hat{\mathcal{H}}_{\text{CEF}}(\mathbf{L})=\sum_{l,m}B_{l}^{m}\hat{O}_{l}^{m}$
(2)

The $B_{l}^{m}\,(\hat{O}_{l}^{m})$ are the Stevens CEF parameters (operators).
For transition metal ions, the Stevens operators are written in terms of the
orbital angular momentum operators $L_{z},L_{+},L_{-}$ [35]. Point charge
calculations reveal that only eight CEF parameters are nonzero, due to the
$C_{2}$ point group symmetry of the Co^2+ average local environments. CEF
levels are inherently magnetic, and thus can be probed using INS. In Sec.
IV.1, we use INS data to determine the Steven’s parameters.

In the intermediate (Russell-Saunders) spin-orbit coupling scheme, the spin-
orbit interaction is given by

$\hat{\mathcal{H}}_{\text{SOC}}(\mathbf{L},\mathbf{S})=\lambda\mathbf{L}\cdot\mathbf{S}=\lambda\sum_{\alpha}\hat{L}_{\alpha}\hat{S}_{\alpha},$
(3)

where $\lambda=$-22.32\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ for the free
Co^2+ ion [34] and for $\alpha=(x,y)$, the orbital angular momentum operators
$\hat{L}_{\alpha}$ are expressed in terms of the usual raising and lowering
operators: $\hat{L}_{x}=\frac{1}{2}(\hat{L}_{+}+\hat{L}_{-})$ and
$\hat{L}_{y}=\frac{1}{2i}(\hat{L}_{+}-\hat{L}_{-})$, and likewise for
$\hat{S}_{x}$ and $\hat{S}_{y}$.

Diagonalizing the CEF Hamiltonian keeping only terms consistent with cubic
symmetry (i.e. a perfect octahedron) results in a ground state orbital triplet
$(^{4}T_{1})$ that is well-separated from the excited multiplets [29], as seen
in Fig. 2. From there, SOC and crystalline distortions further lower the
degeneracy resulting in a ground state doublet $(J_{\text{eff}}=\frac{1}{2})$.
In CoNb2O6, the ground state doublet is well isolated in energy (by
$E\approx$30\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$\approx$350\text{\,}\mathrm{K}$$),
as observed in the INS data presented in Sec. IV.1, allowing for the
description of this system as an effective spin-1/2 system at sufficiently low
temperatures. The components of the $g$-tensor can then be extracted in zero
field as the matrix elements of the magnetic dipole moment operator in the
subspace of the ground state wave functions $|\phi_{n}\rangle$ ($n=1,2$):

$g_{i}=2\langle\phi_{n}|\hat{L}_{i}+2\hat{S}_{i}|\phi_{n}\rangle,$

Strictly speaking, the “$g$-tensor” is not actually a tensor (it does not
transform properly under rotations), and these $g$-values are those obtained
from the square root of the eigenvalues of the $g$-tensor and correspondingly,
the principal axes are the eigenvectors of $G$ (see Ref. [32]).

###  III.2  Inelastic neutron scattering response

The single-ion dynamic structure factor (DSF) associated with Eqn. (1) at
constant momentum transfer $|\mathbf{Q}|$ is calculated in the following form
[33]:

$S(E)=C\sum_{n,m,\alpha}\frac{e^{-\beta
E_{n}}}{Z}\frac{\Gamma_{(n,m)}|\langle\phi_{n}|\hat{L}_{\alpha}+2\hat{S}_{\alpha}|\phi_{m}\rangle|^{2}}{(E_{n}-E_{m}-E)^{2}+\Gamma_{(n,m)}^{2}},$
(4)

where $\beta=1/k_{B}T$, $Z$ is the partition function, $|\phi_{n}\rangle$ is
the wave function associated with energy eigenvalue $E_{n}$ from the single-
ion calculations, and $\Gamma_{(n,m)}$ is a Lorentzian half-width at half-
maximum (HWHM) that can account for spectral broadening due to dispersion,
finite lifetimes in excited states, inhomogeneous broadening, and approximate
instrumental resolution. The actual observed intensity is related to the DSF
through
$I(|\mathbf{Q}|,E)=\frac{k_{f}}{k_{i}}f^{2}(|\mathbf{Q}|)S(|\mathbf{Q}|,E)$
[36], where $k_{f},k_{i}$ are the final and incident wave numbers of neutrons
that scatter with corresponding energy transfer $\Delta
E=\frac{\hbar^{2}}{2m}(k_{i}^{2}-k_{f}^{2})$ and $f(|\mathbf{Q}|)$ is the
magnetic form factor. At constant $|\mathbf{Q}|$, the form factor is also
constant and is assimilated into the constant $C$ in Eqn. (4), which acts as
an overall scale factor.

###  III.3  EPR Hamiltonian

The spin Hamiltonian used in EPR spectral analysis is approximated by two
contributions due to the electron Zeeman interaction and the hyperfine
coupling between the $3d^{7}$ electrons and the ^59Co nuclear spin ($I=7/2$).
It is important to emphasize that here these interactions are written in terms
of the effective spin $S_{\text{eff}}$ rather than the “bare” spin. For
systems with strong single-ion anisotropy, the interaction terms are written
in the following tensor form [37],

$\hat{\mathcal{H}}_{\text{spin}}(\mathbf{S}_{\text{eff}},\mathbf{I}\,)=-\mu_{B}{\mathbf{S}_{\text{eff}}\cdot\hat{g}\cdot{\mathbf{B}}+\hbar{\mathbf{S}_{\text{eff}}}\cdot\hat{A}\cdot{\mathbf{I}}}$
(5)

where $\mathbf{B}$ is the magnetic field responsible for driving the EPR
transition, and $\hat{g},\,\hat{A}$ are the $g$ and hyperfine coupling
tensors, respectively. The eigenvalues of these tensors are readily extracted
from fitting to EPR spectral data. For this work, the absorption lineshape is
assumed to be Lorentzian.

##  IV Results and Discussion

###  IV.1 Single-ion energy levels from INS

Figure 3:  Dynamic structure factor $S(|\mathbf{Q}|,E)$ data, measured via INS
on a CoNb2O6 powder sample. (a) Data at $T=5$ K using $E_{i}=60$ meV neutrons.
The absence of any low-lying CEF levels below
$E=$30\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ indicates the ground state is
a well-isolated doublet at temperatures well below those at which the first
excited state is thermally populated $(T\ll$350\text{\,}\mathrm{K}$)$. (b)
Data at $T=5$ K using $E_{i}=250$ meV neutrons, displaying the transition
energies of the first four excited doublets. Prominent phonon contributions
are seen at high $|\mathbf{Q}|$, spanning from the elastic line to just above
$E=100$ meV. (c) Constant-$|\mathbf{Q}|$ cut through the data in panel (b),
along with the calculated fit based on the single-ion model. Note that the fit
includes a peak at $E=105$ meV that is not visible in the data; the calculated
intensity is relatively small, so the peak is likely disguised by the phonon
background.

The single-ion levels of Co^2+ were observed as peaks corresponding to
transitions in the single-ion DSF acquired from INS experiments on the SEQUOIA
time-of-flight spectrometer. These were fit to a CEF Hamiltonian, and the
resulting calculated energy levels ranged from 0 to
$902\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$. Only a few of these modes were
observed via INS due to weak CEF level intensities and large background at the
higher energy transfers; a complete list of all calculated and observed
energies is found in Table 1, along with the calculated intensities
(normalized to the $30\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ peak) for each
of the modes. An asterisk denotes the mode was not observed, while a dagger
indicates that there were some weak magnetic signatures in constant-$E$ cuts
(see Fig. 6). The modes and their intensities that were observed, along with
the $g$-values from EPR, provide sufficient information to obtain the Steven’s
parameters from fits. The DSF at $T=$5 K for incident neutron energies
$E_{i}=60,\,$250\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ is presented in Fig.
3, along with an intensity versus energy transfer cut and fit for
$E_{i}=$250\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$. All data presented in
this article are background subtracted to remove the contributions from the
sample environment, including the aluminum sample can. These “empty-can”
datasets were collected at all specified incident neutron energies and sample
environment temperatures. The remaining background results from sample-
specific features, such as coherent and incoherent phonon scattering. All data
presented here were corrected by a factor of $k_{f}/k_{i}$ which means that
the quantity being displayed in Fig. 3 is proportional to
$f^{2}(|\mathbf{Q}|)S(|\mathbf{Q}|,E)$ (in arbitrary units).

Table 1:  Calculated single-ion Co^2+ energy levels in CoNb2O6 from the fit to the INS data shown in Fig. 3 compared to observed approximate energy levels. Note that not all of the transitions listed were observed in the data, as is further discussed in the main text. Calculated and observed intensities (relative to the $30\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ peak) are displayed in the far right column, and are separated by horizontal lines that indicate incident neutron energies used in the calculation ($E_{i}=250,700,$2500\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$, from top to bottom). Unobserved modes have an asterisk but a dagger indicates there were weak magnetic signatures in constant-$E$ cuts (Fig. 6 in App. B). Eigenstate | $E_{m}$ | Obs. | Calc. Int. | Obs. Int.
---|---|---|---|---
$\phi_{1},\phi_{2}$ | $0$ | 0 | - | -
$\phi_{3},\phi_{4}$ | $30.4$ | $29.8(8)$ | 1 | 1
$\phi_{5},\phi_{6}$ | $51.9$ | $49(2)$ | 0.30 | 0.0181
$\phi_{7},\phi_{8}$ | $104$ | $*\dagger$ | 0.024 | *
$\phi_{9},\phi_{10}$ | $133$ | $133(1)$ | 0.076 | 0.0039
$\phi_{11},\phi_{12}$ | $150$ | * | 0.0007 | *
$\phi_{13},\phi_{14}$ | $377$ | $*$ | 0.060 | *
$\phi_{15},\phi_{16}$ | $390$ | $*$ | 0.056 | *
$\phi_{17},\phi_{18}$ | $467$ | $*$ | 0.052 | *
$\phi_{19},\phi_{20}$ | $477$ | $*$ | 0.052 | *
$\phi_{21},\phi_{22}$ | $493$ | $*$ | 0.041 | *
$\phi_{23},\phi_{24}$ | $511$ | $*$ | 0.027 | *
$\phi_{25},\phi_{26}$ | $898$ | $*\dagger$ | 0.040 | *
$\phi_{27},\phi_{28}$ | $902$ | $*\dagger$ | 0.039 | *

Our data analysis consisted of fitting the single-ion model described in Sec.
III.2 to the $E_{i}=$250\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ data, which
showed the clearest features. A Gaussian peak at the elastic line with full-
width at half-maximum (FWHM) taken from the instrument resolution function was
also added to account for both coherent and incoherent nuclear elastic
scattering. Then, the instrument resolution function was further used to
provided the FWHM in fitting the inelastic peaks to a Lorentzian profile.
Overall, nine parameters were fit, including the eight CEF parameters
$(B_{l}^{m})$ and the overall scale factor $C$. The data contained 7 clear
pieces of information; 3 peak positions and 4 peak intensities (the intensity
of the peak near 100 meV is an upper limit since the position of said peak is
difficult to discern); thus the fit would be underconstrained. To remedy this,
we added a penalty for the distance of the calculated principal $g$-values
from the INS model, to those measured via EPR (reported in Sec. IV.3). Fits
were accomplished by first employing a particle swarm optimization method
using the point charge calculation as a staring point, then once the algorithm
found a “global” minimum, we used nonlinear least-squares to explore the
region around the minimum and estimate the uncertainties of the parameters. As
seen in Fig. 3(c), there is a broad signal underlying the CEF peaks that we
did decided not to attempt to model phenomenologically, since it still
represents constraint on peak positions and intensities. This “background” may
have contributed to the optimizer’s tendency to get stuck in local minima
without identifying the strong peaks in the data. To improve the tendency to
find the peaks, we added an additional penalty for the distance between the
nearest calculated transitions and the observed peak positions. This was
accomplished by first fitting the peaks individually to extract a peak
position, and then incorporating a least-squares term to match the calculated
peak centers to the observed (this is partially accomplished by fitting the
full intensity profile, but this gives it an extra boost for matching peak
positions). A similar procedure used in an example fit program that is
provided by the software package PyCrystalField [38]. The final results of the
fitted parameters (except for the scale factor) are presented in Table 2. The
first excited level was detected using an incident neutron energy of
$E_{i}=$60\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ was at
$E=$30\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$, and is typical for
octahedrally-coordinated Co^2+ compounds [28, 30, 31]. In such systems it is
attributable to the transition between $J_{\text{eff}}=1/2$ and $3/2$
multiplets [30, 39]. The level was observed at both temperatures
$(T=5,\,$200\text{\,}\mathrm{K}$)$, and displayed a decrease in intensity as
momentum transfer $|\mathbf{Q}|$ increases, which is characteristic of
magnetic excitations.

Table 2: Stevens parameters extracted from fitting to the inelastic neutron scattering data (the Co^2+ single-ion DSF, Eqn. (4)).  Parameter |  | Value$($\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$)$ |
---|---|---|---
$B_{2,0}$ |  | $-1.14(6)$ |
$B_{2,1}$ |  | $0.1(2)$ |
$B_{2,2}$ |  | $6(2)$ |
$B_{4,0}$ |  | $-0.78(4)$ |
$B_{4,1}$ |  | $-0.55(5)$ |
$B_{4,2}$ |  | $0.15(5)$ |
$B_{4,3}$ |  | $0.330(1)$ |
$B_{4,4}$ |  | $3.29(4)$ |

Some spectral features calculated from the fit deviated from the observed
modes in various ways. For example, the mode at
$E=$51.9\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ is observed to be wider and
much less intense than calculated but is appropriately positioned. We
attribute this to possible structural disorder, as has been detailed by Sarkis
et al. in recent AC susceptibility work on CoNb2O6 [40]. Such inhomogeneities
would lead to a range of CEF levels and would broaden observed CEF
transitions.

The $51.9\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ peak corresponds to
transitions only between the ground state and excited states, as the first
excited level is not appreciably thermally populated, since it lies at energy
$E=$348\text{\,}\mathrm{K}$/k_{B}$ above the ground state and orbital
selection rules forbid intra-multiplet transitions. This however does not
exclude the presence of mixing between CEF and exchange interactions, as the
dominant FM exchange energy
($J_{0}\approx$3\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ [19]) is only about
an order of magnitude smaller than the CEF splittings and can possibly lead to
dispersion that can widen the signal. In this case (which is seen in both
Co^2+systems studied in Ref. [28]), these dispersive effects would be
invisible to us since the energies are smaller than the SEQUOIA instrument
resolution for the incident neutron energies used (found in App. C).

There is then the absence of a clear peak near
$E=$104\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ which is suggested by the
single-ion calculation. Looking with a loving eye there is perhaps a slight
peak visible there, and signatures of a magnetic excitation were seen in
constant-$E$ cuts of the
$E_{i}=$150\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$,T=5,$200\text{\,}\mathrm{K}$$
data; the same is true for the highest energy modes (see App. B). While the
masking of lower-energy excitations could potentially be simply explained as
phonon background, the higher-energy cases are not so simple. In this regime
the modes could possibly be concealed by more complex background signals such
as those found due to multi-phonon processes and recoil scattering. This
observed recoil scattering (see Fig. 6 in Appendix B) is plausibly attributed
to the presence of hydrogen, because sharp peaks were found in the
$E_{i}=$2500\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ data at energies that
are consistent with the various O-H bending and stretching modes displayed in
water molecules and hydroxyl groups and suggests that some water had
accumulated on the sample’s surface.

###  IV.2 Calculated CEF ground state and $g$-tensor

The spin and orbital composition of the CEF ground state doublet
$|\phi_{1,2}\rangle$ obtained from the analysis of the INS data is listed in
Table 4 located in Appendix A. The largest contributions come from the
$L=0,1,3$ subspaces. The approximate wavefunctions for the ground state
doublet in the $|L_{z},S_{z}\rangle$ basis are:

$|\phi_{1,2}\rangle\approx\mp 0.63|\pm 3,\pm\frac{3}{2}\rangle\pm
0.53|0,\pm\frac{1}{2}\rangle\pm 0.38|\mp 1,\pm\frac{3}{2}\rangle\pm...$ (6)

The unabbreviated form of Eqn.(6) is what is used to calculate the $g$-tensor
of the system from the INS results, the principal values of which are located
in Table 3. After INS data were fit, the obtained Stevens parameters were
adjusted within their respective 95% confidence intervals to find how much the
resulting $g$-values vary, the range of which is roughly equivalent to the 95%
CI (2$\sigma$) of the $g$-values. From the ground state wavefunction we can
also find the principal axes of the $g$-tensor relative to the
crystallographic directions, and find that the strongest value ($g_{z}=6.4$)
is along an axis in the ac-plane, approximately $37^{\circ}$ from c along a.

Table 3: Parameters from Eqn. (5) extracted by fitting to EPR signal acquired from the dilute sample and include the principal $g$-values ($g_{\alpha}$, unitless and located above the line), the hyperfine coupling values $(A_{\alpha})$, and the Lorentzian linewidth $\gamma$. For the $g$-values, the values obtained via EPR are displayed next to those calculated from the fit to the INS data (and their reasonable range) as described in Sec. LABEL:sec:SIG. The furthest right column shows the percent difference between the two methods. Parameter | EPR | INS | % Difference
---|---|---|---
$g_{x}$ | 3.33 | 3.9(6) | 16%
$g_{y}$ | 3.01 | 3.3(8) | 9.2%
$g_{z}$ | 6.39 | 6.3(6) | 1.4%
$A_{x}$ | $123\text{\,}\mu\mathrm{e}\mathrm{V}$ |  |
$A_{y}$ | $71.5\text{\,}\mu\mathrm{e}\mathrm{V}$ |  |
$A_{z}$ | $521\text{\,}\mu\mathrm{e}\mathrm{V}$ |  |
$\gamma$ | $5.3\text{\,}\mathrm{m}\mathrm{T}$ |  |

This is reasonably close to the observed ordered moment direction of
$31^{\circ}$ in the ac-plane [41, 10]. We note that the observed ordered
moment direction is not necessarily required to be the same as the $g$-tensor
principal axes; ultimately, it is a combination of the exchange anisotropy and
the $g$-tensor anisotropy that dictates the moment direction; one clear
example of this is the splayed ferromagnetism found in systems like Yb2Ti2O7
and Yb2Sn2O7 [42, 43]. As for the final results of the fitted Stevens
parameters, it was found that the confidence interval for $B_{2}^{1}$
straddled the value of zero, and setting the parameter equal to zero had
negligible effect on the $g$-tensor and the calculated intensity pattern. In
fact, when $B_{2}^{1}$ is entirely removed from the fitting routine it results
in the same energy spectrum, but slightly changes the $g$-values (all three
differ from the reported results by less than 6%) and rotates the principal
axes by $3^{\circ}$ in the ac-plane, as well as reducing the error on the
parameters by at least an order of magnitude. This could all indicate that
$B_{2}^{1}$ is an extraneous parameter that is not necessary to describe this
system. We stress that this does not necessarily imply a higher point group
symmetry; the next highest symmetry ($C_{2v}$) would require additional
parameters to vanish.

###  IV.3 EPR spectral analysis

The low-temperature X-band EPR spectrum is displayed in Fig. 4. This
measurement was done on a powder sample, so the directions of the principal
axes of the $\hat{g}$ and $\hat{A}$ could not be directly measured, but the
principal values $(g_{ii},A_{ii})$ extracted from the fit are displayed in
Table 3. The existence of 3 separate principal $g$-values is consistent with a
rhombic $g$-tensor. This is in contrast to an axial system
$(g_{z}=g_{\parallel},\,g_{x}=g_{y}=g_{\perp})$ which is usually associated
with Ising-like single-ion systems. Two of the $g$-values are close but not
equal, leading to an overlap of spectral features associated with the values.
There is a noticeable difference in the clarity of the hyperfine-split
features between the experiment and simulated signals, which is likely due to
insufficient instrument resolution. Microwave power and modulation amplitudes
were adjusted in an effort to resolve the splittings. Note the two features at
$B=150$ and $275\text{\,}\mathrm{m}\mathrm{T}$ in the experimental data which
are not accounted for in the simulation and the source(s) remain unkown.

The principal $g$-values from the EPR spectrum follow the same general trend
as those determined from the CEF calculations ($g_{z}>g_{y}>g_{x}$), and the
values for $g_{z}$ (associated with the Ising moment direction) differ by
approximately $0.6\%$. The reason for the larger differences between the
values of $g_{x}$ and $g_{y}$ is likely compositional, i.e., the INS
measurements were performed on pure CoNb2O6, and EPR spectra were acquired
from the extremely dilute sample. However, we find it important to make a
certain point more salient: the "most correct" $g$-values are likely the ones
found via EPR, as the dilute sample provided an environment that would be much
closer to true single-ion behavior. However, it is also true that the sample
is primarily MgNb2O6, which will have slightly different structural parameters
than CoNb2O6, so these g-values are still an approximation to the “true”
single-ion behavior of CoNb2O6.

As mentioned previously, the EPR spectrum was approximated using a
pseudospin-1/2 model, which means the zero-field splitting (ZFS) terms usually
present in $S>1/2$ systems (like Co^2+ when the full $S=3/2$ $L=3$ manifold is
included [37]) were omitted from Eqn. (5). Additionally, the ^59Co nucleus has
a large quadrupolar moment $(Q/|e|=0.42\times
10^{-24}\,$\mathrm{c}\mathrm{m}^{-2}$$), so a quadrupole interaction term
should also be included. However, the absence of both ZFS and quadrupole terms
had negligible effect on the simulation. This is consistent with another
octahedrally-coordinated Co^2+ compound, discussed in the Supplementary
Information of Ref. [44]. Hyperfine interaction strengths are
$\sim$\text{\,}\mu\mathrm{e}\mathrm{V}$$ and therefore cannot be corroborated
with our INS data due to having insufficient energy resolution. A Lorentzian
FWHM of $\gamma=$5.3\text{\,}\mathrm{m}\mathrm{T}$$ accounted for the
linewidth due to inhomogeneous broadening from unresolved hyperfine features.

Figure 4: Field-swept X-band EPR response at $T=5$ K, compared to a simulated
system. The g-values can be estimated from the center positions of the major
features. The spacing between the jagged peaks correspond to hyperfine
splittings. The cause of the spectral features at approximately
$B=$150.275\text{\,}\mathrm{m}\mathrm{T}$$ (indicated by blue arrows) remain
unknown.

##  V Conclusions

We have presented inelastic neutron scattering and electron paramagnetic
resonance results on the Ising material CoNb2O6. The main results are the
measurements and analysis of the Co^2+ single-ion properties of this system;
primarily the single-ion energy levels, the ground state CEF wavefunction, and
the $g$-values. The ground state of Co^2+ in these octahedrally-coordinated
environments is a Kramers doublet which is well-separated from the first
excited states by over $300\text{\,}\mathrm{K}$. These results show that the
system is well-described by an intermediate spin-orbit coupled Hamiltonian.
Further calculations are needed to explore the Kitaev-like nature of the
interactions in CoNb2O6, which can be accomplished via similar methods to
those used in Ref.[22] and using the CEF ground state wavefunctions and
principal $g$-values that we have provided.

###### Acknowledgements.

This research was supported by the U.S. Department of Energy, grant no. DE-
SC0018972. The authors would like to gratefully acknowledge the discussions
with Joe Zadrozny and Richard Miller. The neutron scattering data were reduced
using the Mantid software package [45], while data analysis utilized
PyCrystalField and MATLAB [46, 38]. The authors also acknowledge the Molecular
and Materials Analysis Center of the Analytical Resources Core at Colorado
State University for instrument access and training. A portion of this
research at Oak Ridge National Laboratory’s Spallation Neutron Source was
sponsored by the Scientific User Facilities Division, Office of Basic Energy
Sciences, U.S. Department of Energy.

##  Appendix A  CEF ground state wavefunctions and g-tensor orientation

The wavefunctions for Co^2+ are expressed as linear combinations of
eigenstates in the intermediate coupling $(\lvert L_{z},S_{z}\rangle)$ basis
as follows:

$\lvert\phi_{n}\rangle=\sum_{L_{z},S_{z}}C^{(n)}_{L_{z},S_{z}}\lvert\,L_{z},S_{z}\rangle$
(7)

Where $L_{z},S_{z}$ are the eigenvalues of the orbital and spin angular
momentum operators $\hat{L}_{z}$ and $\hat{S}_{z}$, respectively. The
$C^{(n)}_{l,s}$ are the amplitudes of the corresponding eigenstates, which can
be found in Table 4. A plot of the absolute square amplitudes is shown in Fig.
5, and demonstrates the highly symmetric nature of the ground state. Knowing
the composition of the ground state wavefunctions in the pseudospin-1/2 regime
is one key to understanding the possible Kitaev nature of the microscopics in
coordinated Co^2+ complexes.

Table 4: List of amplitudes $C^{(n)}$ corresponding to the single-ion Co^2+ ground state doublet wavefunctions $(n=1,2)$ in the $\lvert L_{z},S_{z}\rangle$ basis. See Eqn. 7. Lefthand column is the index of the contributions from the eigenstates of the $|L_{z},S_{z}\rangle$ basis, and is important for understanding Fig. 5. Index | $L_{z}$ | $S_{z}$ | $C_{L_{z},S_{z}}^{(1)}$ | $C_{L_{z},S_{z}}^{(2)}$
---|---|---|---|---
1 | -3 | -3/2 | 0.0498 | 0.6290
2 | -3 | -1/2 | -0.2516 | -0.1213
3 | -3 | 1/2 | 0.0185 | 0.0765
4 | -3 | 3/2 | -0.0841 | -0.0048
5 | -2 | -3/2 | -0.0324 | -0.0483
6 | -2 | -1/2 | 0.0165 | 0.0234
7 | -2 | 1/2 | -0.0161 | -0.0168
8 | -2 | 3/2 | 0.0278 | -0.0255
9 | -1 | -3/2 | -0.0012 | -0.0764
10 | -1 | -1/2 | 0.0497 | -0.0093
11 | -1 | 1/2 | 0.0731 | -0.2433
12 | -1 | 3/2 | 0.3791 | 0.0044
13 | 0 | -3/2 | 0.0458 | 0.0878
14 | 0 | -1/2 | -0.0016 | -0.5293
15 | 0 | 1/2 | 0.5308 | 0.0404
16 | 0 | 3/2 | 0.0844 | -0.0390
17 | 1 | -3/2 | -0.0257 | -0.3799
18 | 1 | -1/2 | 0.2499 | 0.0927
19 | 1 | 1/2 | -0.0132 | -0.0506
20 | 1 | 3/2 | 0.0766 | 0.0048
21 | 2 | -3/2 | 0.0278 | 0.0299
22 | 2 | -1/2 | -0.0156 | 0.0148
23 | 2 | 1/2 | -0.0221 | 0.0147
24 | 2 | 3/2 | -0.0459 | 0.0286
25 | 3 | -3/2 | 0.0019 | 0.0840
26 | 3 | -1/2 | -0.0752 | 0.0125
27 | 3 | 1/2 | -0.1017 | 0.2428
28 | 3 | 3/2 | -0.6270 | -0.0000

Furthermore, it is important to explain the coordinate system regarding the
$g$-tensor principal axes as calculated from the single-ion Hamiltonian (Eqn.
1). As stated in the main text, the fitting routine used starting parameters
from a point charge calculation in which the coordinate axes were aligned
along the crystallographic axes, as in, ($x\parallel a,y\parallel b,z\parallel
c$). This means that the final fit parameters are also expressed in that
coordinate system, where the $g$-tensor is not diagonal but rotated an angle
$\theta=37^{\circ}$ about the b-axis, which is consistent with the observed
moment direction but differs by a few degrees (See Fig. 5 below). We wish to
reiterate that the observed moment direction is a combination of both
$g$-tensor and exchange anisotropy.

Figure 5: (a) A plot of the absolute square amplitudes of the CEF ground state
doublet. The horizontal axis is made up of integers that correspond to the
indices of contributions from the eigenstates of the $|L_{z},S_{z}\rangle$
basis as found in Table 4. The wavefunctions are highly symmetric. Stars
indicate the greatest contributions. (b) Depiction of the $g$-tensor principal
axes (black vectors) as calculated by the methods described in Sec. III.1
referenced to the crystallographic directions, where two axes are rotated
approximately $\theta=37^{\circ}$ from the c-axis in the ac-plane. The third
principal axis is constrained by symmetry to lie along b.

##  Appendix B Signatures of unobserved magnetic modes and recoil scattering

As mentioned in the main text, not all modes were observed. Modes unobserved
in constant-$|\bf{Q}|$ cuts such as the one shown in Fig. 3 can be identified
as being magnetic by observing their $|\bf{Q}|$-dependence at constant energy
transfer; the magnetic form factor $f(|\bf{Q}|)$ causes the intensity to drop
off with increasing momentum transfer. This is seen in Fig. 6(a,b,d), which
displays constant-$E$ cuts of the
$E_{i}=150,$2500\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ data to expose the
magnetic origins of the $104,898$, and
$901\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ modes.

Figure 6: Magnetic origins of several modes unobserved in constant-$|\bf{Q}|$
cuts. (a) Constant-$E$ cuts of the data with incident neutron energy
$E_{i}=$150\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ at
$T=$5\text{\,}\mathrm{K}$$, demonstrating the magnetic nature of the
$104\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$ mode. (b) Similar cuts for the
same incident neutron energy but at $T=$200\text{\,}\mathrm{K}$$,
corroborating the results in (a). (c) Slice of the
$E_{i}=$2500\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$$ data taken at
$T=$5\text{\,}\mathrm{K}$$. The arrow indicates where strong recoil scattering
is observed (as mentioned in Sec. IV.1), and likely masks much of the
intermediate modes. (d) Cuts of the data in panel (c) showing the magnetic
signatures of the excitations near
$$900\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$.$

##  Appendix C Instrumental configurations used on SEQUOIA

A list of configurational settings for the SEQUOIA spectrometer used in this
research is listed in Table 5.

Table 5: Instrumental configuration settings for the SEQOUIA spectrometer for all incident neutron energies $(E_{i})$ and sample environment temperatures. Includes the frequency settings for the $T_{0}$ and high flux/resolution Fermi choppers (FC1/FC2), as well as the elastic energy resolution $\delta E$. $E_{i}($\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$)$ | $T($\text{\,}\mathrm{K}$)$ | $T_{0}($\text{\,}\mathrm{H}\mathrm{z}$)$ | FC $($\text{\,}\mathrm{H}\mathrm{z}$)$ | $\delta E($\text{\,}\mathrm{m}\mathrm{e}\mathrm{V}$)$
---|---|---|---|---
60 | 5, 200 | 90 | FC2, 420 | 1.61
150 | 5, 200 | 90 | FC2, 600 | 3.98
250 | 5, 200 | 120 | FC2, 600 | 6.95
700 | 5 | 150 | FC1, 600 | 50.73
2500 | 5 | 180 | FC1, 600 | 334.75

## References

  * Dutta  _et al._ [2015] A. Dutta, G. Aeppli, B. K. Chakrabarti, U. Divakaran, T. F. Rosenbaum, and D. Sen,  _Quantum Phase Transitions in Transverse Field Spin Models: From Statistical Physics to Quantum Information_ (Cambridge University Press, 2015).
  * Oshiyama  _et al._ [2020] H. Oshiyama, N. Shibata, and S. Suzuki, J. Phys. Soc. Jpn. 89, 104002 (2020).
  * Silevitch  _et al._ [2019] D. M. Silevitch, C. Tang, G. Aeppli, and T. F. Rosenbaum, Nat Commun 10, 1 (2019).
  * Kjäll  _et al._ [2014] J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Phys. Rev. Lett. 113, 107204 (2014).
  * Gingras and Henelius [2011] M. J. P. Gingras and P. Henelius, J. Phys.: Conf. Ser. 320, 012001 (2011).
  * Cui  _et al._ [2019] Y. Cui, H. Zou, N. Xi, Z. He, Y. X. Yang, L. Shu, G. H. Zhang, Z. Hu, T. Chen, R. Yu, J. Wu, and W. Yu, Phys. Rev. Lett. 123, 067203 (2019).
  * Coldea  _et al._ [2010] R. Coldea, D. A. Tennant, E. M. Wheeler, E. Wawrzynska, D. Prabhakaran, M. Telling, K. Habicht, P. Smeibidl, and K. Kiefer, Science 327, 177 (2010).
  * Note [1] It should be noted that the 3D transition actually occurs at $\mu_{0}H_{c}=$5.5\text{\,}\mathrm{T}$$, as observed by inelastic neutron scattering measurements [7]. However, the 1D QCP at zero temperature can be inferred from properties of the quantum critical fan [12].
  * Kobayashi  _et al._ [2014] S. Kobayashi, S. Hosaka, H. Tamatsukuri, T. Nakajima, S. Mitsuda, K. Prokeš, and K. Kiefer, Phys. Rev. B 90, 060412(R) (2014).
  * Kobayashi  _et al._ [2016] S. Kobayashi, S. Mitsuda, S. Hosaka, H. Tamatsukuri, T. Nakajima, H. Koorikawa, K. Prokeš, and K. Kiefer, Phys. Rev. B 94, 134427 (2016).
  * Liang  _et al._ [2015] T. Liang, S. M. Koohpayeh, J. W. Krizan, T. M. McQueen, R. J. Cava, and N. P. Ong, Nat. Commun. 6, 7611 (2015).
  * Kinross  _et al._ [2014] A. W. Kinross, M. Fu, T. J. Munsie, H. A. Dabkowska, G. M. Luke, S. Sachdev, and T. Imai, Phys. Rev. X 4, 031008 (2014).
  * Munsie [2017] T. Munsie,  _Synthesis and Characterization of Constrained Magnetism in Niobates_ , PhD Thesis, McMaster University (2017).
  * Pullar [2009] R. C. Pullar, J. Am. Ceram. Soc. 92, 563 (2009).
  * Maartense  _et al._ [1977] I. Maartense, I. Yaeger, and B. M. Wanklyn, Solid State Communications 21, 93 (1977).
  * Heid  _et al._ [1995] C. Heid, H. Weitzel, P. Burlet, M. Bonnet, W. Gonschorek, T. Vogt, J. Norwig, and H. Fuess, Journal of Magnetism and Magnetic Materials 151, 123 (1995).
  * Kobayashi  _et al._ [2000] S. Kobayashi, S. Mitsuda, and K. Prokes, Phys. Rev. B 63, 024415 (2000).
  * Robinson  _et al._ [2014] N. J. Robinson, F. H. L. Essler, I. Cabrera, and R. Coldea, Phys. Rev. B 90, 174406 (2014).
  * Fava  _et al._ [2020] M. Fava, R. Coldea, and S. A. Parameswaran, Proceedings of the National Academy of Sciences of the United States of America 117, 25219 (2020).
  * Morris  _et al._ [2021] C. M. Morris, N. Desai, J. Viirok, D. Huvonen, U. Nagel, T. Room, J. W. Krizan, R. J. Cava, T. M. McQueen, S. M. Koohpayeh, R. K. Kaul, and N. P. Armitage, Nature Phys. 17, 832 (2021).
  * Kramers and Wannier [1941] H. A. Kramers and G. Wannier, Phys. Rev. 60, 252 (1941).
  * Liu  _et al._ [2020] H. Liu, J. Chaloupka, and G. Khaliullin, Physical Review Letters 125, 047201 (2020).
  * Vivanco  _et al._ [2020] H. K. Vivanco, B. A. Trump, C. M. Brown, and T. M. McQueen, Phys. Rev. B 102, 224411 (2020).
  * Zhong  _et al._ [2020] R. Zhong, T. Gao, N. P. Ong, and R. J. Cava, Sci. Adv. 6, eaay6953 (2020).
  * Granroth  _et al._ [2010] G. E. Granroth, A. I. Kolesnikov, T. E. Sherline, J. P. Clancy, K. A. Ross, J. P. Ruff, B. D. Gaulin, and S. E. Nagler, Journal of Physics: Conference Series 251, 12058 (2010).
  * Sarvezuk  _et al._ [2011] P. W. C. Sarvezuk, E. J. Kinast, C. V. Colin, M. A. Gusmão, J. B. M. da Cunha, and O. Isnard, Journal of Applied Physics 109, 07E160 (2011).
  * Stoll and Schweiger [2006] S. Stoll and A. Schweiger, J. Magn. Reson. 178(1), 42 (2006).
  * Ross  _et al._ [2017] K. A. Ross, J. M. Brown, R. J. Cava, J. W. Krizan, S. E. Nagler, J. A. Rodriguez-Rivera, and M. B. Stone, Phys. Rev. B 95, 144414 (2017).
  * Sarte  _et al._ [2019] P. M. Sarte, M. Songvilay, E. Pachoud, R. A. Ewings, C. D. Frost, D. Prabhakaran, K. H. Hong, A. J. Browne, Z. Yamani, J. P. Attfield, E. E. Rodriguez, S. D. Wilson, and C. Stock, Phys. Rev. B 100, 075143 (2019).
  * Tomiyasu  _et al._ [2011] K. Tomiyasu, M. K. Crawford, D. T. Adroja, P. Manuel, A. Tominaga, S. Hara, H. Sato, T. Watanabe, S. I. Ikeda, J. W. Lynn, K. Iwasa, and K. Yamada, Phys. Rev. B 84, 054405 (2011).
  * Holden  _et al._ [1971] T. M. Holden, W. J. L. Buyers, E. C. Svensson, R. A. Cowley, M. T. Hutchings, D. Hukin, and R. W. H. Stevenson, J. Phys. C: Solid State Phys. 4, 2127 (1971).
  * Abragam and Bleaney [2012] A. Abragam and B. Bleaney,  _Electron Paramagnetic Resonance of Transition Ions_ (OUP Oxford, 2012).
  * Jensen and Mackintosh [1991] J. Jensen and A. Mackintosh,  _Rare Earth Magnetism_ (Clarendon Press, Oxford, UK, 1991).
  * Kanamori [1957] J. Kanamori, Prog. Theor. Phys. 17, 177 (1957).
  * Hutchings [1964] M. T. Hutchings, Solid State Phys. 16, 224 (1964).
  * Squires [1996] G. L. Squires,  _Introduction to the Theory of Thermal Neutron Scattering_ (Dover Publications, 1996).
  * Piwowarska  _et al._ [2019] D. Piwowarska, P. Gnutek, and C. Rudowicz, Appl. Magn. Reson. 50, 797 (2019).
  * Scheie [2018] A. Scheie, Pycrystalfield, https://github.com/asche1/PyCrystalField (2018).
  * Buyers  _et al._ [1971] W. J. L. Buyers, T. M. Holden, E. C. Svensson, R. A. Cowley, and M. T. Hutchings, J. Phys. C: Solid State Phys. 4, 2139 (1971).
  * Sarkis  _et al._ [2021] C. L. Sarkis, S. Säubert, V. Williams, E.S. Choi, T.R. Reeder, H.S. Nair, and K.A. Ross, Physical Review B 104, 214424 (2021).
  * Scharf  _et al._ [1979] W. Scharf, H. Weitzel, I. Yaeger, I. Maartense, and B. M. Wanklyn, Journal of Magnetism and Magnetic Materials 13, 121 (1979).
  * Gaudet  _et al._ [2016] J. Gaudet, K. A. Ross, E. Kermarrec, N. P. Butch, G. Ehlers, H. A. Dabkowska, and B. D. Gaulin, Phys. Rev. B 93, 064406 (2016).
  * Yaouanc  _et al._ [2013] A. Yaouanc, P. Dalmas de Réotier, P. Bonville, J. A. Hodges, V. Glazkov, L. Keller, V. Sikolenko, M. Bartkowiak, A. Amato, C. Baines, P. J. C. King, P. C. M. Gubbens, and A. Forget, Phys. Rev. Lett. 110, 127207 (2013).
  * Gómez-Coca  _et al._ [2014] S. Gómez-Coca, A. Urtizberea, E. Cremades, P. J. Alonso, A. Camón, E. Ruiz, and F. Luis, Nat. Commun. 5, 4300 (2014).
  * Arnold  _et al._ [2014] O. Arnold, J. C. Bilheux, J. M. Borreguero, A. Buts, S. I. Campbell, L. Chapon, M. Doucet, N. Draper, R. Ferraz Leal, M. A. Gigg, V. E. Lynch, A. Markvardsen, D. J. Mikkelson, R. L. Mikkelson, R. Miller, K. Palmen, P. Parker, G. Passos, T. G. Perring, P. F. Peterson, S. Ren, M. A. Reuter, A. T. Savici, J. W. Taylor, R. J. Taylor, R. Tolchenov, W. Zhou, and J. Zikovsky, NIM 764, 156 (2014).
  * MATLAB [2010] MATLAB,  _version 7.10.0 (R2010a)_ (The MathWorks Inc., Natick, Massachusetts, 2010).