The Curious Case of NiRh$_2$O$_4$: A Spin-Orbit Entangled Diamond Lattice Paramagnet

TL;DR

This paper investigates NiRh$_2$O$_4$, revealing it as a spin-orbit entangled diamond lattice paramagnet with frustrated exchange interactions, using ab initio and model Hamiltonian approaches to explain its magnetic properties and excitations.

Contribution

It introduces a comprehensive theoretical framework combining ab initio calculations and model Hamiltonians to explain the magnetic behavior of NiRh$_2$O$_4$, emphasizing spin-orbit coupling and exchange frustration.

Findings

Identification of a spin-orbit entangled non-magnetic ground state.

Prediction of dispersive gapped magnetic modes and dark states.

Demonstration that exchange frustration suppresses exciton condensation.

Abstract

Motivated by the interest in topological quantum paramagnets in candidate spin-$1$ magnets, we investigate the diamond lattice compound NiRh$_2$O$_4$ using {\it ab initio} theory and model Hamiltonian approaches. Our density functional study, taking into account the unquenched orbital degrees of freedom, shows stabilization of $S=1$ and $L=1$ state. We highlight the importance of spin-orbit coupling, in addition to Coulomb correlations, in driving the insulating gap, and uncover frustrating large second-neighbor exchange mediated by Ni-Rh covalency. A single-site model Hamiltonian incorporating the large tetragonal distortion is shown to give rise to a spin-orbit entangled non-magnetic ground state, largely accounting for the entropy, magnetic susceptibility, and inelastic neutron scattering results. Incorporating inter-site exchange within a slave-boson theory, we show that exchange frustration can suppress exciton condensation. We capture the dispersive gapped magnetic modes, uncover `dark states' invisible to neutrons, and make predictions for future experiments.

Equations

$$H=-\delta L_z^2+\lambda \vec{L}\cdot \vec{S}$$

$$E_0^0[2]$$

$$E_2^0[2]$$

$${\Delta }_0=0$$

$${\Delta }_1\approx 2\frac{{\lambda }^2}{\delta }$$

$${\Delta }_2\approx \lambda +\frac{{\lambda }^2}{\delta }$$

$${\Delta }_3\approx 2\lambda +2\frac{{\lambda }^2}{\delta }$$

$${\Delta }_4\approx \delta +\lambda +3\frac{{\lambda }^2}{\delta }$$

$${\Delta }_5\approx \delta +\lambda +4\frac{{\lambda }^2}{\delta }$$

$$H_{\mathrm{loc}}$$

$$H_{\mathrm{ex}}$$

$$H_{\mathrm{con}}$$

$$H_{\mathrm{tot}}$$

$$E_{\mathbf{k}}^{\sigma }$$

$${\tilde{E}}_{\mathbf{k}}^{\sigma }$$

Full text

*Introduction. — * Symmetry protected topological phases of quantum matter, e.g., two dimensional (2D) and 3D topological insulators Hasan and Kane (2010); Qi and Zhang (2011), Weyl semimetals Yan and Felser (2017), and topological superconductors Qi and Zhang (2011), have been extensively discussed in the context of electronic systems. Following these remarkable discoveries, interacting spins and bosons have also been theoretically proposed to support symmetry-protected topological ground states with conventional bulk excitations but unusual gapless or gapped edge states Vishwanath and Senthil (2013); Levin and Gu (2012); Pollmann et al. (2012); Senthil and Levin (2013); Senthil (2015). Recently, there has been an exciting proposal that certain $S\!\!=\!\!1$ spin models on the diamond lattice may realize a time-reversal symmetry protected topological quantum paramagnet Wang et al. (2015), a stable 3D analogue of the $S\!\!=\!\!1$ Haldane chain Haldane (1983); Affleck et al. (1987).

This has led to a renewed interest in candidate spinel materials AB2O4 with A-site spins living on the diamond lattice. Previous studies of A-site magnetic spinels, such as MnSc2S4 ($S\!\!=\!\!5/2$) and CoAl2O4 ($S\!\!=\!\!3/2$), revealed degenerate spin spirals driven by frustration Fritsch et al. (2004); Tristan et al. (2005); Bergman et al. (2007); Bernier et al. (2008); Gao et al. (2016); Oitmaa (2019). On the other hand, FeSc2S4 shows weak Néel order in proximity to a non-magnetic ground state induced by spin-orbit coupling (SOC) Chen et al. (2009a, b); Plumb et al. (2016). The search for $S\!\!=\!\!1$ topological paramagnets has recently led to an intense investigation of NiRh2O4 using a variety of tools Chamorro et al. (2018).

NiRh2O4 is an unusual example of spin-$1$ 3d ions on the tetrahedrally coordinated A site, which is structurally stabilized by placing 4d Rh3+ ion at the octahedrally coordinated B-site. While NiRh2O4 is cubic at high temperature Blasse and Schipper (1963); Chamorro et al. (2018), it transforms into a tetragonal phase below $T\!\sim\!380$ K. Remarkably, in contrast to expectations from a Jahn-Teller mechanism which would favor $c/a\!<\!1$ and an $S\!\!=\!\!1$ ground state with quenched orbital angular momentum, the tetragonal phase is found to be elongated with $c/a\!\approx\!1.05$. Such a tetragonal distortion, with $c/a\!\!>\!\!1$, leaves the $t_{2}$ states of Ni partially filled, with orbital degrees of freedom unfrozen, allowing spin-orbit coupling (SOC) to play an important role. The mechanism for tetragonal distortion thus relies on SOC-induced orbital ordering, as previously discussed Tchernyshyov (2004); Maitra and Valenti (2007) in the context of the B-site active spinel ZnV2O4.

An early theoretical study Chen (2017) of NiRh2O4 considered a model with antiferromagnetic (AFM) first and second-neighbor Heisenberg exchanges ($J_{1}$ and $J_{2}$), applicable to frustrated spinels, and proposed that the non-magnetic ground state might arise from large single-ion anisotropy $DS_{z}^{2}$, with $D\!\!>\!\!0$ favoring local $S_{z}\!\!=\!\!0$. A pseudospin functional renormalization group study of the $J_{1}$-$J_{2}$ model Buessen et al. (2018) found that while the $S\!\!=\!\!1$ case favors a quantum spiral spin liquid ground state, the impact of tetragonal distortion or large $D/J_{1}\!\!\gtrsim\!\!8$ is to respectively favor Néel order or the $S_{z}\!\!=\!\!0$ ground state. Both studies effectively ignored orbital degrees of freedom. More recently, it was proposed Li and Chen (2018) that strong SOC with a tetrahedral crystal field could support a $J_{\rm eff}\!\!=\!\!0$ state at $d^{8}$ filling, generalizing the idea of $J_{\rm eff}\!\!=\!\!0$ insulators for $d^{4}$ filling in an octahedral crystal field Khaliullin (2013); Akbari and Khaliullin (2014); Svoboda et al. (2017); however, this might be overwhelmed by other energy scales (e.g., distortions or inter-site exchange) given weak SOC for Ni2+. On the experimental front, the inelastic neutron scattering (INS) results Chamorro et al. (2018) on NiRh2O4 were analyzed using spin-wave theory of an AFM state despite the absence of Néel order.

A satisfactory theoretical description of NiRh2O4 is thus lacking. Here, we combine first-principles density functional theory (DFT) and a model Hamiltonian study to unravel the curious case of NiRh2O4, explaining existing data and making predictions for future experiments.

*Density functional theory. — * We have carried out a first-principles study of NiRh2O4 in full-potential all electron approach of linear augmented plane wave (FLAPW) method Blaha et al. (2001), muffin-tin orbital method Andersen and Jepsen (1984); Andersen and Saha-Dasgupta (2000), as well as in pseudo-potential plane wave basis Kresse and Furthmüller (1996) with projected augmented potential (PAW) Blöchl (1994). The exchange-correlation functional was chosen to be generalized gradient approximation (GGA) Perdew et al. (1996), supplemented with onsite Hubbard correction GGA+$U$ Anisimov et al. (1993). Calculational details may be found in the Supplementary Material (SM) sm .

The electronic structure of NiRh2O4, calculated within GGA+$U$ ($U_{\rm Ni}$=5 eV, $J_{H}$=1 eV) resulted in half-metallic solutions for both the high temperature cubic and the low temperature tetragonal phases. Calculations show the spin splitting at Ni site to be large ($\approx\!\!1$ eV) while that at Rh site is an order of magnitude smaller ($\approx\!\!0.1$ eV), in accordance with the nominal magnetic and non-magnetic character of Ni2+ and Rh3+ respectively. In the high-symmetry cubic phase (see SM sm for details), the octahedral crystal field around Rh splits the 4$d$ states into $t_{2g}$ and $e_{g}$ with a large splitting $\sim\!\!3$ eV, while the tetrahedral crystal field around Ni splits the 3$d$ states into $e$ and $t_{2}$ with a relatively smaller splitting $\approx\!\!0.6$ eV. The $d$ states of high spin Ni are thus fully occupied in the up-spin channel; in the down-spin channel, the Ni $t_{2}$ states admixed with Rh $t_{2g}$ and O $p$ states cross the Fermi level (EF). The Rh $t_{2g}$ states are mostly occupied, except for the mixing with Ni states in down spin channel, while Rh $e_{g}$ states are empty. This is in accordance with nominal valence of Ni2+ with 2 holes in $t_{2}$ manifold, and low-spin nominally $d^{6}$ occupancy of Rh. This general picture remains valid also in the tetragonal phase as shown in Fig. 1. The tetragonal distortion, however introduces additional splitting among the cubic symmetry split states. This splits the Ni $t_{2}$ states with Ni $d_{xy}$ level positioned above Ni $d_{xz}/d_{yz}$ with splitting of $\approx$ 0.1 eV. One of the two holes of Ni thus occupies the down spin $d_{xy}$ level, while the other hole occupies the down spin doubly degenerate $d_{xz}/d_{yz}$ levels. This leaves the GGA+$U$ solution half-metallic even in the tetragonal phase, as shown in Fig. 1(a). The crystal and spin splittings at the tetragonal phase is shown in Fig. 1(c), which further highlights the energetic proximity of Ni $t_{2}$ and Rh $t_{2g}$ states in down-spin channel, driving the high degree of mixing between the two. This mixing gives rise to a small nonzero magnetic moment $\approx$ 0.06-0.07 $\mu_{B}$ at the otherwise nonmagnetic, low-spin, nominally $d^{6}$ Rh site, while the Ni moment is found to be 1.5 -1.6$\mu_{B}$. The remaining moment lives on O sites, giving rise to a net moment of 2 $\mu_{B}$/f.u in both cubic and tetragonal phases.

Given the active orbital degrees of freedom at Ni site, we next explore the effect of SOC. Within the GGA+$U$+SOC approach, the orbital state at Ni is derived from the $d_{xz}\pm id_{yz}$ orbitals. Due to partial occupancy of both orbitals, Ni develops a large orbital moment of $\sim\!\!1.0\mu_{B}$, supporting formation of a $S\!\!=\!\!1$, $L\!\!=\!\!1$ state. Repeating the calculation within GGA+SOC scheme, leads to a significantly smaller estimate of Ni orbital moment of $\approx$ 0.1 $\mu_{B}$, due to inability of GGA to capture the orbital polarization effect foo (a). While GGA+SOC splits the partially occupied orbitally degenerate states in down spin channel, this splitting is insufficient to open an insulating gap. This situation is similar to that discussed in case of FeCr2S4 Sarkar et al. (2009). The Coulomb correlation within GGA+$U$+SOC is thus crucial to produce a renormalized, large, orbital polarization ani which drives the system insulating, with a $\sim\!0.25$eV charge gap, as shown in Fig. 1(c).

We next estimate the Ni-Ni magnetic exchange from the knowledge of the effective hopping strengths and onsite energies in the Wannier basis of Ni-$t_{2}$ only low-energy Hamiltonian (see SM for details). The dominant AFM interactions in cubic phase turn out to be between four nearest-neighbor (NN) Ni sites ($J_{1}$), which belong to two different face-centered cubic (fcc) sublattices of the diamond lattice, and twelve next-nearest neighbor (NNN) Ni sites ($J_{2}$), which belong to the same fcc sublattice. The tetragonal distortion splits the twelve NNN Ni-Ni interactions into four in-plane ($J^{{}_{{}^{\prime}}}_{2}$) and eight out of plane ($J^{{}_{{}^{\prime\prime}}}_{2}$) interactions (see Fig. 2). The substantial mixing between Ni and Rh states, makes the Ni-O-Rh-O-Ni superexchange paths strong, as seen from the overlap of Wannier functions in Fig. 2 (see encircled part). The calculated exchanges are $J_{1}\!\!\approx\!\!1.2$ meV foo (b), with $J^{{}_{{}^{\prime}}}_{2},J^{{}_{{}^{\prime\prime}}}_{2}\!\!\approx\!0.4J_{1}$, showing strong magnetic frustration.

Single-site model. — Armed with the DFT results, we construct an effective single-site Hamiltonian for the $L\!\!=\!\!1$ and $S\!\!=\!\!1$ state, taking into account the tetragonal distortion ($\delta$) and SOC ($\lambda$);

[TABLE]

Based on DFT inputs, we consider the limit $\delta\gg\lambda$, and show that this leads to a simple, yet complete, understanding of the low temperature phenomenology of this distorted spinel.

In the regime $\delta\gg\lambda$, we start by constructing orbital eigenstates with well-defined $L_{z}$, which leads to a ground doublet with $L_{z}=\pm 1$ and an excited orbital singlet with $L_{z}=0$ which is split off by an energy $\delta$.

Next, let us take the spin degrees of freedom into account, which couple via SOC $\lambda\ll\delta$. The dominant SOC coupling is $\lambda L_{z}S_{z}$, which leads to a sequence of states in increasing order of energy which we label by $|L_{z},S_{z}\rangle$:

[TABLE]

with degeneracies shown in square brackets. We can perturbatively treat $\lambda(L^{+}S^{-}+L^{-}S^{+})/2$, since it only couples the low lying states at $E^{0}_{0,1,2}$ to the high energy states at $E^{0}_{3}$. Let us define the symmetric state $|e\rangle=(|+,-\rangle+|-,+\rangle)/\sqrt{2}$. We then find the sequence of states, with energies defined relative to the ground state,

[TABLE]

With these states and energies in hand, and a choice $\lambda\!\!\sim\!\!10$meV and $\delta\!\!\sim\!\!100$meV, we readily obtain a broad-brush understanding of some key experimental observations as summarized below. (The choice of $\delta\!\!\sim\!\!100$meV agrees with the spin-averaged crystal field splitting between $d_{xy}$ and $d_{yz}/d_{xz}$ orbitals from our DFT). We present further arguments against alternative scenarios in the SM sm .

Ground state: We find that the ground state is a non-magnetic singlet. This is consistent with the lack of any magnetic order down to the lowest temperature in this material foo (c). In contrast to previous proposals of non-magnetic $S_{z}=0$ state, our proposed state is a spin-orbit entangled “Schrodinger-cat” type state arising from weak off-diagonal SOC induced splitting of a doublet.

Thermodynamics: Since the gap to the states $|\psi_{4\pm}\rangle,|\psi_{5}\rangle$ are large, we expect to recover only an entropy $S_{\rm low}=R\ln 6$ for $T<300$K, consistent with specific heat measurements Chamorro et al. (2018) carried out up to room temperature (which corresponds to $T\ll\Delta_{4}$). At low temperatures, the state at $\Delta_{1}$ leads to a Schottky peak in $C/T$ at $T\sim 10$K from the level $|\psi_{1}\rangle$ (see SM sm ). It is not clear why this peak has not been observed; one possibility is that it may be affected by defects, which also likely lead to the observed spin freezing for $T\lesssim 6$K. The higher levels $|\psi_{2\pm}\rangle$ lead to a broad Schottky anomaly for $T\sim 30$-$40$K, similar to the experiments.

Neutron scattering: Our results for the local dynamical spin correlation function $S_{\rm loc}(\omega)$ are summarized in Fig. 3(a). The first excited state is nondegenerate, separated by an energy $\Delta_{1}\approx 2\lambda^{2}/\Delta\!\approx\!2$meV. We note that $|\psi_{0}\rangle$ and $|\psi_{1}\rangle$ are connected via $S^{z}$, so $|\psi_{1}\rangle$ should be visible in non-spin-flip scattering, but appears difficult to observe due to the resolution and the background, as well as possibly defects. The second excited state is a doublet $|\psi_{2,\pm}\rangle$ with an energy gap $\Delta_{2}\!\approx\!\lambda\!+\!\lambda^{2}/\Delta$. We propose that it is this doublet state which has been observed as a gapped mode in INS experiments Chamorro et al. (2018). The above parameter choice leads to the gap $\Delta_{2}\approx 11$meV, in crude agreement with the data. Based on our analysis, the states $|\psi_{3\pm}\rangle$ at an energy gap $\Delta_{3}\approx 22$meV and the singlet state $|\psi_{5}\rangle$ at a gap $\Delta_{4}\approx 108$meV are both “dark states”, invisible to neutrons due to vanishing matrix elements. Finally, $|\psi_{4\pm}\rangle$ with a gap $\Delta_{5}\approx 107$meV should be visible but with spectral weight much smaller that of $|\psi_{2\pm}\rangle$. This is a prediction for future INS experiments.

Magnetic susceptibility: The numerically computed single-site magnetic susceptibility $\chi$ can be fitted to an apparent “Curie-Weiss” form $\chi(T)=\chi_{0}+\alpha/(T-T_{0})$, with a negligible background $\chi_{0}\sim 10^{-5}$, an effective “Curie-Weiss” scale $T_{0}\approx 16(2)$K, and $\alpha\approx 0.85(2)$ (see SM) foo (d). In analyzing experiments, we expect $\chi_{0}$ will get lumped together with a background van Vleck type contribution which is conventionally subtracted. Our estimate for $T_{0}$ is small and “ferromagnetic” in sign, so that the $T_{0}^{\rm expt}\approx-11$K observed in experiments Chamorro et al. (2018) must be attributed to weak residual intersite AFM exchanges on the scale of $\sim\!\!1$ meV. Setting the fitted value of $\alpha\!\!\equiv\!\!S_{\rm eff}(S_{\rm eff}+1)/3$, yields an effective spin $S_{\rm eff}\!\!=\!\!1.4$ (or an effective magnetic moment $p_{\rm eff}\!\!\sim\!\!3.6\mu_{B}$), larger than a spin-only value $S\!\!=\!\!1$ as in experiments Chamorro et al. (2018).

Inter-site exchange. — We next incorporate inter-site interactions via a simple $J_{1}$-$J_{2}$ Heisenberg exchange model $H_{\rm ex}\!\!=\!\!\frac{1}{2}\sum_{i,j}J_{ij}{\bf S}_{i}\cdot{\bf S}_{j}$. In order to compute the spin dynamics in the low energy Hilbert space, we introduce, in the spirit of slave-boson theory Sachdev and Bhatt (1990); Li and Chen (2018), four local boson operators, $c^{\dagger}_{0},c^{\dagger}_{1},d^{\dagger}_{\pm}$, which respectively create states $|\psi_{0}\rangle$, $|\psi_{1}\rangle$, and $|\psi_{2\pm}\rangle$. Projecting the Heisenberg model to this Hilbert subspace, and imposing the local completeness constraint $c^{\dagger}_{0}c^{\phantom{\dagger}}_{0}\!+\!c^{\dagger}_{1}c^{\phantom{\dagger}}_{1}\!+\!d^{\dagger}_{\alpha}d^{\phantom{\dagger}}_{\alpha}\!=\!1$ (with an implicit sum on $\alpha=\pm$), we find that the site spin-$1$ operators may be approximated as $S^{z}=(c^{\dagger}_{1}c^{\phantom{\dagger}}_{0}+c^{\dagger}_{0}c^{\phantom{\dagger}}_{1})$ and $S^{\pm}=(c_{0}^{\dagger}\pm c^{\dagger}_{1})d^{\phantom{\dagger}}_{\pm}+d^{\dagger}_{\mp}(c^{\phantom{\dagger}}_{0}\mp c^{\phantom{\dagger}}_{1})$. At mean field level, we replace $c_{0}\!\to\!\langle c_{0}\rangle$, and retain leading powers in $\langle c_{0}\rangle$, to arrive at the Hamiltonian $H_{\rm tot}\!=\!H_{\rm loc}\!+\!H_{\rm ex}\!+\!H_{\rm con}$, where

[TABLE]

The different pieces correspond respectively to the local single-site Hamiltonian, the inter-site exchange Hamiltonian, and the constraint imposed (on average) via the Lagrange multiplier $\mu$. We note that the $c$ and $d$ bosons are decoupled at this order (except for the constraint). We can thus solve this in momentum space separately for these two sectors, leading to

[TABLE]

Here, $\sigma=\pm$, and the excitation energies are given by

[TABLE]

with $\gamma^{\phantom{\dagger}}_{{\bf{k}}}=\sum_{\ell_{1}}e^{i{{\bf{k}}}\cdot\ell_{1}}$ and $\eta^{\phantom{\dagger}}_{{\bf{k}}}=\sum_{\ell_{2}}e^{i{{\bf{k}}}\cdot\ell_{2}}$, where $\ell_{1},\ell_{2}$ are respectively the $4$ nearest-neighbor and $12$ next-neighbor vectors.

We choose $\Delta_{1}\!\!=\!\!1.8$meV and $\Delta_{2}\!\!=\!\!11$meV based on the single-site model, and $J_{1}\!\!=\!\!1.2$meV and $J_{2}/J_{1}\!\!=\!\!0.4$ from our DFT. Using these parameters, we minimize the ground state energy with respect to $\langle c_{0}\rangle^{2}$ while choosing $\mu$ to satisfy the constraint. We find the optimal $\langle c_{0}\rangle^{2}\!\approx\!0.7$ and $\mu\!\approx\!-2.1$meV.

The resulting weighted and powder-averaged dynamic spin structure factor relevant to INS experiments, ${\cal S}({{\bf{Q}}},\omega)=\sum_{\alpha}(1-Q^{2}_{\alpha}/Q^{2})S_{\alpha\alpha}({{\bf{Q}}},\omega)$, including a $1$meV broadening to mimic the experimental resolution but ignoring form factors, is plotted in Fig. 3(b) (see also SM sm ). The upper gapped mode, arising from the $|\psi_{2\pm}\rangle$ states, is in rough agreement with INS observations of a gapped dispersive mode Chamorro et al. (2018); we find that it really consists of two peaks due to two sublattices on the diamond lattice. The lower gapped mode is the “optical branch” of the $|\psi_{1}\rangle$ state. It collapses in energy, with increasing $Q$, from $\sim\!\!8$meV down to $\sim\!\!0.5$meV, and persists as an intense small-gap band, robust against magnetic condensate formation due to frustrating $J_{2}$ exchange. The lower energy “acoustic branch” of the $|\psi_{1}\rangle$ state is also gapped, but it has negligible intensity and is not visible here (see SM sm ). The small-$Q$ behavior depicted here may be partly masked by neutron kinematic constraints.

Summary and discussion. — We have combined DFT and model calculations to address the mystery of NiRh2O4, broadly capturing the existing thermodynamic and INS observations. In light of our work, it may be useful to revisit the low temperature specific heat and low energy INS on higher purity samples, and use INS to probe the predicted high energy crystal field level around $\sim 110$ meV. THz spectroscopy Laurita et al. (2015); Zhang et al. (2018) on NiRh2O4 could help to test our prediction of the “optical” $|\psi_{1}\rangle$ mode at $Q=0$, and infrared spectrocopy could be used to measure the insulating charge gap. It may be possible to use resonant inelastic X-ray scattering at a Ni-edge Lu et al. (2018) to look for the predicted $|\psi_{3\pm}\rangle$ and $|\psi_{5}\rangle$ “dark states” which are invisible to neutrons. Finally, our work suggests that NiRh2O4 does not realize a topological quantum paramagnet. However, it guides future searches by suggesting that tetragonal compression, presumably achievable by application of uniaxial strain, may provide the means to quench orbital angular momentum and suppress SOC effects, potentially stabilizing more exotic phases.

AP acknowledges support through a Discovery Grant from NSERC of Canada. TS-D thanks the Department of Science and Technology, India for the support through a Thematic Unit of Excellence.