# Supplementary material to Hidden covalent insulator and spin excitations in
SrRu2O6

Diana Csontosová  Department of Condensed Matter Physics, Faculty of Science,
Masaryk University, Kotlářská 2, 611 37 Brno, Czechia  Institute for Solid
State Physics, TU Wien, 1040 Vienna, Austria  Jiří Chaloupka  Department of
Condensed Matter Physics, Faculty of Science, Masaryk University, Kotlářská 2,
611 37 Brno, Czechia  Hiroshi Shinaoka  Department of Physics, Saitama
University, Saitama 338-8570, Japan  Atsushi Hariki  Department of Physics and
Electronics, Graduate School of Engineering, Osaka Metropolitan University,
1-1 Gakuen-cho, Nakaku, Sakai, Osaka 599-8531, Japan  Jan Kuneš  Institute for
Solid State Physics, TU Wien, 1040 Vienna, Austria

##  I Model and Computational Method

###  A Density functional calculation and wannierized Hamiltonian

Figure 1: Band structure of SrRu2O6 along $\Gamma-M-K-\Gamma-A$ path in the
first Brillouin zone obtained by DFT calculation. The red lines show $t_{2g}$
bands used in DMFT procedure.

Our DFT+DMFT calculations proceed in steps: (1) construction of an effective
model from a converged DFT calculation and (2) solution of the DMFT self-
consistent equation Georges  _et al._ (1996) followed by solution of the
Bethe-Salpeter equation with local two-particle irreducible vertex Zlatic and
Horvatic (1990); Jarrell (1992).

The DFT calculations were performed in the experimental $P{\overline{3}}1m$
structure Hiley  _et al._ (2014, 2015), see Fig. 1 of the main paper, using
the Wien2K package Blaha  _et al._ with the generalized gradient approximation
(GGA) Perdew  _et al._ (1996). The Wannier functions were constructed using
the WIEN2WANNIER Mostofi  _et al._ (2014) and WANNIER90 Kuneš  _et al._ (2010)
codes. For the $t_{2g}$-only case, we compared models with and without the
spin-orbit coupling (SOC). However, in order to get the orientation of
magnetic moments consistent with experiment (easy axis coincident with the
$c$-axis) we rotated the spin quantization axis to the direction of $c$-axis,
i. e. $\hkl(111)$ direction, in the presence of SOC. Below we show the local
Hamiltonians without ($H^{\mathrm{nsoc}}_{\mathrm{loc}}$) and with
($H^{\mathrm{soc}}_{\mathrm{loc}}$) SOC in the
$\\{d_{xy\uparrow},d_{yz\uparrow},d_{zx\uparrow},d_{xy\downarrow},d_{yz\downarrow},d_{zx\downarrow}\\}$
basis with respect to the local orbital coordinates

$\begin{pmatrix}\overrightarrow{x}_{\text{loc}},\overrightarrow{y}_{\text{loc}},\overrightarrow{z}_{\text{loc}}\end{pmatrix}=\pm\begin{pmatrix}\overrightarrow{x}_{\text{gl}},\overrightarrow{y}_{\text{gl}},\overrightarrow{z}_{\text{gl}}\end{pmatrix}\begin{pmatrix}0.0000000&0.1417008&-0.1417008\\\
-0.1626620&0.0818110&0.0818110\\\ 0.1156982&0.1156982&0.1156982\end{pmatrix}.$

The $\pm$ signs correspond to the two sublattices with Ru at $(1/3,2/3,1/2)$
and $(2/3,1/3,1/2)$, respectively. The global coordinates follow the standard
WIEN2k convention with $\overrightarrow{z}_{\text{gl}}$ parallel to the
hexagonal $c$-axis. The spin are set globally with quantization ($z$) axes
along the hexagonal $c$-axis. The local Hamiltonians on the two sites are
identical:

$\displaystyle
H^{\mathrm{nsoc}}_{\mathrm{loc}}=\begin{pmatrix}-61.3&-58.0&-58.2&0.0&0.0&0.0\\\
-58.0&-60.9&-58.0&0.0&0.0&0.0\\\ -58.2&-58.0&-61.4&0.0&0.0&0.0\\\
0.0&0.0&0.0&-61.3&-58.0&-58.2\\\ 0.0&0.0&0.0&-58.0&-60.9&-58.0\\\
0.0&0.0&0.0&-58.2&-58.0&-61.4\end{pmatrix}\times 10^{-3},$

$\displaystyle
H^{\mathrm{soc}}_{\mathrm{loc}}=\begin{pmatrix}-62.9&-61.2+27.8\mathrm{i}&-61.5-27.7\mathrm{i}&0.0&23.8+41.4\mathrm{i}&47.8\\\
-61.2-27.8\mathrm{i}&-62.7&-61.2+27.8\mathrm{i}&-23.8-41.4\mathrm{i}&0.0&23.8-41.4\mathrm{i}\\\
-61.5+27.7\mathrm{i}&-61.2-27.8\mathrm{i}&-62.9&-47.8&-23.8+41.4\mathrm{i}&0.0\\\
0.0&-23.8+41.4\mathrm{i}&-47.8&-62.9&-61.2-27.8\mathrm{i}&-61.5+27.7\mathrm{i}\\\
23.8-41.4\mathrm{i}&0.0&-23.8-41.4\mathrm{i}&-61.2+27.8\mathrm{i}&-62.7&-61.2-27.8\mathrm{i}\\\
47.8&23.8+41.4\mathrm{i}&0.0&-61.5-27.7\mathrm{i}&-61.2+27.8\mathrm{i}&-62.9\end{pmatrix}\times
10^{-3}.$

Using band structure obtained by DFT we further follow the standard DMFT
procedure, where the Hubbard model is mapped on the Anderson impurity model
(AIM). The AIM with the atomic interaction $H_{\rm atom}$ at Ru site

$\displaystyle H_{\rm atom}=$
$\displaystyle\sum_{\gamma\sigma,\gamma^{\prime}\sigma^{\prime}}h_{\gamma\sigma,\gamma^{\prime}\sigma^{\prime}}d_{\gamma\sigma}^{\dagger}d_{\gamma^{\prime}\sigma^{\prime}}+U\sum_{\gamma}n_{\gamma\uparrow}n_{\gamma\downarrow}+\sum_{\gamma>\gamma^{\prime},\sigma}[(U-2J_{\mathrm{H}})n_{\gamma\sigma}n_{\gamma^{\prime}-\sigma}$
$\displaystyle+$
$\displaystyle(U-3J_{\mathrm{H}})n_{\gamma\sigma}n_{\gamma^{\prime}\sigma}]+J_{\mathrm{H}}\sum_{\gamma\neq\gamma^{\prime},\sigma}(d_{\gamma-\sigma}^{\dagger}d_{\gamma^{\prime}\sigma}^{\dagger}d_{\gamma\sigma}^{\phantom{\dagger}}d_{\gamma^{\prime}-\sigma}^{\phantom{\dagger}}-d_{\gamma-\sigma}^{\dagger}d_{\gamma\sigma}^{\dagger}d_{\gamma^{\prime}\sigma}^{\phantom{\dagger}}d_{\gamma^{\prime}-\sigma}^{\phantom{\dagger}}),$

is solved using the continuous-time Quantum Monte Carlo method (CT-QMC) in the
hybridization expansion formalism Werner  _et al._ (2006); Gull  _et al._
(2011). $h_{\gamma\sigma,\gamma^{\prime}\sigma^{\prime}}$ describes the one-
particle interaction including the SOC if present. In practice, local basis
diagonalizing $h_{\gamma\sigma,\gamma^{\prime}\sigma^{\prime}}$ is employed.
Here, $\gamma$ and $\sigma$ represent the $t_{2g}$ orbitals and spin,
respectively, with $d_{\gamma\sigma}^{\dagger}$
($d_{\gamma\sigma}^{\phantom{\dagger}}$) being the corresponding creation
(annihilation) operators and $n_{\gamma\sigma}$ the number operators.

After self-consistency is achieved, physical quantities of the system are
obtained from the AIM with the renormalized bath and the lattice with the
converged self-energy $\Sigma(i\omega_{n})$. Besides the static one-particle
observables such as occupation numbers or ordered spin moments, the one-
particle spectral densities are obtained by analytic continuation of the self-
energy with maximum entropy method Wang  _et al._ (2009); Jarrell and
Gubernatis (1996); Kaufmann and Held (2021).

###  B Bethe-Salpeter equation

Once we have one-particle correlation function, we run two-particle CT-QMC
calculation in order to get local two-particle correlation function
$\chi_{ijs\nu,kls^{\prime},\nu^{\prime}}(\omega)$, where $i,j,k,l$ are the
spin-orbital indices, $s,s^{\prime}$ are the sublattice indices, and
$\nu,\nu^{\prime}$ are fermionic indices, which represent the imaginary time
evolution, e.g. Matsubara frequencies or imaginary time. The generalize
susceptibility $\chi_{ijs,kls^{\prime}}(\omega)$ is then obtained by summing
over fermionic frequencies

$\displaystyle\chi_{ijs,kls^{\prime}}(\omega)=-\frac{1}{\beta}\sum_{\nu,\nu^{\prime}}\chi_{ijs\nu,kls^{\prime},\nu^{\prime}}(\omega),$
(1)

where $\beta$ is inverse temperature. Since the final quantity which we want
to obtain is lattice dynamic susceptibility, we solve the Bethe-Salpeter
equation (BSE) first for impurity (2) from which we get local two-particle
irreducible vertex $\Gamma^{s}_{mn\nu_{1},pq\nu_{2}}(\omega)$ Kuneš (2011).
That is further inserted to the lattice BSE (3)

$\displaystyle\chi_{ij\nu,kl\nu^{\prime}}(s;\omega)=\chi^{0}_{ij\nu,kl\nu^{\prime}}(s;\omega)+\chi^{0}_{ij\nu,\underline{mn\nu_{1}}}(s;\omega)\Gamma^{s}_{\underline{mn\nu_{1}},\underline{pq\nu_{2}}}(\omega)\chi_{\underline{pq\nu_{2}},kl\nu^{\prime}}(s;\omega),$
(2) $\displaystyle
X_{ijs\nu,kls^{\prime}\nu^{\prime}}(\mathbf{q},\omega)=X^{0}_{ijs\nu,kls^{\prime}\nu^{\prime}}(\mathbf{q},\omega)+X^{0}_{ijs\nu,\underline{mns_{1}\nu_{1}}}(\mathbf{q},\omega)\Gamma^{\underline{s_{1}}}_{\underline{mn\nu_{1}},\underline{pq\nu_{2}}}(\omega)X_{\underline{pqs_{1}\nu_{2}},kls^{\prime}\nu^{\prime}}(\mathbf{q},\omega).$
(3)

The underlined indices are summed over. The lattice
$X^{0}_{ijs\nu,kls^{\prime}\nu^{\prime}}(\mathbf{q},\omega)$ and local bubbles
$\chi^{0}_{ij\nu,kl\nu^{\prime}}(s;\omega)$

$\displaystyle
X^{0}_{ijs\nu,kls^{\prime}\nu^{\prime}}(\mathbf{q},\omega)=-\delta_{\nu\nu^{\prime}}G_{is,ks^{\prime}}(\underline{\mathbf{k}}\\!+\\!\mathbf{q},\nu\\!+\\!\omega)G_{ls^{\prime},js}(\underline{\mathbf{k}},\nu)$
$\displaystyle\chi^{0}_{ij\nu,kl\nu^{\prime}}(s;\omega)=-\delta_{\nu\nu^{\prime}}G_{is,ks}(\underline{\mathbf{k}},\nu\\!+\\!\omega)G_{ls,js}(\underline{\mathbf{k}^{\prime}},\nu).$

are obtained from the 1P propagator

$G_{is,js^{\prime}}(\mathbf{k},\nu)=\left[i\nu-
h_{\mathbf{k}}-\Sigma(\nu)\right]^{-1}_{is,js^{\prime}}.$

As we can see, the fermionic Matsubara frequencies are dummy variables in BSE.
This means that the equation is invariant under their unitary transformation,
which allows us to use the basis of Legendre polynomials. In contrast to the
basis of fermionic Matsubara frequencies, 1P and 2P Green’s functions converge
rather fast in Legendre basis Boehnke  _et al._ (2011). That leads to a
sizable reduction of the computational and storage cost. In the present
calcuations we use mixed basis representation, where Legendre basis is used
for fermionic indices with $N=40$ Legendre polynoms, and Matsubara basis for
bosonic index for the first $N=10$ bosonic frequencies.

In order to compare our results with RIXS measurements, we calculate the
matrix elements of dynamic susceptibilities representing $S=3/2\rightarrow
1/2$ transitions $\langle X_{-\mathbf{q}}X_{\mathbf{q}}\rangle_{\omega}$ and
$\langle S^{x}_{-\mathbf{q}}S^{x}_{\mathbf{q}}\rangle_{\omega}$ corresponding
to transverse spin susceptibility

$\expectationvalue{\tilde{O}_{-\mathbf{q}}\tilde{O}_{\mathbf{q}}}_{\omega}=\int_{0}^{\beta}\\!\\!\\!\\!\\!\mathrm{d}\tau
e^{i\omega_{n}\tau}\expectationvalue{\tilde{O}_{-\mathbf{q}}(\tau)\tilde{O}_{\mathbf{q}}(0)}.$
(4)

Here $\expectationvalue{\tilde{O}}=\frac{1}{Z}\Tr\tilde{O}e^{-\beta H}$
denotes the thermal average and $\tilde{O}=O-\expectationvalue{O}$. Since the
calculations are performed in a two-site unit cell, we calculate the final
correlation function as a linear combination of contributions on the
sublattices $A$ and $B$ Niyazi  _et al._ (2021)

$\displaystyle\langle\tilde{O}_{-\mathbf{q}}\tilde{O}_{\mathbf{q}}\rangle_{\omega}$
$\displaystyle=\langle\tilde{O}_{-\mathbf{q},A}\tilde{O}_{\mathbf{q},A}\rangle_{\omega}+\langle\tilde{O}_{-\mathbf{q},B}\tilde{O}_{\mathbf{q},B}\rangle_{\omega}$
(5)
$\displaystyle+e^{-\mathrm{i}\mathbf{q}\mathbf{s}}\langle\tilde{O}_{-\mathbf{q},A}\tilde{O}_{\mathbf{q},B}\rangle_{\omega}+e^{\mathrm{i}\mathbf{q}\mathbf{s}}\langle\tilde{O}_{-\mathbf{q},B}\tilde{O}_{\mathbf{q},A}\rangle_{\omega},$
(6)

where $\mathbf{s}$ is sublattice vector given by the difference between
positions of ruthenium atoms with value $\left(1/3,1/3,0\right)$. Each term is
a contraction of the generalized lattice susceptibility
$\chi_{ijs,kls^{\prime}}(\mathbf{q},\omega)$ over the spin-orbital indices

$\begin{split}\frac{1}{N}\expectationvalue{\tilde{O}_{-\mathbf{q},s}\tilde{O}_{\mathbf{q},s^{\prime}}}_{\omega}=M_{\underline{ij}}M_{\underline{kl}}X_{\underline{ij}s,\underline{kl}s^{\prime}}(\mathbf{q},\omega),\end{split}$

with matrix elements $M_{ij}$ following from the definition of the operators
introduced in Eqs. (2) and (3) the main paper. Finally, performing analytical
continuation Levy  _et al._ (2017) from the imaginary frequency axis to the
real frequency axis and taking imaginary part of such quantities, we obtain
the dynamic susceptibility spectra.

##  II Results of DMFT calculations

###  A Magnetic moment and easy axis of magnetization

The presence of SOC induce single ion anisotropy (SIA) which further leads to
existence of easy magnetization axis and spin gap appearing in the magnon
spectra. The direction of easy axis is $\hkl(111)$, i. e. along $c$-axis. In
order to check if we got the system with easy magnetization axis in a correct
direction, we rotated magnetic moment around $y$ axis and let the calculation
converge. In Fig. 2 (a) we show all three components of the vector of magnetic
moment for both sites in a unit cell (dashed and full lines). The dotted lines
show the magnetic moment obtained by the calculation with no SOC. It is clear
that the direction of magnetic moment converges to $z$ axis, representing
$\hkl(111)$ direction.

Further, we plot the temperature dependence of magnetic moment for various
values of $J_{\mathrm{H}}$ (see Fig. 2 (b)). Values of magnetic moment are
multiplied by scaling factor 0.52 determined from the comparison of $m$
computed on Wannier functions and within a muffin tin sphere Hariki  _et al._
(2017). Below $\approx 600\,\mathrm{K}$ the magnetization for all
$J_{\mathrm{H}}$ is rather saturated. Consistently with previous study Hariki
_et al._ (2017), the magnetic moment for a low temperature does not differ
much with the change of $J_{\mathrm{H}}$.

Figure 2: (a) The convergence of magnetic moment $m$ in the system with SOC
included to the easy axis in $z$-direction ($c$-axis). (b) Temperature
dependence of magnetic moment for various values of Hund’s coupling
$J_{\mathrm{H}}$. Black dashed line is a polynomial fit of $m(T)$ for
$J_{\mathrm{H}}=0.19\,\mathrm{eV}$.

###  B Self-energies in paramagnetic phase

We show diagonal element of paramagnetic self-energies on the imaginary (see
Fig. 3) and the real frequency axis (see Fig. 4). The self-energies are plot
for various temperatures $T$ and Hund’s coupling $J_{\mathrm{H}}$, describing
the system in Mott and covalent insulating phase. For $J=0.19\,\mathrm{eV}$,
where the system is Mott insulator, we observe diverging imaginary part of
self-energy on the imaginary frequency axis with $\omega_{m}\rightarrow 0$,
which gives us a peak in $\omega=0$ on the real frequency axis. This leads to
the charge gap in the one-particle spectral function. For high temperatures,
above $T=464\,\mathrm{K}$ and $J_{\mathrm{H}}=0.16\,\mathrm{eV}$ the self-
energy possesses Mott-like behaviour. On the other hand, below
$T=580\,\mathrm{K}$ we observe Fermi liquid-like self-energy, where the
imaginary part on the imaginary axis goes to zero as $\omega_{m}$ goes to 0.
This gives us quadratic dependence for small $\omega$.

Figure 3: Diagonal element of self-energies
$\mathrm{Im}\\{\Sigma_{ii}(\omega_{m})\\}$ on imaginary axis for
$J=0.16\,\mathrm{eV}$ and $J=0.19\,\mathrm{eV}$, and various temperatures $T$.
Figure 4: Diagonal element of self-energies
$\mathrm{Im}\\{\Sigma_{ii}(\omega_{m})\\}$ on the real axis for
$J=0.16\,\mathrm{eV}$ and $J=0.19\,\mathrm{eV}$, and various temperatures $T$.
The crossover temperature $T^{\star}$ for $J$=0.16 eV lies between 464 K and
580 K.

###  C Transition from the states with $S=1/2$ to states with $S=3/2$

The matrix element of dynamic susceptibility generated by operator
$X_{\mathbf{R}s}=\left(\operatorname{d}^{\dagger}_{\mathbf{R}s1\uparrow}\operatorname{d}^{\phantom{\dagger}}_{\mathbf{R}s1\downarrow}-\operatorname{d}^{\dagger}_{\mathbf{R}s2\uparrow}\operatorname{d}^{\phantom{\dagger}}_{\mathbf{R}s2\downarrow}\right)+\mathrm{H.c.}$,
which is responsible for $S=1/2\rightarrow 3/2$ transition, is shown for
$J_{\mathrm{H}}=0.16\,\mathrm{eV}$ and $T=464\,\mathrm{K}$ with various
combination of inclusion of SOC and pair-hopping term in the interaction
Hamiltonian in Fig. 5. The colour scale is set to be the same for all three
figures. In the leftmost figure we show the system, where pair-hopping is
included and SOC is not. There is no mixing with the magnon spectra, i. e.
$\Delta S^{z}=\pm 1$ transitions while preserving the total spin $S=3/2$. When
the SOC is included, the mixing with magnon spectra is present (middle and the
rightmost figure). In the case of inclusion of both, SOC and pair-hopping, we
observe even stronger admixture of magnon. In this case, the data on the
imaginary frequency axis are noisier and become more sensitive to the
parameters of the analytical continuation, e. g. default model.

Figure 5: Matrix element of dynamic susceptibility showing the
$S=1/2\rightarrow 3/2$ transitions. Various systems are compared: without SOC
and pair-hopping (left), with SOC and without pair-hopping (middle), and with
both, SOC and pair-hopping, included (right).

###  D The effect of pair-hopping on magnon spectra

As we can see in Fig. 6, the pair-hopping term in the interaction Hamiltonian
has no qualitative effect on the magnon spectra. Since the ground state
quadruplet predominantly consists of the states with one electron per orbital,
the pair-hopping does not affect the transition spectra much. However, we can
observe a decrease of SIA - with pair-hopping $\Delta_{\mathrm{SIA}}=6.6$ meV,
while without pair-hopping $\Delta_{\mathrm{SIA}}=8.6$ meV. This leads to the
reduction of the size of the spin gap, what is well visible at the spectral
weight that is in $\Gamma^{\prime}$ larger in the system where pair-hopping is
included.

Figure 6: Comparison of magnon spectra (a) with and (b) without pair-hopping
included in the interaction Hamiltonian. The white line shows the spectral
weight $\Omega_{\mathbf{q}}$.

##  III Strong coupling expansion

To provide a comparison with a picture based on local $S=3/2$ moments, here we
consider the three-orbital Hubbard model used in the DMFT calculations and
perform a naive derivation of its effective low-energy model by standard
perturbation theory to second order in electronic hopping.

In the absence of the spin-orbit coupling
$\zeta\sum_{i}\mathbf{s}_{i}\cdot\mathbf{l}_{i}$, the ground-state manifold of
$d^{3}$ Ru5+ ions would correspond to pure spin-$3/2$ stabilized by Hund’s
coupling. At finite SOC strength $\zeta$, the four states forming spin-$3/2$
acquire an admixture of higher-energy states and split into two doublets with
$S^{z}\\!+\\!L^{z}=\pm\frac{3}{2},\pm\frac{1}{2}$. We map the resulting
eigenstates to pseudospin $\widetilde{S}=3/2$ and use them as a local basis of
our effective model. Pseudospins $\widetilde{S}=3/2$ are subject to single-ion
anisotropy term $-K(\widetilde{S}^{z})^{2}$ with the ionic contribution to $K$
being related to the above splitting via
$K_{0}=\frac{1}{2}\Delta_{\mathrm{SIA}}=\frac{1}{2}(E_{\pm\frac{1}{2}}-E_{\pm\frac{3}{2}})$.
As illustrated in Fig. 7(a), in our parameter regime $K_{0}$ appears roughly
quadratic to cubic in $\zeta$ which is consistent with the expectations based
on Ref. Liu  _et al._ , 2022. To simplify the calculations, we omit fine
details in the single-electron part of the ionic Hamiltonian obtained by DFT
and replace it by a two-parameter approximation assuming spherical symmetry of
the orbitals. To reproduce the single-electron spectrum of
$H^{\mathrm{soc}}_{\mathrm{loc}}$ given in Sec. I, we set the spin-orbit
coupling strength to $\zeta=0.1\>\mathrm{eV}$ and include the trigonal crystal
field splitting $\Delta=0.2\>\mathrm{eV}$. As in the DMFT calculations, we
take $U=2.7\>\mathrm{eV}$ and either $J_{\mathrm{H}}=0.16\>\mathrm{eV}$ or
$J_{\mathrm{H}}=0.19\>\mathrm{eV}$ to capture the local correlations.

In the honeycomb structure the Ru ions get connected via $90^{\circ}$
metal–oxygen–metal bonding path as well as by a significant direct overlap of
orbitals. We parametrize the hopping amplitudes in these two channels
according to

$\mathcal{H}_{t}(ij)=\left[t(d_{zx}^{\dagger}d_{yz}^{\phantom{\dagger}}+d_{yz}^{\dagger}d_{zx}^{\phantom{\dagger}})-t^{\prime}d_{xy}^{\dagger}d_{xy}^{\phantom{\dagger}}\right]_{ij}+\mathrm{H.c.}\,,$
(7)

where the orbitals are labeled using the main octahedral axes as a reference
frame. The above expression applies to nearest-neighbor Ru–Ru bonds where the
two metal–oxygen–metal links form a plaquette perpendicular to the octahedral
$z$ axis. Hopping rules for the other bond directions are obtained by a cyclic
permutation using the $C_{3}$ symmetry of the honeycomb lattice. Sum over spin
components is implied. The values of $t=0.33\>\mathrm{eV}$ and
$t^{\prime}=0.30\>\mathrm{eV}$ extracted from the DFT calculation are used to
model SrRu2O6.

The exchange interactions among pseudospins $\widetilde{S}=3/2$ are obtained
by eliminating the hopping part of the Hamiltonian via second order
perturbation theory. First we consider a simplified case with neglected
trigonal crystal field splitting and SOC, where the pseudospins become real
spins $S=3/2$. A straightforward calculation leads to a Heisenberg interaction
with the exchange parameter given by

$J=\frac{8}{9}\,\frac{t^{2}+\tfrac{1}{2}t^{\prime 2}}{U+2J_{\mathrm{H}}}\,.$
(8)

A detailed discussion of exchange interactions in a similar setting can be
found in Ref. Stavropoulos  _et al._ , 2021 focused on CrI3. Evaluated for our
parameters, Eq. (8) gives $J\approx 45\>\mathrm{meV}$ which is close to the
value $JS=65\>\mathrm{meV}$ obtained by fitting the magnon dispersion observed
by RIXS on SrRu2O6 Suzuki  _et al._ (2019).

As the next step, we have performed the same perturbative calculation
including the trigonal crystal field and SOC both in the pseudospin
$\widetilde{S}=3/2$ states as well as in the virtual states. On the ionic
level, this leads to only minor changes in the ionic eigenstates compared to
the pure spin-$3/2$ situation. Also the resulting exchange interaction is
still of almost purely Heisenberg form with the above $J$ value remaining
essentially unchanged (the difference is below $1\>\mathrm{meV}$). The main
effect of SOC combined with the trigonal crystal field thus lies in generating
the single-ion anisotropy $-K\widetilde{S}_{z}^{2}$ that comes both from the
ionic splitting $\Delta_{\mathrm{SIA}}$ of the $\widetilde{S}=3/2$ states and
from the exchange processes. In Fig. 7 we compare the two respective
contributions to $K=K_{0}+\Delta K$ that turn out to be of a similar
importance. The magnitude of the total $K$ is consistent with the fits of RIXS
data though the actual value is a bit larger than $KS=1.8\>\mathrm{meV}$ given
in Ref. Suzuki  _et al._ , 2019.

Figure 7:  (a) Contributions to the single-ion anisotropy term
$-K(\widetilde{S}^{z})^{2}$ for varying SOC strength $\zeta$. The ionic
contribution $K_{0}$ is determined from the splitting of the low-energy
quartet via
$\Delta_{\mathrm{SIA}}=2K_{0}=E_{\pm\frac{1}{2}}-E_{\pm\frac{3}{2}}$. $\Delta
K$ is the exchange correction obtained by second order perturbation theory in
electronic hopping. Apart from the parameter values indicated in the figure,
we have used $U=2.7\>\mathrm{eV}$ and the trigonal crystal field
$\Delta=0.2\>\mathrm{eV}$. (b) Exchange correction depending on the hopping
amplitudes. The ratio $t/t^{\prime}$ is kept fixed. Dashed lines indicate the
parameter values used for SrRu2O6.

##  IV Exact diagonalization of the three-orbital Hubbard model on a bond

It is instructive to inspect the potential issues of the pseudospin-$3/2$
picture utilizing a tiny system where an exact solution of the underlying
Hubbard model is available. To this end we diagonalize a two-site version of
our Hubbard model acting on an isolated bond and check whether the low-energy
part of its spectrum can be satisfactorily captured by the effective
Heisenberg+SIA model with the Hamiltonian
$J\boldsymbol{\widetilde{S}}_{1}\cdot\boldsymbol{\widetilde{S}}_{2}-K[(\widetilde{S}_{1}^{z})^{2}+(\widetilde{S}_{2}^{z})^{2}]$.

Figure 8:  (a) Energy levels of our Hubbard model on an isolated bond
depending on the hopping amplitudes. The calculation was performed for
$U=2.7\>\mathrm{eV}$, $J_{\mathrm{H}}=0.19\>\mathrm{eV}$,
$\zeta=0.10\>\mathrm{eV}$, $\Delta=0.20\>\mathrm{eV}$. Both $t$ and
$t^{\prime}$ were employed with the ratio $t/t^{\prime}$ kept fixed at
$t/t^{\prime}=0.33/0.30$. The energies are plotted relative to the ground-
state level. First three excited levels are shown in blue, next five by green,
and next seven by red. This coloring highlights the quasidegeneracy related to
the level structure of the spin-$3/2$ Heisenberg model on the bond and is used
in all the panels (a)–(e). (b) Closer view on the low-energy part of panel
(a). The dashed line shows the Heisenberg exchange value $J$ obtained by
second order perturbation theory. In small $t$ limit, it coincides with the
first excitation energy. (c) Squared total bond pseudospin
$\boldsymbol{\widetilde{S}}_{\mathrm{tot}}=\boldsymbol{\widetilde{S}}_{1}+\boldsymbol{\widetilde{S}}_{2}$
as measured in the selected low-energy eigenstates after projecting them onto
$d^{3}-d^{3}$ configuration subspace. Crosses indicate values obtained after a
further projection to the $\widetilde{S}=3/2$ basis. (d) The same for the
pseudospin at one of the sites. The average
$\langle\widetilde{S}^{2}\rangle=\langle\widetilde{S}^{2}_{1}\rangle=\langle\widetilde{S}^{2}_{2}\rangle$
takes the value $\frac{15}{4}=3.75$ in the ideal pseudospin-$3/2$ situation.
(e) Conditional probability of having both sites populated by pseudospin
$\widetilde{S}=3/2$ provided that there is a $d^{3}-d^{3}$ configuration on
the bond. This probability is again evaluated for the selected low-energy
eigenstates. (f) Probability of encountering various combinations of $d^{3}$
ionic eigenstates in the $d^{3}-d^{3}$ part of the bond ground state. There
are 20 such $d^{3}$ eigenstates that can be grouped to (i) the low-energy
quasidegenerate pseudospin-$3/2$ quartet, (ii) ten states at the relative
energy $\approx 0.5\>\mathrm{eV}$, and (iii) six states around
$1\>\mathrm{eV}$.

Such a comparison is performed in Fig. 8 for the
$J_{\mathrm{H}}=0.19\>\mathrm{eV}$ parameter set. In the $J\gg K$ case that
corresponds to our parameter regime, we should find low-energy eigenstates
having total bond pseudospin $\widetilde{S}_{\mathrm{tot}}=0$, $1$, $2$, and
$3$ at the relative energies $(E_{\widetilde{S}_{\mathrm{tot}}}-E_{0})/J=0$,
$1$, $3$, and $6$, respectively. The low-energy spectrum shown in Fig.
8(a),(b) indeed reveals this kind of structure, comprising $1$-, $3$-, $5$-,
and $7$-fold quasidegenerate levels. A comparison to the Heisenberg exchange
$J$ obtained for the same parameters by the second-order perturbation method
of Sec. III indicates that the perturbative approach is well applicable up to
about $t\approx 0.2\>\mathrm{eV}$. Above this point, a quantitative
discrepancy of pertubative $J$ and the first excitation energy sets in and
rapidly grows. This problem concerns the lower states, the energy distance of
$\widetilde{S}_{\mathrm{tot}}=2$ and $3$ levels is consistent with the
perturbative $J$ in the whole range shown.

One possible interpretation of the distorted level structure at larger $t$ is
that higher-order processes introduce new interaction terms beyond the
Heisenberg one e.g. a biquadratic interaction. However, even though one may
formally reproduce the observed energy level shifts by extending the
pseudospin model, the source of the discrepancy above $t\approx
0.2\>\mathrm{eV}$ lies actually in the inadequacy of the pseudospin
$\widetilde{S}=3/2$ picture itself, as it is illustrated by the remaining
panels of Fig. 8. A first indication comes from the total bond pseudospin
$\boldsymbol{\widetilde{S}}_{\mathrm{tot}}=\boldsymbol{\widetilde{S}}_{1}+\boldsymbol{\widetilde{S}}_{2}$
in the form of the average $\langle\widetilde{S}_{\mathrm{tot}}^{2}\rangle$
which is expected to take the values
$\widetilde{S}_{\mathrm{tot}}(\widetilde{S}_{\mathrm{tot}}+1)=0,2,6,12$ for
any kind of isotropic pseudospin $\widetilde{S}=3/2$ model. In Fig. 8(c) we
plot this quantity evaluated for the 16 lowest states resembling the
Heisenberg model spectrum. We took only the $d^{3}-d^{3}$ part of the
eigenstates to cut off $d^{2}-d^{4}$ and other charge excitations beyond the
effective model. Similarly to the previous observation, above $t\approx
0.2\>\mathrm{eV}$ the average $\langle\widetilde{S}_{\mathrm{tot}}^{2}\rangle$
starts to to deviate strongly from its expected value, namely for the ground
state and to a less extent for the triplet excited states. The expected values
of $\langle\widetilde{S}_{\mathrm{tot}}^{2}\rangle$ are only recovered after a
projection to the $\widetilde{S}=3/2$ subspace. Note that the deviations at
small $t$ are due to the change of the parameter regime – the anisotropy $K$
takes over $J\propto t^{2}/U$ in that case and the behavior of
$\langle\widetilde{S}_{\mathrm{tot}}^{2}\rangle$ changes completely. For
example, the ground-state $\langle\widetilde{S}_{\mathrm{tot}}^{2}\rangle$
shows a hump at low $t$ which is consistent with
$\langle\widetilde{S}_{\mathrm{tot}}^{2}\rangle=3\,[1-1/\sqrt{1+(4K/3J)^{2}}]$
obtained for the ground state of the Heisenberg+SIA model on the bond. The
departure from the pseudospin $\widetilde{S}=3/2$ picture starting around
$t\approx 0.2\>\mathrm{eV}$ is also well visible in Fig. 8(d) showing the on-
site average
$\langle\widetilde{S}_{1}^{2}\rangle=\langle\widetilde{S}_{2}^{2}\rangle$ that
measures the pseudospin length. In Fig. 8(e), we take another point of view
and consider directly the conditional probability of having two pseudospins
$\widetilde{S}=3/2$ in the projected $d^{3}$–$d^{3}$ part of the eigenstates.
This probability starts to drop significantly around $t\approx
0.2\>\mathrm{eV}$ suggesting again the incompleteness of the pseudospin
picture.

Based on the previous analysis, our target values of $t=0.33\>\mathrm{eV}$ and
$t^{\prime}=0.30\>\mathrm{eV}$ belong already to the parameter regime
characterized by a large participation of the higher-lying multiplet states.
Motivated by the question whether some particular higher states are prominent
in their participation, we take the corresponding bond ground state, project
it to $d^{3}-d^{3}$ configuration subspace and evaluate the probabilities to
encounter various combinations of the multiplet states. The result presented
in Fig. 8(f) shows a rather homogeneous involvement of the multiplet states
around $0.5\>\mathrm{eV}$ as well those around $1\>\mathrm{eV}$, even though a
bit weaker in the latter case. Dominant in the probability diagram are
naturally the four possible combinations of pseudospin $\widetilde{S}=3/2$
that may participate in $\widetilde{S}_{\mathrm{tot}}=0$ state, i.e.
$|+\frac{3}{2}\rangle\otimes|-\frac{3}{2}\rangle$,
$|-\frac{3}{2}\rangle\otimes|+\frac{3}{2}\rangle$,
$|+\frac{1}{2}\rangle\otimes|-\frac{1}{2}\rangle$,
$|-\frac{1}{2}\rangle\otimes|+\frac{1}{2}\rangle$. Interestingly, the off-
diagonal blocks combining the pseudospin $\widetilde{S}=3/2$ states with
higher multiplet states are nearly vanishing, giving the $\widetilde{S}=3/2$
sector some kind of independence. It may be understood as a consequence of
“near” spin conservation that is not too disturbed by the presence of SOC.
Under spin conservation the singlet bond ground state cannot contain a spin
$S=3/2$ state on one site together with any of the higher states at the other
site since they all have total spin $S=1/2$. Even though the single-bond
results need to be considered with care, we may speculate that this type of
sector separation may help the $\widetilde{S}=3/2$ pseudospins to retain some
of their features, e.g. the Heisenberg-like dynamics, despite the activation
of the higher multiplet states. In such case, the corresponding “dilution” of
the pseudospin-$3/2$ system could be possibly accounted for simply by a
renormalization of its parameters.

## References

  * Georges  _et al._ (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
  * Zlatic and Horvatic (1990) V. Zlatic and B. Horvatic, Solid State Commun. 75, 263 (1990).
  * Jarrell (1992) M. Jarrell, Phys. Rev. Lett. 69, 168 (1992).
  * Hiley  _et al._ (2014) C. I. Hiley, M. R. Lees, J. M. Fisher, D. Thompsett, S. Agrestini, R. I. Smith, and R. I. Walton, Angew. Chem. Int. Ed. 53, 4423 (2014).
  * Hiley  _et al._ (2015) C. I. Hiley, D. O. Scanlon, A. A. Sokol, S. M. Woodley, A. M. Ganose, S. Sangiao, J. M. De Teresa, P. Manuel, D. D. Khalyavin, M. Walker, M. R. Lees, and R. I. Walton, Phys. Rev. B 92, 104413 (2015).
  * (6) P. Blaha, K. Schwarz, G. Madsen, D. Kvasnicka, and J. Luitz,  _WIEN2k, An Augmented Plane Wave + Local Orbitals Program for Calculating Crystal Properties (Karlheinz Schwarz, Techn. Universitat Wien, Austria, 2001), ISBN 3-9501031-1-2_.
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Mostofi  _et al._ (2014) A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. 185, 2309 (2014).
  * Kuneš  _et al._ (2010) J. Kuneš, R. Arita, P. Wissgott, A. Toschi, H. Ikeda, and K. Held, Comput. Phys. Commun. 181, 1888 (2010).
  * Werner  _et al._ (2006) P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006).
  * Gull  _et al._ (2011) E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Rev. Mod. Phys. 83, 349 (2011).
  * Wang  _et al._ (2009) X. Wang, E. Gull, L. de’ Medici, M. Capone, and A. J. Millis, Phys. Rev. B 80, 045101 (2009).
  * Jarrell and Gubernatis (1996) M. Jarrell and J. Gubernatis, Phys. Rep. 269, 133 (1996).
  * Kaufmann and Held (2021) J. Kaufmann and K. Held, (2021), arXiv:2105.11211 .
  * Kuneš (2011) J. Kuneš, Phys. Rev. B 83, 085102 (2011).
  * Boehnke  _et al._ (2011) L. Boehnke, H. Hafermann, M. Ferrero, F. Lechermann, and O. Parcollet, Phys. Rev. B 84, 075145 (2011).
  * Niyazi  _et al._ (2021) A. Niyazi, D. Geffroy, and J. Kuneš, Phys. Rev. B 104, 075152 (2021).
  * Levy  _et al._ (2017) R. Levy, J. LeBlanc, and E. Gull, Computer Physics Communications 215, 149 (2017).
  * Hariki  _et al._ (2017) A. Hariki, A. Hausoel, G. Sangiovanni, and J. Kuneš, Phys. Rev. B 96, 155135 (2017).
  * Liu  _et al._ (2022) X. Liu, D. Churchill, and H.-Y. Kee, Phys. Rev. B 106, 035122 (2022).
  * Stavropoulos  _et al._ (2021) P. P. Stavropoulos, X. Liu, and H.-Y. Kee, Phys. Rev. Research 3, 013216 (2021).
  * Suzuki  _et al._ (2019) H. Suzuki, H. Gretarsson, H. Ishikawa, K. Ueda, Z. Yang, H. Liu, H. Kim, D. Kukusta, A. Yaresko, M. Minola, J. A. Sears, S. Francoual, H.-C. Wille, J. Nuss, H. Takagi, B. J. Kim, G. Khaliullin, H. Yavaş, and B. Keimer, Nat. Mater 18, 563 (2019).