# Curvature Dependent Diffusion Flow on Surface with Thickness

Naohisa Ogawa<EMAIL_ADDRESS>Hokkaido Institute of Technology, Sapporo
006-8585 Japan

###### Abstract

Particle diffusion in a two dimensional curved surface embedded in $R_{3}$ is
considered. In addition to the usual diffusion flow, we find a new flow with
an explicit curvature dependence. New diffusion equation is obtained in
$\epsilon$ (thickness of surface) expansion. As an example, the surface of
elliptic cylinder is considered, and curvature dependent diffusion coefficient
is calculated.

###### pacs:

87.10.-e, 02.40.Hw, 02.40.Ma, 82.40.Ck

##  I Motivation

The particle motion on a given curved surface is old but interesting problem
in wide range of physics. Especially the diffusion process of particles on
such a manifold is still an open problem, and related to various kinds of
phenomena.

For example the motion of protein on cell membrane has great importance in
biophysics. There are several research papers discussing on this problem. Some
of them are treating this problem by using usual diffusion equation with
curved coordinate, and discuss the curvature (Gauss curvature) dependence of
its solution diffusion_equation . Other of them use the Langevin equation on
curved surface and calculating the curvature dependence of diffusion
coefficient langevin_equation .

The quantum mechanics of particle motion on such a curved manifold is also
considered by many authors. This problem is usually explained by the
Schroedinger equation with Laplace-Beltrami operator. However, when we treat
the curved surface as embedded one in 3 dimensional Euclidean space, situation
is changed and then we have a quantum potential term related to the curvature
additional to the kinetic operator da Costa ,ogawa_fujii ,fujii .

Another example is in larger scale physics in which our consideration is
devoted. Patterns of animal skins are well expressed by the reaction diffusion
equation Turing . But the patterns are different for each parts even in one
individual. For example, Char fish, the side part has white spot pattern, but
the back part has labyrinth pattern. (For these two patterns, see for example
Shoji_Iwasa .) One of the reasons might come from the curvature difference
between side part and back part. If the diffusion is influenced by the
curvature, this difference of patterns might be explained. Furthermore, the
cross section of fish has form of ellipsoid and the surface can be
approximated as the one of elliptic cylinder. In two dimensional space, we
have only two kinds of curvature, one is Gauss curvature and other is mean
curvature. Both are constructed from second fundamental tensor by taking
determinant or trace. Gauss curvature can also be constructed only by metric
tensor and its derivatives, but this is not the case for the mean curvature.
The elliptic cylinder, in which we have much interest, has zero Gauss
curvature and non-zero mean curvature. Therefore to explain the pattern change
of Char fish, solution of the diffusion equation should depend on mean
curvature. This is impossible if we start from usual diffusion equation
because it depends only on metric but not on second fundamental tensor.
Therefore we need some new diffusion equation, which bring not only Gauss
curvature but also mean curvature. In this article, we discuss how to
construct such curvature dependent diffusion equation.

##  II Coordinate and Metric

The simple extension of diffusion equation in Euclidean space to Riemannian
space can be done by changing Laplacian with Cartesian coordinate to the one
with curved coordinate, i.e. Laplace Beltrami operator. This coordinate change
is not enough for our purpose, however. The way of construction of new
diffusion equation in this paper is the followings. We re-identify the two
dimensional diffusion as the limiting process from three dimensional
diffusion. We place the curved surface $\Sigma$ in three dimensional Euclidean
space $R_{3}$, and we put two similar copies of $\Sigma$, called
$\tilde{\Sigma}$ and $\Sigma^{\prime}$ at a small distance of $\epsilon/2$.
Our particles can only move between these two surfaces, and later we take a
limit $\epsilon\to 0$. We look for the form of diffusion equation in this
limit. The coordinates we use hereafter is the followings. (See fig.1)

$\vec{X}$ is the Cartesian coordinate in $R_{3}$. $\vec{x}$ is the Cartesian
coordinate which specifies only the points on $\Sigma$. $q^{i}$ is the curved
coordinate on $\Sigma$. (Small Latin indices $i,j,k,\cdots$ runs from 1 to 2.)
$q^{0}$ is the coordinate in $R_{3}$ normal to $\Sigma$. Further by using the
normal unit vector $\vec{n}(q^{1},q^{2})$ on $\Sigma$ at point
$(q^{1},q^{2})$, we can identify any points between two surfaces
$\Sigma^{\prime}$ and $\tilde{\Sigma}$ by the following thin-layer
approximation fujii .

$\vec{X}(q^{0},q^{1},q^{2})=\vec{x}(q^{1},q^{2})+q^{0}\vec{n}(q^{1},q^{2}),$
(1)

where $-\epsilon/2\leq q^{0}\leq\epsilon/2$.

Figure 1: Embedding and Coordinate

From this relation we can obtain the curvilinear coordinate system between two
surfaces ($\subset R_{3}$) by the coordinate $q^{\mu}=(q^{0},q^{1},q^{2})$,
and metric $G_{\mu\nu}$. (Hereafter Greek indices $\mu,\nu,\cdots$ runs from 0
to 2.)

$G_{\mu\nu}=\frac{\partial\vec{X}}{\partial
q^{\mu}}\cdot\frac{\partial\vec{X}}{\partial q^{\nu}}.$ (2)

Each part of $G_{\mu\nu}$ is the following.

$G_{ij}=g_{ij}+q^{0}(\frac{\partial\vec{x}}{\partial
q^{i}}\cdot\frac{\partial\vec{n}}{\partial
q^{j}}+\frac{\partial\vec{x}}{\partial
q^{j}}\cdot\frac{\partial\vec{n}}{\partial
q^{i}})+(q^{0})^{2}\frac{\partial\vec{n}}{\partial
q^{i}}\cdot\frac{\partial\vec{n}}{\partial q^{j}},$ (3)

where

$g_{ij}=\frac{\partial\vec{x}}{\partial
q^{i}}\cdot\frac{\partial\vec{x}}{\partial q^{j}}$ (4)

is the metric on $\Sigma$. Hereafter indices $i,j,k\cdots$ are lowered or
rised by $g_{ij}$ and its inverse $g^{ij}$. We also obtain

$G_{0i}=G_{i0}=0,~{}~{}G_{00}=1.$ (5)

We can proceed the calculation by using the new variables. We first define the
tangential vector to $\Sigma$ by

$\vec{B}_{k}=\frac{\partial\vec{x}}{\partial q^{k}}.$ (6)

Note that $\vec{n}\cdot\vec{B}_{k}=0$. Then we obtain two relations.

Gauss equation:

$\frac{\partial\vec{B}_{i}}{\partial
q^{j}}=-\kappa_{ij}\vec{n}+\Gamma^{k}_{ij}\vec{B}_{k},$ (7)

Weingarten equation:

$\frac{\partial\vec{n}}{\partial q^{j}}=\kappa_{j}^{m}\vec{B}_{m},$ (8)

where

$\Gamma^{k}_{ij}\equiv\frac{1}{2}g^{km}(\partial_{i}g_{mj}+\partial_{j}g_{im}-\partial_{m}g_{ij}).$

$\kappa_{ij}$ is called Euler-Schauten tensor, or second fundamental tensor
defined as

$\kappa_{ij}=\frac{\partial\vec{n}}{\partial q^{i}}\cdot\vec{B}_{j}.$ (9)

The second fundamental tensor $\kappa_{ij}$ is the projection of
$\partial\vec{n}$ into the surface. Furthermore, the mean curvature is given
by

$\kappa=g^{ij}\kappa_{ij},$ (10)

and Ricci scalar curvature $R$ is obtained by

$R/2=\det(g^{ik}\kappa_{kj})=\det(\kappa^{i}_{j})=\frac{1}{2}(\kappa^{2}-\kappa_{ij}\kappa^{ij}).$
(11)

Then we have the formula for metric of curvilinear coordinate in a
neighborhood of $\Sigma$.

$G_{ij}=g_{ij}+2q^{0}\kappa_{ij}+(q^{0})^{2}\kappa_{im}\kappa^{m}_{j}.$ (12)

Under the inversion $q^{0}\to-q^{0}$, we have $\kappa_{ij}\to-\kappa_{ij}$ as
well as $\vec{n}\to-\vec{n}$ from
$\vec{n}=\partial_{0}\vec{X}/\mid\partial_{0}\vec{X}\mid$. Therefore $G_{ij}$
is invariant under $q^{0}\to-q^{0}$.

Now we have the total metric tensor such as,

$G_{\mu\nu}=\left(\begin{array}[]{cc}1&~{}~{}~{}0{}{}{}\\\
0&G_{ij}\end{array}\right).$ (13)

##  III Embedding of Diffusion field

Let us denote 3 dimensional diffusion field as $\phi^{(3)}$, and Laplacian as
$\Delta^{(3)}$. Then we have the equation with normalization condition

$\displaystyle\frac{\partial\phi^{(3)}}{\partial t}=D\Delta^{(3)}\phi^{(3)},$
(14) $\displaystyle 1$ $\displaystyle=$
$\displaystyle\int\phi^{(3)}(q^{0},q^{1},q^{2})\sqrt{G}~{}d^{3}q,$ (15)

where $D$ is the diffusion constant, and $G=\det(G_{\mu\nu})=\det(G_{ij})$.
Our aim is to construct the effective two dimensional diffusion equation from
3D equation above.

$\displaystyle\frac{\partial\phi^{(2)}}{\partial
t}=D\Delta^{(eff)}\phi^{(2)},$ (16) $\displaystyle 1$ $\displaystyle=$
$\displaystyle\int\phi^{(2)}(q^{1},q^{2})\sqrt{g}~{}d^{2}q,$ (17)

where $\phi^{(2)}$ is the two dimensional diffusion field, $g=\det(g_{ij})$,
and $\Delta^{(eff)}$ is unknown effective 2D diffusion operator which might
not be equal to simple 2D Laplace Beltrami operator.

From two normalization conditions, we obtain

$\displaystyle 1$ $\displaystyle=$
$\displaystyle\int\phi^{(3)}(q^{0},q^{1},q^{2})\sqrt{G}~{}d^{3}q,$
$\displaystyle=$
$\displaystyle\int[\int_{-\epsilon/2}^{\epsilon/2}dq^{0}(\phi^{(3)}\sqrt{G/g})]~{}\sqrt{g}~{}d^{2}q,$
$\displaystyle=$ $\displaystyle\int\phi^{(2)}(q^{1},q^{2})\sqrt{g}~{}d^{2}q.$

Therefore we obtain the relation,

$\phi^{(2)}(q^{1},q^{2})=\int_{-\epsilon/2}^{\epsilon/2}\tilde{\phi}^{(3)}dq^{0},$
(18)

where

$\tilde{\phi}^{(3)}\equiv\phi^{(3)}\sqrt{G/g}.$ (19)

We multiply $\sqrt{G/g}$ to equation (14) and integrate by $q^{0}$, then we
obtain

$\frac{\partial\phi^{(2)}}{\partial
t}=D\int_{-\epsilon/2}^{\epsilon/2}\tilde{\Delta}^{(3)}\tilde{\phi}^{(3)}dq^{0},$
(20)

where

$\tilde{\Delta}^{(3)}\equiv\sqrt{G/g}~{}~{}\Delta^{(3)}\sqrt{g/G}.$ (21)

Next we analyze new operator $\tilde{\Delta}^{(3)}$. From the form of Laplace
Beltrami operator

$\Delta^{(3)}=G^{-1/2}\frac{\partial}{\partial
q^{\mu}}G^{1/2}G^{\mu\nu}\frac{\partial}{\partial q^{\nu}},$

we have

$\displaystyle\tilde{\Delta}^{(3)}$ $\displaystyle=$ $\displaystyle
g^{-1/2}\frac{\partial}{\partial
q^{\mu}}G^{1/2}G^{\mu\nu}\frac{\partial}{\partial q^{\nu}}(g/G)^{1/2}$ (22)
$\displaystyle=$ $\displaystyle\tilde{\Delta}^{(2)}+\tilde{\Delta}^{(1)},$

where

$\tilde{\Delta}^{(2)}\equiv g^{-1/2}\frac{\partial}{\partial
q^{i}}G^{1/2}G^{ij}\frac{\partial}{\partial q^{j}}(g/G)^{1/2},$ (23)

and

$\tilde{\Delta}^{(1)}\equiv\frac{\partial}{\partial
q^{0}}G^{1/2}\frac{\partial}{\partial q^{0}}G^{-1/2}.$ (24)

Then our diffusion equation has form

$\frac{\partial\phi^{(2)}}{\partial
t}=D\int_{-\epsilon/2}^{\epsilon/2}\tilde{\Delta}^{(2)}\tilde{\phi}^{(3)}dq^{0}.$
(25)

The contribution from $\tilde{\Delta}^{(1)}$ vanishes because

$\displaystyle\int_{-\epsilon/2}^{\epsilon/2}\tilde{\Delta}^{(1)}\tilde{\phi}^{(3)}dq^{0}$
$\displaystyle=$ $\displaystyle
g^{-1/2}\int_{-\epsilon/2}^{\epsilon/2}\frac{\partial}{\partial
q^{0}}(G)^{1/2}\frac{\partial}{\partial q^{0}}\phi^{(3)}~{}dq^{0}$ (26)
$\displaystyle=$ $\displaystyle
g^{-1/2}[~{}(G)^{1/2}\frac{\partial\phi^{(3)}}{\partial
q^{0}}]\mid_{-\epsilon/2}^{\epsilon/2}=0.$

The last equality is the requirement that diffusion flow does not pass through
the surface: $\Sigma^{\prime}$ and $\tilde{\Sigma}$.

Now we calculate r.h.s of (25) up to ${\cal O}(\epsilon^{2})$. Since we have

$\tilde{\phi}^{(3)}={\cal O}(\epsilon^{-1}),$ (27)

from (18), we need to expand $\tilde{\Delta}^{(2)}$ up to ${\cal
O}(\epsilon^{2})$. The following relations are useful

$\displaystyle G_{ij}$ $\displaystyle=$ $\displaystyle
g_{ij}+2q^{0}\kappa_{ij}+(q^{0})^{2}\kappa_{im}\kappa^{m}_{j},$ (28)
$\displaystyle G^{ij}$ $\displaystyle=$ $\displaystyle
g^{ij}-2q^{0}\kappa^{ij}+3(q^{0})^{2}\kappa^{i}_{m}\kappa^{mj}+{\cal
O}(\epsilon^{3}),$ (29) $\displaystyle G_{~{}~{}}$ $\displaystyle=$
$\displaystyle g~{}\\{1+2q^{0}\kappa+(q^{0})^{2}(\kappa^{2}+R)+{\cal
O}(\epsilon^{3})\\},$ (30) $\displaystyle G^{1/2}$ $\displaystyle=$
$\displaystyle g^{1/2}\\{1+q^{0}\kappa+\frac{1}{2}(q^{0})^{2}R+{\cal
O}(\epsilon^{3})\\},$ (31)

where $R=\kappa^{2}-\kappa_{ij}\kappa^{ij}$ is used.

Then the operator $\tilde{\Delta}^{(2)}$ can be expanded as follows

$\tilde{\Delta}^{(2)}=\Delta^{(2)}+q^{0}\hat{A}+(q^{0})^{2}\hat{B}+{\cal
O}(\epsilon^{3}),$ (32)

where,

$\hat{A}=-g^{-1/2}\frac{\partial}{\partial
q^{i}}g^{1/2}(2\kappa^{ij}\frac{\partial}{\partial
q^{j}}+g^{ij}\frac{\partial\kappa}{\partial q^{j}}),$ (33)
$\displaystyle\hat{B}$ $\displaystyle=$ $\displaystyle
g^{-1/2}\frac{\partial}{\partial
q^{i}}g^{1/2}(3\kappa^{im}\kappa_{m}^{j}\frac{\partial}{\partial q^{j}}$ (34)
$\displaystyle+$
$\displaystyle\frac{1}{2}g^{ij}\frac{\partial(\kappa^{2}-R)}{\partial
q^{j}}+2\kappa^{ij}\frac{\partial\kappa}{\partial q^{j}}).$

Then our two dimensional effective diffusion equation up to ${\cal
O}(\epsilon)$ is,

$\displaystyle\frac{\partial\phi^{(2)}}{\partial t}$ $\displaystyle=$
$\displaystyle D\Delta^{(2)}\phi^{(2)}$ (35) $\displaystyle+$ $\displaystyle
D\hat{A}\int_{-\epsilon/2}^{\epsilon/2}q^{0}\tilde{\phi}^{(3)}dq^{0}$
$\displaystyle+$ $\displaystyle
D\hat{B}\int_{-\epsilon/2}^{\epsilon/2}(q^{0})^{2}\tilde{\phi}^{(3)}dq^{0}+{\cal
O}(\epsilon^{3}).$

To proceed the $q^{0}$ integration, we suppose there is no diffusion flow in
normal direction in layer , that is,

$0=\frac{\partial\phi^{(3)}}{\partial q^{0}}=g^{1/2}\frac{\partial
G^{-1/2}\tilde{\phi}^{(3)}}{\partial q^{0}}.$ (36)

Solution is,

$\displaystyle\tilde{\phi}^{(3)}$ $\displaystyle=$
$\displaystyle\frac{1}{N}(G/g)^{1/2}\phi^{(2)}(q^{1},q^{2}),$ (37)
$\displaystyle N$ $\displaystyle\equiv$
$\displaystyle\int_{-\epsilon/2}^{\epsilon/2}(G/g)^{1/2}dq^{0}.$ (38)

Each integration can be explicitly performed, and we obtain

$\displaystyle N~{}~{}~{}~{}~{}$ $\displaystyle=$
$\displaystyle\epsilon+\frac{R}{24}\epsilon^{3}+{\cal O}(\epsilon^{5}),$ (39)
$\displaystyle<q^{0}>~{}$ $\displaystyle=$
$\displaystyle\frac{\kappa\epsilon^{2}}{12}+{\cal O}(\epsilon^{4}),$ (40)
$\displaystyle<(q^{0})^{2}>$ $\displaystyle=$
$\displaystyle\frac{\epsilon^{2}}{12}+{\cal O}(\epsilon^{4}),$ (41)

where we have used the definition

$<f(q^{0})>\equiv\frac{1}{N}\int_{-\epsilon/2}^{\epsilon/2}f(q^{0})(G/g)^{1/2}dq^{0}.$
(42)

We obtain the final form of equation up to ${\cal O}(\epsilon^{2})$ as

$\displaystyle\frac{\partial\phi^{(2)}}{\partial t}$ $\displaystyle=$
$\displaystyle
D\Delta^{(2)}\phi^{(2)}+\tilde{D}(\hat{A}\kappa+\hat{B})\phi^{(2)}$ (43)
$\displaystyle=$ $\displaystyle
D\Delta^{(2)}\phi^{(2)}+\tilde{D}g^{-1/2}\frac{\partial}{\partial
q^{i}}~{}g^{1/2}$ $\displaystyle\times$
$\displaystyle\\{(3\kappa^{im}\kappa_{m}^{j}-2\kappa\kappa^{ij})\frac{\partial}{\partial
q^{j}}-\frac{1}{2}g^{ij}\frac{\partial R}{\partial
q^{j}}\\}\phi^{(2)},~{}~{}~{}$

where $\tilde{D}=\frac{\epsilon^{2}}{12}D$.

We give two comments here. First, ${\cal O}(\epsilon)$ term disappears. Since
$\epsilon$ has the dimension of length, it always appears with curvature
$\kappa$. Therefore the 1st order term, if it exists, it contains 1st order of
curvature $\kappa$. But this curvature depends on unphysical choice of normal
unit vector $\vec{n}$, and so it does not appear.

Second, additional potential term disappears. In quantum mechanics, the
similar embedding techniques leads to the appearance of additional potential
term written by curvature. But in our classical case we have no such terms.
Because in diffusion equation, potential term breaks probability conservation
law, i.e.

$\partial\phi/\partial t=(D\Delta+V(x))\phi,$ $\to~{}\frac{d}{dt}\int
d^{3}x~{}\phi=\int d^{3}xV(x)\phi\neq 0.$

The normal diffusion flow can be written in general coordinate,

$J_{N}^{i}=-Dg^{ij}\frac{\partial\phi^{(2)}}{\partial q^{j}}.$ (44)

The anomalous diffusion flow is,

$J_{A}^{i}=-\tilde{D}\\{(3\kappa^{im}\kappa_{m}^{j}-2\kappa\kappa^{ij})\frac{\partial\phi^{(2)}}{\partial
q^{j}}-\frac{1}{2}g^{ij}\frac{\partial R}{\partial q^{j}}\phi^{(2)}\\}.$ (45)

The Diffusion equation is written as

$\displaystyle-\frac{\partial\phi^{(2)}}{\partial t}$ $\displaystyle=$
$\displaystyle\nabla_{i}(J_{N}^{i}+J_{A}^{i}),$ (46) $\displaystyle=$
$\displaystyle g^{-1/2}\frac{\partial}{\partial
q^{j}}~{}g^{1/2}(J_{N}^{i}+J_{A}^{i}),$

where $\nabla_{i}$ is the covariant derivative. By using a suitable boundary
condition, we can prove

$\frac{d}{dt}\int\phi^{(2)}g^{1/2}d^{2}q=0.$

##  IV Properties of curvature dependent flow

The anomalous flow equals to zero for the flat surface. The last term in
equation (45) shows that curvature gradient generate the flow without particle
density gradient. From the signature of this term, this flow goes from the
smaller Ricci scalar point to the larger Ricci scalar point. Ricci scalar can
take the negative, zero, and positive values. (Ricci scalar $R$ is related to
Gauss curvature by $R/2=\det[\kappa^{i}_{j}]$.) Let us work with the
coordinate which satisfies

$g_{ij}=\delta_{ij},~{}~{}\kappa^{i}_{j}=~{}\mbox{diag}[1/r_{1},~{}1/r_{2}],$

at point $P$, where $r_{i}$ is the curvature radius along the $q^{i}$
coordinate and it takes positive or negative value for convex or concave. (The
metric can be diagonalized by choosing the two coordinates as to satisfy
orthogonality, and it can be normalized by using the re-parametrization. The
second fundamental tensor is diagonalized by rotation of coordinate.)

Then at point $P$, we have $R=\frac{2}{r_{1}r_{2}}$ and,

  * •

$R>0$ if the surface is convex or concave.

  * •

$R=0$ if the surface is essentially flat.

  * •

$R<0$ if the surface is hyperbolic.

Therefore the flow goes from hyperbolic or flat points to convex or concave
points with positive larger Ricci scalar value.

Next we consider the first term in (45). We have positive or negative value
for

$f^{ij}\equiv 3\kappa^{im}\kappa_{m}^{j}-2\kappa\kappa^{ij},$

depending on the value of curvature. In our coordinate, we can immediately
write it in the simple form

$f^{ij}=\delta^{ij}(\frac{1}{(r_{i})^{2}}-\frac{2}{r_{1}r_{2}}).$ (47)

When the surface is hyperbolic ($R<0$),

$f^{11}=\frac{1}{(r_{1})^{2}}+\frac{2}{\mid
r_{1}r_{2}\mid}>0,~{}~{}f^{22}=\frac{1}{(r_{2})^{2}}+\frac{2}{\mid
r_{1}r_{2}\mid}>0,$

usual diffusion occurs (See fig. 2).

Figure 2: Wave packet on hyperbolic surface diffuses in two directions.

When the surface is convex or concave ($R>0$),

$f^{11}=\frac{1}{(r_{1})^{2}}-\frac{2}{\mid r_{1}r_{2}\mid}\\\ =\frac{\mid
r_{2}\mid-2\mid r_{1}\mid}{\mid r_{1}\mid^{2}\mid r_{2}\mid},$
$f^{22}=\frac{1}{(r_{2})^{2}}-\frac{2}{\mid r_{1}r_{2}\mid}\\\ =\frac{\mid
r_{1}\mid-2\mid r_{2}\mid}{\mid r_{2}\mid^{2}\mid r_{1}\mid}.$

In this case, we have three possibilities. One possibility is that both are
negative, if

$1/2<\mid\frac{r_{2}}{r_{1}}\mid<2.$

Then we have no diffusion but concentration occurs (See fig. 3).

Figure 3: Wave packet on convex surface concentrates.

The second possibility is that one of two is positive and the other is
negative, if

$\mid\frac{r_{2}}{r_{1}}\mid<1/2,~{}~{}\mbox{or}~{}\mid\frac{r_{2}}{r_{1}}\mid>2.$

Then we have diffusion in one direction, but concentration in another
direction (See fig. 4).

The third possibility for $R>0$ is,

$\mid\frac{r_{2}}{r_{1}}\mid=1/2,~{}~{}\mbox{or}~{}\mid\frac{r_{2}}{r_{1}}\mid=2.$

This is critical point, where the flow stops for larger curvature direction
and flow concentrates for smaller curvature direction.

Figure 4: Wave packet on convex surface with one direction curvature is over
two times higher than another. Packet diffuses in higher curvature direction
and concentrates in smaller curvature direction.

When the Ricci scalar is zero ($R=0$), for example $r_{2}=\infty$, $f^{22}=0$
and $f^{11}>0$ follows. The diffusion occurs only in $q^{1}$ direction but not
in another direction. (See fig. 5).

Figure 5: Wave packet on elliptic cylinder. Packet diffuses only in curved
direction but not in another direction.

In this way, this anomalous diffusion flow has much varieties depending on the
curvature.

##  V One Example: Elliptic Cylinder

Let us consider one simple example where diffusion coefficient depends on
curvature. The surface of elliptic cylinder is the case of $R=0$ just as
figure 5, but the surface has non zero mean curvature.

Ellipsoid is given by the equation

$(\frac{x}{a})^{2}+(\frac{y}{b})^{2}=1.$ (48)

Figure 6: Elliptic cylinder

Then any points on cylinder are specified by curved coordinate $\theta$ and
$z$;

$x=a\cos\theta,~{}~{}y=b\sin\theta.$ (49)

Another choice of coordinate instead of $\theta$ is,

$du=\sqrt{dx^{2}+dy^{2}}=f(\theta)~{}d\theta.$ (50)

where $f(\theta)$ is defined as

$f(\theta)\equiv\sqrt{a^{2}\sin^{2}\theta+b^{2}\cos^{2}\theta}.$ (51)

The length of $u$ is given by

$\displaystyle u(\phi)$ $\displaystyle=$
$\displaystyle\int_{0}^{\phi}f(\theta)d\theta$ (52) $\displaystyle=$
$\displaystyle b\int_{0}^{\phi}\sqrt{1-k^{2}\sin^{2}\theta}d\theta\equiv
bE(k,\phi),$

with $k=\sqrt{1-(a/b)^{2}},~{}a\leq b.$ $E(k,\phi)$ is the Elliptic integral
of the second kind. (See appendix)

The total length of $u$ is given by

$U\equiv 4bE(k,\pi/2).$

We use the normalized value for $u$ hereafter such that,

$\tilde{u}=u/U,~{}~{}0\leq\tilde{u}\leq 1.$ (53)

The normal unit vector $\vec{n}$ is given as

$\displaystyle\vec{n}$ $\displaystyle=$
$\displaystyle(\frac{x}{a^{2}\sqrt{(x^{2}/a^{4})+(y^{2}/b^{4})}},\frac{y}{b^{2}\sqrt{(x^{2}/a^{4})+(y^{2}/b^{4})}})$
(54) $\displaystyle=$
$\displaystyle\frac{1}{f(\theta)}(b\cos\theta,a\sin\theta).$

$\frac{\partial\vec{n}}{\partial\theta}=\frac{1}{f}(-b\sin\theta,a\cos\theta)-\frac{\partial_{\theta}f}{f}\vec{n}.$
(55)
$\vec{B}_{\theta}=\frac{\partial\vec{x}}{\partial\theta}=(-a\sin\theta,b\cos\theta).$
(56)

Then we obtain the second fundamental tensor.

$\kappa_{\theta\theta}=\vec{B}_{\theta}\cdot\frac{\partial\vec{n}}{\partial\theta}=\frac{ab}{f}.$
(57)

Then we collect all the necessary quantities as follows

$\displaystyle g_{\theta\theta}=f^{2},~{}~{}g_{zz}=1,~{}~{}g_{\theta z}=0,$
$\displaystyle\kappa_{\theta\theta}=\frac{ab}{f},~{}~{}\kappa_{zz}=\kappa_{\theta
z}=0,~{}~{}\kappa=\frac{ab}{f^{3}}.$ (58)

Then we obtain the total diffusion equation expressed by the parameters
$\theta$ and $z$.

$\frac{\partial\phi^{(2)}}{\partial
t}=(\frac{1}{f}\frac{\partial}{\partial\theta})D_{\theta}(\frac{1}{f}\frac{\partial}{\partial\theta})\phi^{(2)}+D\frac{\partial^{2}}{\partial
z^{2}}\phi^{(2)},$ (59)

where the effective diffusion coefficient depends on mean curvature.

$D_{\theta}=D(1+\frac{\epsilon^{2}\kappa^{2}}{12})=D(1+\varepsilon^{2}(b\kappa)^{2}),$

where $\varepsilon\equiv\epsilon/(2\sqrt{3}b)$.

Under this equation, we obtain the following particle number conservation law.

$\frac{d}{dt}\int
dz\int_{0}^{2\pi}d\theta~{}f(\theta)~{}\phi^{(2)}(\theta,z)~{}=0$

with suitable Neumann boundary condition.

By using the variable $\tilde{u}$ instead of $\theta$ we have simple
dimensionless equation,

$\frac{\partial\phi^{(2)}}{\partial\tau}=\frac{\partial}{\partial\tilde{u}}(1+V)\frac{\partial}{\partial\tilde{u}}\phi^{(2)}+\frac{\partial^{2}}{\partial\eta^{2}}\phi^{(2)}$
(60)

where
$\tau=tD/U^{2},~{}\eta=z/U,~{}V=\varepsilon^{2}(b\kappa)^{2},~{}U=4bE(k,\pi/2)$.

By using the approximation of elliptic function given in appendix, curvature
dependent potential $V$ can be written as function of $\tilde{u}$. The
simulation can be done as usual diffusion equation. For $0\leq\tilde{u}\leq 1$
and $0\leq\eta\leq 4$ using periodic boundary condition, since the length of
$\eta$ is larger than one of $\tilde{u}$, the diffusion in $u$ direction
occurs fastly and then the diffusion in $\eta$ direction follows slowly just
like $\phi\sim a+\sum_{k}b(k)e^{-k^{2}\tau}\cos(k\eta)$. The $u$-directional
diffusion can not occur uniformly, because at $\tilde{u}=0.25,$ and $0.75$ the
diffusive coefficient is higher than other points. Therefore the slope of
diffusion field is small especially at these two points during the diffusion
process. (fig. 7)

Figure 7: Snap shot of diffusion process starting from the wave packet
$\phi=\sin^{10}(\pi\tilde{u})$ as the initial condition. At two points (0.25
and 0.75), diffusion occurs quickly and its slope is smaller than other.

Figure 8: Mean curvature as a function of $\tilde{u}$ when $b/a=2$

##  VI Conclusion

We have discussed on the diffusion equation on curved surface embedded in
$R_{3}$. We obtained the new diffusion equation up to ${\cal O}(\epsilon^{2})$
which includes anomalous diffusive flow additional to the usual one. This
anomalous flow depends on the second fundamental tensor, and it has not only
diffusion but also concentration properties depending on the curvature of its
manifold.

At the point with negative Ricci scalar $R<0$, surface is hyperbolic, and
diffusion in both direction occurs. (fig.2)

When Ricci scalar is positive $R>0$, we have three possibilities. $r_{i}$
appearing below is curvature radius in each direction ($i=1,2$).

  * •

Concentration in both direction (fig.3), when

$1/2<\mid\frac{r_{2}}{r_{1}}\mid<2.$

  * •

Concentration in smaller curvature direction, and diffusion in higher
curvature direction (fig.4), when

$\mid\frac{r_{2}}{r_{1}}\mid<1/2,~{}~{}\mbox{or}~{}\mid\frac{r_{2}}{r_{1}}\mid>2.$

  * •

Concentration in smaller curvature direction, and no flow in higher curvature
direction, when

$\mid\frac{r_{2}}{r_{1}}\mid=1/2,~{}~{}\mbox{or}~{}\mid\frac{r_{2}}{r_{1}}\mid=2.$

When Ricci scalar is zero $R=0$, surface is essentially flat, but we have
finite curvature radius in one direction. Then we have diffusion only in this
direction (fig.5).

In the case of surface of elliptic cylinder, we gave a concrete form of
equation and we showed the curvature dependent diffusion coefficient.

$D_{u}=D(1+\frac{\epsilon^{2}\kappa^{2}}{12}),~{}~{}~{}D_{z}=D,$

where $\kappa$ is the mean curvature. In this case curvature dependence is
simply included into diffusion coefficient. However this is not true in
general case, where situation is much more complicated, and this can be seen
from the form of anomalous flow.

The application to the pattern formation by reaction diffusion using this
obtained equation is not yet finished. This will be done in further
publication.

##  VII Appendix

We approximate the elliptic integral of the second kind.

$E(k,\phi)\equiv\int_{0}^{\phi}\sqrt{1-k^{2}\sin^{2}\theta}~{}d\theta$ (61)

with

$k=\sqrt{1-(a/b)^{2}}.$

Under the expansion in powers of $k^{2}$, we obtain the power series of
Elliptic integral of the second kind.

$E(k,\phi)=\phi-\sum_{n=1}^{\infty}\frac{k^{2n}(2n-3)!!}{n!~{}2^{n}}\int_{0}^{\phi}\sin^{2n}\theta
d\theta.$ (62)

Since the integration part can be expanded by $\phi$ and $\sin 2n\phi$, we
have

$E(k,\phi)=a_{0}\phi+\sum_{n=1}^{\infty}a_{n}\sin 2n\phi.$ (63)

with the relation

$\displaystyle a_{0}$ $\displaystyle=$ $\displaystyle\frac{2}{\pi}E(k,\pi/2),$
(64) $\displaystyle a_{n}$ $\displaystyle=$
$\displaystyle(-1)^{n}\frac{2E(k,\pi/2)}{n\pi}$ (65)
$\displaystyle+\frac{4}{\pi}\int_{0}^{\pi/2}\sin
2n\phi~{}E(k,\phi)d\phi.~{}~{}(n\geq 1)$

For the real Char fishes, $b/a$ takes values $1.5\sim 2.5$. Then the value of
$k$ takes $0.75\sim 0.92$. When $b/a=2$, each values of $a_{n}$ is given
numerically

$a_{0}=0.771,~{}~{}a_{1}=0.123,~{}~{}a_{2}=-0.00506,~{}~{}a_{3}=0.000558.$

Now we have for $u$ given in (52),

$u/b=E(k,\phi)=a_{0}\phi+a_{1}\sin 2\phi+a_{2}\sin 4\phi+\cdots.$ (66)

And we rewrite it by using normalized $u$,

$\phi=2\pi\tilde{u}-\frac{a_{1}}{a_{0}}\sin 2\phi-\frac{a_{2}}{a_{0}}\sin
4\phi-\cdots,$ (67)

where $\tilde{u}=u/(4bE(k,\pi/2)).$

The iteration method up to order $(a_{1}/a_{0})^{1}$ gives

$\phi=2\pi\tilde{u}-\frac{a_{1}}{a_{0}}\sin(4\pi\tilde{u}).$ (68)

Then we take the derivative by $u$ in both hand sides.

$\frac{1}{f}=\frac{1}{ba_{0}}(1-\frac{2a_{1}}{a_{0}}\cos(4\pi\tilde{u})),$
(69)

where the following relation is used.

$\frac{du}{d\phi}=f(\phi)\equiv\sqrt{a^{2}\sin^{2}\phi+b^{2}\cos^{2}\phi}.$

Then the mean curvature given by (58) is obtained as function of $u$.

$b\kappa=\frac{ab^{2}}{f^{3}}=\frac{1}{\beta(a_{0})^{3}}(1-\frac{2a_{1}}{a_{0}}\cos(4\pi\tilde{u}))^{3},$
(70)

where $\beta=b/a$. This function is shown in figure 8.

## References

  * (1) J. Faraudo, J. Chem. Phys, 116 (2002) 5831-5841; J. Balakrishnan, arXiv:physics/0308089, 25 Aug 2003.
  * (2) A. Naji and F. Brown, J. Chem. Phys. 126 (2007) 235103; E. Reister and U. Seifert, arXive:cond-mat/0503568, 23 Mar 2005.
  * (3) R. C. T. da Costa, Phys. Rev. 23 (1981) 1982; J. Tolar, 1988 Lecture Notes in Physics 313, ed. H. D. Doever, J. D. Henning and T. D. Raev, (Springer-Verlag, Berlin, Heidelberg) 268.
  * (4) N. Ogawa, K. Fujii, and K. P. Kobushkin, Prog. Theor. Phys. 83 (1990) 894; N. Ogawa, K. Fujii, N. M. Chepilko, and K. P. Kobushkin, Prog. Theor. Phys. 85 (1991) 1189; N. Ogawa, Prog. Theor. Phys. 87 (1992) 513.
  * (5) K. Fujii and N. Ogawa, Prog. Theor. Phys. 89 (1993) 575.
  * (6) A. M. Turing, Phil. Trans. R. Soc. London B 237 (1952) 37-72; H. Meinhardt, Models of Biological Pattern Formation., Academic Press, London (1982); J. D. Murray, Mathematical Biology. 2nd ed. Springer, New York (1989).
  * (7) H. Shoji, Y. Iwasa and S. Kondo, J. Theor. Biol. 224 (2003) 339-350; H. Shoji and Y. Iwasa, J. Theor. Biol. 237 (2005) 104-116.


11footnotetext: See Supplemental Material

# Electronic structure of the frustrated diamond lattice magnet NiRh2O4

B. Zager  Department of Physics, Brown University, Providence, Rhode Island
02912, United States  J. R. Chamorro  Department of Chemistry, Johns Hopkins
University, Baltimore, MD, USA  Institute for Quantum Matter, William H.
Miller III Department of Physics and Astronomy, Johns Hopkins University,
Baltimore, MD, USA  L. Ge  School of Physics, Georgia Institute of Technology,
Atlanta, Georgia 30332, USA  V. Bisogni  National Synchrotron Light Source II,
Brookhaven National Laboratory, Upton, New York 11973, USA  J. Pelliciari
National Synchrotron Light Source II, Brookhaven National Laboratory, Upton,
New York 11973, USA  J. Li  National Synchrotron Light Source II, Brookhaven
National Laboratory, Upton, New York 11973, USA  G. Fabbris  Advanced Photon
Source, Argonne National Laboratory, Lemont, Illinois 60439, USA  T. M.
McQueen  Department of Chemistry, Johns Hopkins University, Baltimore, MD, USA
Institute for Quantum Matter, William H. Miller III Department of Physics and
Astronomy, Johns Hopkins University, Baltimore, MD, USA  Department of
Materials Science and Engineering, Johns Hopkins University, Baltimore, MD,
USA  M. Mourigal  School of Physics, Georgia Institute of Technology, Atlanta,
Georgia 30332, USA  K. W. Plumb  Department of Physics, Brown University,
Providence, Rhode Island 02912, United States

###### Abstract

The $A$-site spinel NiRh2O4 is the only known realization of a spin-1 diamond
lattice magnet and is predicted to host unconventional magnetic phenomena
driven by frustrated nearest and next-nearest neighbor exchange as well as
orbital degeneracy. Previous works found no sign of magnetic order but found a
gapped dispersive magnetic excitation indicating a possible valence bond
magnetic ground state. However, the presence of many competing low energy
degrees of freedom and limited empirical microscopic constraints complicates
further analysis. Here, we carry out resonant inelastic x-ray scattering
(RIXS) and x-ray absorption spectroscopy (XAS) to characterize the local
electronic structure of NiRh2O4. The RIXS data can be partly described by a
single-ion model for tetrahedrally coordinated Ni2+ and indicates a tetragonal
distortion $\Delta t_{2}\\!=\\!70$ meV that splits the $t_{2}$ orbitals into a
high energy orbital singlet and lower energy orbital doublet. We identify
features of the RIXS spectra that are consistent with a Rh-Ni two-site
excitation indicating strong metal-metal hybridization mediated by oxygen in
NiRh2O4. We also identify signatures of electron-phonon coupling through the
appearance of phonon sidebands that dress crystal field excitations. These
results establish the key energy scales relevant to the magnetism in NiRh2O4
and further demonstrate that covalency and lattice dynamics play essential
roles in controlling the magnetic ground states of $A$-site spinels.

##  I Introduction

Materials with frustrated magnetic interactions provide a platform for
realizing novel phases which avoid conventional symmetry-breaking order [1].
Such phases are acutely sensitive to a hierarchy of competing energy scales
involving spin, charge, orbital, and lattice degrees of freedom [2]. This
enables precise tuning of the collective orders in these systems to explore
new phenomena and develop new technologies [3]. Along this line, there has
recently been a significant interest in novel phases driven by strong spin
orbit coupling (SOC) that has motivated significant work on $4d$/$5d$
transition metal compounds [4]. In some cases, however, the relatively weak
SOC in $3d$ compounds may still become relevant. In particular, when there is
an orbital degeneracy, SOC may compete with the Jahn-Teller (JT) effect to
produce a spin-orbital-lattice entangled state [5].

The spinel structure comprises an intensely studied class of materials with
both fundamental significance and widespread applications. The general formula
$AB_{2}X_{4}$ contains tetrahedrally coordinated $A^{2+}$ ions and
octahedrally coordinated $B^{3+}$ ions which occupy diamond and pyrochlore
sublattices respectively. In the case of magnetic $B$ ions, the system forms a
prototypical geometrically frustrated pyrochlore magnet. On the other hand,
for magnetic $A$ site ions, the bipartite diamond lattice may be frustrated in
the presence of competing nearest and next-nearest neighbor antiferromagnetic
exchange [6]. Previously studied diamond lattice magnets include the $A$-site
spinels Co3O4 [7], $A$Al2O4 ($A=\;$Mn, Fe, Co) [8, 9], $A$Sc2S4 ($A=\;$Mn, Fe)
[10, 11], MnSc2Se4 [12], and $A$Rh2O4 ($A=\;$Co, Cu, Ni) [13, 14], as well as
the lanthanides LiYbO2 [15] and NaCeO2 [16]. These materials host a variety of
magnetic phenomena, ranging from long-range ordered states to disordered spin
liquid and spin glass states. In many cases, especially in the presence of an
orbital degree of freedom, these materials lie near a quantum critical point,
with multiple competing phases exhibiting strong sensitivity to disorder and
external perturbations [17, 18, 19, 20, 21, 22, 23, 24, 25]. It is thus
essential to characterize the microscopic energy scales associated with the
spin, lattice, and orbital degrees of freedom in these materials.

Another topic of recent interest is frustrated magnetism in spin-1 materials
[26, 27, 28], which support more degrees of freedom than spin-$\frac{1}{2}$,
while still being sensitive to quantum fluctuations. In particular, recent
studies have predicted spin-1 diamond lattice antiferromagnets to host
unconventional magnetic phenomena, namely topological paramagnetism [29],
spiral spin liquid phases [30], quantum critical phenomena [31], and excitonic
magnetism [32].

The spinel NiRh2O4 contains spin-1 Ni2+ ions on the $A$ sites and nonmagnetic
Rh3+ ions on the $B$ sites, thus realizing the only known spin-1 on a diamond
lattice. Early studies identified a cubic to tetragonal structural transition
[33] and an apparent antiferromagnetic transition [34]. However, a recent
study on high quality samples found no sign of magnetic order down to 0.1 K
[14], implying that chemical disorder may have stabilized ordering found in
the original studies. In that study, specific heat and x-ray diffraction
measurements found a small entropy loss associated with the structural
distortion at $T=440$ K, indicating only a partial lifting of orbital
degeneracy. This finding is consistent with the $p_{\text{eff}}=3.3\mu_{B}$
paramagnetic moment that is significantly larger than the pure spin-1 value of
$2.83\mu_{B}$, implying an orbital contribution to the magnetism. The Curie-
Weiss temperature of $\Theta_{\text{CW}}=-11$ K indicates a large frustration
parameter $f=\Theta_{\text{CW}}/T_{N}>100$. Inelastic neutron scattering (INS)
found gapped dispersive magnetic excitations suggestive of a valence bond
ground state. However, an incomplete knowledge of the electronic ground state
has limited modeling of the INS data. Recent theoretical work predicted that
spin-orbit coupling, crystal fields, and correlations generate a spin-orbital
singlet ground state in NiRh2O4 [32, 35]. Such a state would explain the
absence of magnetic ordering and key features of the magnetization, specific
heat, and neutron data. But measurements of the high energy excitations
characterizing the orbital configurations of NiRh2O4 are required to confirm
this picture.

In this paper, we use a combination of resonant inelastic x-ray (RIXS), x-ray
absorption spectroscopy (XAS), and inelastic neutron scattering to provide a
detailed account of the electronic structure of NiRh2O4. We observe crystal
field excitations for Ni2+ in a distorted tetrahedral environment that are
well-described by a single-ion model including Coulomb interaction, spin-orbit
coupling, crystal field splitting. We find additional electronic excitations
corresponding to two-site Rh to Ni charge transfer. We also find that the
orbital excitations are coupled to optical phonons giving rise to distinct
phonon sidebands that dress $d$-$d$ excitations. These results provide a
detailed description of the local electronic structure of a novel frustrated
magnet along with the key energy scales required for understanding the
magnetic ground state and low-energy excitations. We have additionally shown
that both covalency and lattice dynamics play essential roles in this material
and should be considered in any realistic model for the magnetism. These
insights provide guidance for exploring novel effects in NiRh2O4 and other
magnetically frustrated spinels using pressure, magetic fields, and chemical
substitution.

This paper is organized as follows. In section II, we describe the
experimental details. In section III, we present the RIXS, XAS, and INS
results and describe our attempts to model the data. We first use a single ion
model for Ni2+ in a distorted tetrahedral crystal field, then a minimal two-
site Ni-Rh hopping model, and finally we incorporate electron phonon coupling
in order to capture the lineshapes of the RIXS spectra. Our findings are
summarized in section IV.

Figure 1: Overview of RIXS data measured at $T\\!=\\!40$ K and
$2\theta\\!=\\!90$°. (a) Local coordination of Ni and Rh in NiRh2O4. Rh is
octahedrally coordinated while Ni2+ $3d$ orbitals are split by tetrahedral
crystal field with tetragonal distortion $\Delta t_{2}\\!>\\!0$. (b) Measured
RIXS intensity vs. energy loss and incident energy measured in $\pi$
polarization. The XAS measured in total fluorescence yield is plotted in
white. (c) RIXS spectrum in both $\pi$ and $\sigma$ polarization at
$E_{in}\\!=\\!853.4$ eV, indicated by the dashed line in (b). The inset shows
the low-energy part of the spectrum and an empirical fit.

##  II Experiment

Polycrystalline NiRh2O4 was synthesized following the methods described in
[14]. The samples contain 4% nonmagnetic Rh2O3 by mass, which is not expected
to contribute any observable effects to our measurements.

Resonant inelastic x-ray scattering (RIXS) measurements at the Ni $L_{3}$ edge
were performed at the soft inelastic x-ray scattering (SIX) beamline 2-ID at
the National Synchrotron Light Source II (NSLS-II) at Brookhaven National
Laboratory. Measurements were carried out at 40 K and 25 K using an incident
x-ray polarization in both linear horizontal ($\pi$) and linear vertical
($\sigma$) geometries. The scattering angle was fixed at 90° to minimize the
contribution from Thomson elastic scattering and giving a momentum transfer of
$Q\\!=\\!0.61$ Å-1. The incident x-ray energy was varied across the Ni $L_{3}$
edge ($\sim$853 eV). The combined energy resolution was determined to be 31
meV based on the full width at half maximum (FWHM) of the elastic signal from
carbon tape.

X-ray absorption spectroscopy (XAS) measurements at the O $K$ edge were
performed at NSLS-II beamline 2-ID (SIX), measured in fluorescence yield (FY).
XAS measurements at the Rh $L$ edges at Advanced Photon Source 4-ID-D,
measured in total electron yield (TEY).

Inelastic neutron scattering (INS) measurements were performed on the Fine
Resolution Fermi chopper spectrometer (SEQUOIA) at the Spallation Neutron
Source, Oak Ridge National lab. The incident neutron energy was fixed at 160
meV using the coarse resolution chopper (FC2) rotating at 600 Hz. The same 4 g
sample used in [14] was loaded in an aluminum can, and held at 3.6 K for the
measurement. Scattering contributions from the sample environment were removed
by subtraction of an empty can data set during data reduction.

##  III Results

Fig. 1(b) shows the RIXS intensity as a function of incident energy
$E_{\text{in}}$ and energy loss $\Delta E$ measured at 40 K with incident
$\pi$ polarization. The XAS, shown in white, contains a weak pre-edge feature
at 852.9 eV, a main peak split into two features at 853.4 eV and 853.8 eV, and
a satellite peak at 856.4 eV. The RIXS spectrum shows an elastic line at
$\Delta E\\!=\\!0$ eV, five Raman-like features at 0.065 eV, 0.25 eV, 0.5 eV,
1.1 eV, and 1.6 eV, a broad charge-transfer (CT) background between 2 and 4
eV, and a fluorescence line (FL) at constant scattered photon energy.

Fig. 1(b) shows RIXS scans at the main resonance of
$E_{\text{in}}\\!=\\!853.4$ eV, indicated by the dotted line in Fig. 1(a), in
both $\pi$ and $\sigma$ incident polarization at 40 K. The inset shows the low
energy region near the elastic line. The elastic intensity due to Thomson
scattering is expected to be strongly suppressed in the $\pi$-polarized data
for the 90$\degree$ scattering angle configuration of these measurements.
However, the $\pi$ and $\sigma$-polarized data show similar intensity around
$\Delta E\\!=\\!0$, suggesting unresolved low energy excitations. The elastic
line cannot be fit by a single resolution-limited Voigt function, further
indicating the presence of unresolved low energy contributions. This can be
explained by the $\sim$12 meV dispersive magnetic excitation observed in
inelastic neutron scattering [14]. By including an additional resolution-
limited Voigt function at 12 meV, we obtain an adequate fit to the quasi-
elastic line, as shown in the inset of Fig. 1(b) for the $\pi$-polarized data,
which also includes contributions from two overlapping higher energy peaks. A
similar fit can be obtained for the $\sigma$-polarized data, with a slightly
larger contribution from the elastic peak likely originating from Thomson
scattering that contributes in that geometry.

###  III.1 Single-ion model

Figure 2: Comparison of the two crystal field schemes used for modeling the
RIXS data. Scenario 1 ($\Delta t_{2}\\!>\\!0$) is shown in (a) and scenario 2
($\Delta t_{2}\\!<\\!0$) is shown in (b). The pink markers show the
eigenvalues. The orange markers in the inset indicate which eigenvalues are
dipole-allowed transitions from the ground state. (c) Comparison of the two
schemes for modeling the Ni $L_{3}$ XAS (TFY). The calculated XAS spectra are
normalized to the intensity at the maximum of the measured spectrum.

To identify the features in the measured spectra, we calculate the RIXS cross
section for a single Ni2+ ion including Coulomb interaction, spin-orbit
coupling (SOC), and crystal field (CF) splitting. Calculations were
implemented using the EDRIXS package [36], further details are described in
the Supplemental Material [37]. We include onsite Coulomb interactions as
parameterized by the Slater integrals: $F_{dd}^{2}$, and $F_{dd}^{4}$ describe
the direct Coulomb repulsion between $d$ electrons, $F_{pd}^{2}$ describes
direct Coulomb repulsion between $d$ electrons and the $2p$ core hole, and
$G_{pd}^{1}$ and $G_{pd}^{3}$ describe the Coulomb exchange between $d$
electrons and the $2p$ core hole. Initial values for the above parameters are
calculated for a free Ni2+ ion in $3d^{8}$ and $2p^{5}3d^{9}$ configurations
by the Hartree-Fock method using Cowan’s code [38]. To account for the reduced
intra-atomic repulsion due to covalency effects in the solid, we include
empirical scale factors on $F_{dd}$, $F_{pd}$, and $G_{pd}$ that are
determined by fitting to the RIXS data.

The crystal-field is parameterized by the tetrahedral splitting $10Dq\\!<\\!0$
and the splitting due to the tetragonal distortion $\Delta
t_{2}\\!=\\!E_{b_{2}}\\!-\\!E_{e^{\prime}}$, $\Delta
e\\!=\\!E_{b_{1}}\\!-\\!E_{a_{1}}$. The tetrahedral crystal field ($T_{d}$
point group) splits the $d^{8}$ configuration into $e^{4}t_{2}^{4}$. The
tetragonal distortion ($T_{d}\rightarrow D_{2d}$) further splits the $e$
doublet into $a_{1}$ ($d_{z^{2}}$) and $b_{1}$ ($d_{x^{2}-y^{2}}$) singlets,
and the $t_{2}$ triplet into a $b_{2}$ ($d_{xy}$) singlet and $e^{\prime}$
($d_{xy}$, $d_{yz}$) doublet, as shown in Fig. 1(a). These symmetry
considerations leave the sign and magnitude of the splittings unconstrained.
There have been no reported direct measurements of the crystal field splitting
in NiRh2O4, and there are conflicting reports of the sign of the tetragonal
splitting. This splitting cannot be constrained only from knowledge of the
crystal structure because expectations from electrostatic considerations are
often inaccurate due to additional effects from covalency and spin-orbit
coupling [2]. As we will show below, both effects are significant in NiRh2O4.

The low temperature structure reported in [14] is tetragonal with elongated
NiO4 tetrahedra and compressed RhO6 octahedra. This led the authors to propose
a crystal field scheme with the $t_{2}$ levels split into a lower $b_{2}$
state and upper $e^{\prime}$ states ($\Delta t_{2}\\!<\\!0$), following the
model from [39], as expected from electrostatic considerations. However, the
DFT calculation in [35] suggests that the $t_{2}$ levels are split into lower
$e^{\prime}$ states and an upper $b_{2}$ state ($\Delta t_{2}\\!>\\!0$). We
find that although a single ion model is not sufficient to explain the full
RIXS spectra for NiRh2O4, it can capture many of the essential features, and
enables us to distinguish between the two possible scenarios for the
tetragonal splitting and place constraints on its magnitude.

We constrain our model to maintain consistency with material trends for
insulating Ni compounds and with other spectroscopic measurements on NiRh2O4
[39, 40]. The peaks around 0.5 eV and 1 eV fix the value of $10Dq\approx 0.55$
eV, as they correspond to $e^{4}t_{2}^{4}\rightarrow e^{3}t_{2}^{5}$ and
$e^{4}t_{2}^{4}\rightarrow e^{2}t_{2}^{6}$ transitions. This $10Dq$ value is
consistent with DFT [35] and the values for tetrahedrally coordinated Ni2+ in
other compounds [39]. We can also constrain the value of atomic spin-orbit
coupling $\lambda$ by the requirement of a dipole-allowed level near 11 meV to
agree with the neutron scattering data. Finally, the peak at 1600 meV
corresponds to an excitation within the Hund’s multiplet and fixes the value
of $F_{dd}$. To match the energy of the highest energy peak, $F_{dd}$ must be
set to 0.5, giving a Hund’s coupling of
$J_{H}=\frac{1}{14}(F_{dd}^{2}+F_{dd}^{4})=0.71$ eV. This parameter also sets
the energy of intra-$t_{2}$ spin-flip excitations ($S=1\rightarrow S=0$)
between 1.1 and 1.3 eV. These excitations have minimal intensity in the
computed RIXS cross section and we do not expect to observe them above the
signal-to-noise of our data, see Fig. 2(a) and 2(b). The 50% reduction in the
intra-atomic Coulomb interaction $F_{dd}$, compared with the atomic values,
suggests strong Ni-O covalency [41] and already indicates the inadequacy of a
single ion model for NiRh2O4, we address this later in section III.2. The
Slater integrals $F_{pd}$ and $G_{pd}$ determine the intermediate state
energies and do not affect the RIXS peak energies. However, they do influence
the RIXS intensities and energy splitting of the main XAS peak.

Table 1: Parameters for the single-ion model using $\Delta t_{2}>0$ and $\Delta t_{2}<0$. $10Dq$, $\Delta t_{2}$, $\Delta e$, and $\lambda$ are in meV. $F_{dd}$, $F_{pd}$, and $G_{pd}$ are dimensionless. | $10Dq$ | $\Delta t_{2}$ | $\Delta e$ | $\lambda$ | $F_{dd}$ | $F_{pd}$ | $G_{pd}$
---|---|---|---|---|---|---|---
$\Delta t_{2}>0$ | -580 | 70 | 56 | 13 | 0.5 | 0.7 | 0.75
$\Delta t_{2}<0$ | -530 | -50 | -40 | 27 | 0.5 | 0.7 | 0.75

We have tested models for tetrahedral compression and elongation and find that
they both capture many features of the data; table 1 shows the best parameters
for the two scenarios (1) $\Delta t_{2}\\!>\\!0$ and (2) $\Delta
t_{2}\\!<\\!0$, plotted in Fig. 2. $\Delta e$ is underdetermined by our data
and mainly contributes to the splitting of the peaks above 500 meV. Here, we
assume that $\Delta e$ has the same sign as $\Delta t_{2}$. We also expect
$|\Delta e|<|\Delta t_{2}|$, due to the reduced Ni-O hybridization of the $e$
states compared to the $t_{2}$ states in tetrahedral symmetry. We find that
the constraint $\Delta e=0.8\Delta t_{2}$ provides an adequate agreement with
the data as shown in Fig. 2. Both models reproduce the observed peaks at 65
meV, 0.5 eV, and 1.5 eV and while scenario 1 better reproduces the small
dichroism of the 0.5 eV feature, scenario 2 more faithfully reproduces the
negligible dichroism of the 1 eV and 1.6 eV peaks. The high energy
excitations, $\Delta E\\!>\\!500$ meV, are far above the insulating gap in
NiRh2O4 and are thus coupled to delocalized states. This effect is not
captured in the single-ion model and may explain the Fano-like lineshape.
However, we found that a more careful examination of the low energy
excitations for $E<100$ meV enables a distinction between crystal field
models.

The low energy part of the spectrum is shown in the insets of Fig. 2(a) and
2(b), with the orange markers indicating energy levels with nonvanishing
neutron cross section. Since photons emitted with energies close to the
absorption edge are more likely to be reabsorbed, we expect strong self
absorption effects near the elastic line and our model should predict a larger
quasi-elastic intensity than what is observed. The observed intensities are
also likely modified by intersite magnetic exchange interactions [35] that are
not included in our single site model. Nevertheless, a careful comparison of
the observed linear dichroism at low energy transfers reveals that scenario 1,
$\Delta t_{2}\\!>\\!0$, is more consistent with our data. The low energy
subspace of scenario 1 is equivalent to the model from [35], and yields
excited states at 3∗, 11(2)∗, 22(2), and 65(3)∗ meV. The low energy subspace
of scenario 2 is equivalent to the model from [39], yielding excited states at
11(2)∗, 38∗, 55, 65(2)∗, 72, and 75 meV. Asterisks indicate those states with
nonvanishing neutron intensity and parentheses indicate the state degeneracy,
disregarding any splitting of less than 1 meV. The calculation for scenario 1
shows that the dominant contribution to the quasi-elastic RIXS intensity comes
from the ground state and a 3 meV excitation for both polarizations [Fig.
2(a)]. This is consistent with the nearly equivalent quasielastic lines we
measured for each polarization. However, for scenario 2, the dominant
contribution comes from the ground state for $\sigma$ polarization and from
the 11 meV excitation for $\pi$ polarization [Fig. 2(b)]. Although the
energies of these excitations fall within our experimental resolution,
scenario 2 should result in a more pronounced difference in the quasielastic
lineshape between $\sigma$ and $\pi$ polarizations and that is not consistent
with our data.

To provide an additional check for our single ion model, we compare with the
measured x-ray absorption spectrum (XAS) at the Ni $L_{3}$ edge measured in
total fluorescence yield (TFY) shown in Fig. 2(c). The main XAS peak is split
into lower and upper peaks at 853.4 and 853.8 eV, corresponding to the states
$2p^{5}e^{4}t_{2}^{5}$ and $2p^{5}e^{3}t_{2}^{6}$. In the single-ion model,
this splitting depends on $10Dq$, $F_{pd}$, and $G_{pd}$. Although the XAS
lineshape is known to be distorted by TFY measurements, the calculated
spectrum shows better qualitative agreement with our data for scenario 1. In
particular, a model with $\Delta t_{2}>0$ more faithfully reproduces the
relative intensities of the 853.4 and 853.8 eV peaks, and captures additional
observed intensity at 854.7 eV that is not predicted by a model with $\Delta
t_{2}<0$ [Fig 2 (c)]. However, we find that for either scenario, the single-
ion model cannot capture the large, $\sim 3$ eV, splitting between the main
peak and satellite peak at 856.4 eV.

Based on the above considerations, we conclude that scenario 1, splitting the
$t_{2}$ levels into lower energy $d_{xz},d_{yz}$ orbitals and a higher energy
$d_{xy}$ orbital, is the correct crystal field scheme in NiRh2O4. However,
there are many notable discrepancies between the single ion model and our
data. First, a single-ion model completely fails to reproduce the intense 250
meV RIXS peaks for any reasonable set of parameters [37], and second, it does
not accurately capture the broad asymmetric lineshape of the 0.5 eV peak. Both
of these features in the RIXS spectra arise because of two distinct effects
that cannot be accounted for in a single-ion description. First, there is
strong Ni-Rh hybridization, and second, there is strong electron-phonon
coupling. The evidence for each of these effects and a detailed discussion of
their respective influence on the RIXS spectra and magnetism in NiRh2O4 is
discussed below in sections III.2 and III.3.

###  III.2 Ni -Rh hybridization

Figure 3: (a) XAS at the O $K$ edge for NiRh2O4 (TFY), with NiO (TFY) and
ZnRh2O4 (TEY) for comparison (NiO data from [42] and ZnRh2O4 data from [43])
(b) XAS measured at the Rh $L_{3}$ edge and (c) the Rh $L_{2}$ edge.

We now extend our model to include the influence of electronic hybridization
between Ni and Rh sites. Such metal-metal hybridization is supported both by
experiment and ab-initio calculations. Density functional theory calculations
[35] have found strong Ni-Rh hybridization is mediated by the intermediate
oxygens. Furthermore, the calculation found an insulating gap of 250 meV, with
the highest occupied and lowest unoccupied states having mostly Rh and Ni
character respectively. This suggests that the lowest interband transition
consists of Rh$\rightarrow$Ni excitations.

We have carried out x-ray absorption (XAS) measurements at the O $K$ and Rh
$L$ edges in order to characterize the electronic states at the Rh and O sites
in NiRh2O4. Fig. 3(a) shows the O $K$ edge XAS for NiRh2O4 measured in total
fluorescence yield (TFY), with the spectra for NiO (TFY) and ZnRh2O4 (TEY)
shown for comparison. NiO provides a comparison to bonding in an NiO6
octahedron while ZnRh2O4, being isostructural to the cubic phase of NiRh2O4,
provides a comparison to bonding in a RhO6 octahedron in the absence of an
unfilled neighboring $3d$ shell.

The pre-edge region contains a small peak at 530.5 eV and a large peak at 533
eV, corresponding to O $2p$ states hybridized with empty metal $d$ states. The
region above 535 eV corresponds to Ni $4sp$ and Rh $5sp$ states. By comparing
to the projected density of states from DFT calculations [35], we can assign
the small pre-edge peak at 530.5 eV to the unfilled Ni $t_{2}$ states and the
large peak at 533 eV to the unfilled Rh $e_{g}$ states. The pre-edge peak
intensity is determined by both the number of empty metal states and the
degree of hybridization [44]. In NiO and ZnRh2O4, where the metal sites
provide two empty states per O site, the pre-edge peak intensities are
comparable. In NiRh2O4, there are 2 empty Rh states and 0.5 empty Ni states
per O site. This small increase in the number of available states alone cannot
explain the significant enhancement of pre-edge peak intensity. The intensity
can thus be explained by an increased hybridization due to the cooperative
influence of the Ni and Rh. This can be expected based on the large inductive
effect of Rh3+ [45]. It may also be enhanced due to the long range exchange
interactions. In the five-site exchange pathway $A$-O-$B$-O-$A$ via a
nonmagnetic $B$ cation, the dominant contribution is thought to involve the
empty states at the $B$ sites [46, 47, 48].

Table 2: Table of peak positions and widths in eV extracted from the Rh $L$ edge XAS by fitting a Lorentzian + arctangent lineshape. | $E_{L3}$ | $E_{L2}$ | $\Gamma_{L3}$ | $\Gamma_{L2}$
---|---|---|---|---
Rh acetate | 3005.5(1) | 3147.7(1) | 2.0(1) | 1.8(2)
RhCl${}_{3}\negthinspace\cdot\negthinspace x$H2O | 3005.7(1) | 3148.0(1) | 1.8(1) | 1.6(2)
NiRh2O4 | 3006.1(1) | 3148.3(1) | 2.0(1) | 1.7(2)
Rh2O${}_{3}\negthinspace\cdot\negthinspace 5$H2O | 3006.3(1) | 3148.5(1) | 1.9(1) | 1.7(2)
Ag3LiRh2O6 | 3007.0(1) | 3149.4(1) | 3.0(2) | 2.4(3)

Fig. 3(b) and Fig. 3(c) show the XAS at the Rh $L_{3}$ and $L_{2}$ edges
respectively, for NiRh2O4 and reference samples with known valence: Rh acetate
(2+), RhCl${}_{3}\negthinspace\cdot\negthinspace x$H2O (3+),
Rh2O${}_{3}\negthinspace\cdot\negthinspace 5$H2O (3+), Ag3LiRh2O6 (4+) [49].
The sharp white line peaks correspond to transitions from $2p$ core levels to
empty $4d$ states. The Rh2+ and Rh3+ spectra contain a single peak at each
edge, while the Rh4+ spectrum contains a shoulder on the low energy side,
corresponding to the empty $t_{2g}$ state. As expected, this shoulder is
suppressed at the $L_{2}$ edge [50]. To quantify the Rh valence, we obtain the
white line peak position for each compound from a fit to a Lorentzian plus
arctangent lineshape [51, 52]. The results of this fit are summarized in table
2. The reference compounds for Rh3+ show a 0.6 eV (0.5 eV) difference in peak
position at the $L_{3}$ ($L_{2}$) edge. The higher peak position of Rh2O3
compared to RhCl3 can be attributed to the larger covalency of the Rh-O bond
compared to the Rh-Cl bond, giving a more delocalized charge density around
the Rh site in Rh2O3 [53]. We find that the peak position for NiRh2O4 lies
between these two compounds with no significant differences between lineshapes
and peak widths. This confirms that despite the strong hybridization, Rh
maintains the charge distribution of the 3+ oxidation state with no signs of
charge disproportionation, which often occurs in mixed 3d/4d compounds [54,
55, 56].

Figure 4: Calculated RIXS spectrum from the two-site model with $t=30$ meV,
$\Delta=-600$ meV compared with measured data. Magenta ticks indicate
predicted excitation energies.

Having confirmed the strong Ni-Rh hybridization and Rh oxidation state in
NiRh2O4, we assign the RIXS peak at 250 meV to a two-site orbital excitation
from Rh $t_{2g}$ to Ni $t_{2}$, corresponding to the transition $\text{Ni
}3d^{8}+\text{Rh }4d^{6}\rightarrow\text{Ni }3d^{9}+\text{Rh }4d^{5}$. This
assignment is corroborated by our effective two site model discussed below.
Although such metal-metal charge transfer (MMCT) is well known from optical
studies of spinel ferrites [57, 58], there are few reports of these features
seen in RIXS [59, 60, 61].

A full microscopic description of MMCT in NiRh2O4 is considerably complicated
by the 12 fold nearest neighbor Ni-Rh coordination with intermediate Ni-O-Rh
bond angles that prohibit any strict orthogonality constraints on hopping
pathways. In order to capture the essential features of Ni-Rh hybridization
present in the RIXS spectra, we construct a minimal model by adding a single
set of filled Rh $t_{2g}$ states to the single-ion model from the previous
section (see Supplemental Material [37] for details). In this two-site model,
the Ni and Rh sites are each subject to Coulomb, crystal field, and spin-orbit
interactions, with on-site energy difference $\Delta$ and hopping $t$. Since
multiple sites are involved, we now include the monopole part of the Coulomb
interaction at each site $U_{\text{Ni}}=6$ eV and $U_{\text{Rh}}=3$ eV. We
assume an equal hopping from each Ni $t_{2}$ to each Rh $t_{2g}$ and neglect
any hopping from the Ni $e$ levels. The effective crystal field parameters
obtained from a single-ion model must also be adjusted in the presence of
hopping. As shown in Fig. 4, we find good agreement between data and model for
the parameters $\Delta=-600$ meV, $t=30$ meV, $10Dq=-550$ meV, $\lambda=13$
meV, and $\Delta t_{2}=40$ meV, with the same constraint $\Delta e=0.8\Delta
t_{2}$. We also set the Rh tetragonal splitting $\Delta t_{2g}=0$ meV, but the
results are mostly unchanged for nonzero $\Delta t_{2g}$. This corresponds to
scenario 1 discussed above, but similar results are obtained for scenario 2
(see Supplemental Material [37]). We emphasize that the two-site model does
not preserve the symmetry of the Ni site, making it ineffective for comparing
the two scenarios. This is evident from the inset of Fig. 4, where all low-
energy levels have non-negligible dipole character.

Despite the simplicity, the minimal model provides a robust qualitative
description of the data, reproducing the energy of all observed RIXS
excitations over a 2 eV range of energy transfers with effective parameters
describing the approximate energy scales of the Ni-Rh hybridization. A more
detailed approach should incorporate the empty Rh $e_{g}$ levels, which may be
essential to the long range exchange interaction, and might consider a double
cluster model [62, 63], or symmetry-adapted Rh orbitals [64]. In addition,
through comparison of high energy inelastic neutron scattering data that is
sensitive to optical phonons, with the RIXS spectra, we find that lattice
vibrations enter as an essential energy scale coupled to the electronic states
in NiRh2O4. A consideration of electron phonon coupling and vibronic
excitations is thus necessary to capture the broad asymmetric lineshape of the
RIXS spectra, as discussed in the following section.

###  III.3 Electron-phonon coupling

In orbitally degenerate systems, there is a tendency for strong electron-
lattice coupling. In many cases, the degeneracy may be lifted by a static
lattice distortion via the Jahn-Teller (JT) mechanism [40, 65]. Another
possibility is the formation of vibronic modes via the dynamical JT effect,
where orbital degeneracy is broken by coupling to lattice vibrations [66].
When the JT distortion energy is comparable to the spin-orbit coupling, the
system may host a set of spin-orbital-lattice entangled states [5]. In
NiRh2O4, the tetragonal distortion does not fully lift the orbital degeneracy
[14]. The weak tetrahedral crystal field splitting (in comparison to
octahedral) enables the JT energy to be comparable to SOC. Based on these
considerations, we expect lattice dynamics to play a key role in the low-lying
spin-orbital excitations.

Here we consider the effects of electron-phonon coupling on the RIXS spectrum
of NiRh2O4. Although the RIXS spectra may contain contributions from optical
phonons, their precise energies are obscured by the relatively coarse energy
resolution on the scale of the phonon energies, and the coincidence of optical
phonons with low energy electronic excitations. In order to more precisely
quantify optical phonon energies in NiRh2O4, we have re-examined the inelastic
neutron scattering data for energies up to 100 meV, covering the full phonon
bandwidth. The high energy inelastic neutron scattering data is shown in Fig.
5(a).

Figure 5: (a) Inelastic neutron scattering spectra $I(Q,E)$ of NiRh2O4 at
$T=3.6$ K. Collective spin-orbit excitations are visible around 10 meV, and
the highest energy optical phonons appear between 60 and 70 meV. (b) Constant
energy cuts through $I(Q,E)$ integrated over $E\pm 5$ meV, dashed line shows
the momentum dependence of the Ni2+ form-factor $|f(Q)|^{2}$, solid line is a
fit to $C|f(Q)|^{2}+AQ^{2}+B$ as described in the text. (c) Constant momentum-
transfer cut integrated over $Q=6\pm 1$ Å-1, solid line is a fit to three
Lorentzians as described in the text.

Below 20 meV, the previously reported collective spin-orbit excitation is
visible. As discussed above, this feature was not directly resolved in the
RIXS spectra, but accounts for the broadening of the elastic line (inset of
Fig. 1(b) and is captured by the single-ion model. At higher energies, we
observe optical phonons, centered around 65 meV, coincident with the
intra-$t_{2}$ excitation in the RIXS spectrum. The quadratically increasing
intensity with increasing momentum transfer as shown in Fig. 5(b) indicates
that this signal originates primarly from scattering by phonons, but we also
find a component attributable to magnetic scattering. By fitting this cut to
$I(Q)=C|f(Q)|^{2}+AQ^{2}+B$, where $f(Q)$ is the Ni2+ magnetic form-factor,
$A,B$, and $C$ are constants, we obtain a good description of the data with
the parameters $A=0.04$, $B=0.2$, $C=0.78$. Optical phonon energies were
extracted directly from the constant momentum transfer cut in Fig. 5(c). We
found that including three Lorentzian functions centered at 57.2(6), 65.3(5),
and 72(1) meV and with energy linewidths of $\Gamma=$ 7.6(1.7), 7.5(2.5), and
8.5(2.0) meV respectively was necessary to adequately describe the data. The
phonon linewidths are significantly broadened over the instrumental resolution
of $\sim 4.3$-meV Gaussian FWHM at 65 meV. These high energy phonons originate
from vibrations of Ni coordinating O tetrahedra, and we expect six distinct
modes in this energy range for a cubic cell that are further split in the
tetragonal phase [67, 40, 68]. Although the broad phonon lineshapes we
observed may be accounted for by unresolved phonon mode splittings, it may
also indicate phonon damping caused by coupling to other electronic
excitations. Indeed, the coinciding energy of these phonons and the
intra-$t_{2}$ crystal field excitation measured by RIXS suggests the
possibility of a hybridized vibrational-electronic or “vibronic” excitation in
NiRh2O4.

Such electron phonon coupling occurs because the charge distribution of the
excited state on Ni repels the surrounding oxygen ions. For the intra $t_{2}$
excitation at 65 meV, the extra electron in the $d_{xy}$ orbital is partially
screened by the resulting $d_{xz}/d_{yz}$ hole. However, the $e\rightarrow
t_{2}$ excitation at 500 meV similarly leaves an extra charge in the $t_{2}$
levels, with a hole in the $e$ levels. In this case, we expect the screening
to be less effective, and thus the $e\rightarrow t_{2}$ excitation should
couple even more strongly to phonons. We can use these considerations to model
the lineshape of the 500 meV RIXS peak as a vibronic excitation. Although the
single-ion model indicates at least two states comprise this peak and the
neutron data shows at least three phonon modes may be involved in the
coupling, we will consider a tractable model that includes only a single
electronic excitation coupled to a single phonon as such a model is sufficient
to capture the essential features of our data. Within this simplified model,
we treat the main peak at energy $E_{0}\\!=\\!496$ meV as a bare $d$-$d$
excitation, or zero-phonon line. This bare dd-excitation is dressed by
additional phonon sidebands corresponding to a $d\text{-}d$ excitation plus
$n$ phonons of energy $E_{\text{ph}}$. We assume a Lorentzian lineshape of
FWHM width $\Gamma$ for each peak separated by energy $E_{\text{ph}}$ with
relative intensities given by a Poisson distribution [69]

$I_{n}=e^{-g}\frac{g^{n}}{n!},$ (1)

where $g$ is the dimensionless electron-phonon coupling, which can be
interpreted as the mean number of phonons emitted by the excitation.

Figure 6: Poisson model for phonon sidebands fit to the 500 meV RIXS peak,
including Gaussian fits to the overlapping neighboring peaks. The phonon
sideband cross-section was modeled using equation 1 and Huang-Rhys parameter
$g=0.46(3)$, $\omega_{0}=71(2)$ meV, $E_{0}=497(2)$ meV, and $\Gamma=60(2)$
meV.

Fig. 6 shows a fit of this minimal model to the 500 meV RIXS peak for
$E_{0}=496(1)$ meV, $E_{\text{ph}}=68(2)$ meV, $\Gamma=70(2)$ meV, and
$g=0.45(3)$ providing an excellent description of the data. The large value of
$\Gamma$ compared to the phonon energies suggests that there are many
overlapping excitations within each phonon sideband. This is consistent with
the single-ion model for scenario 1 that gives two crystal field excitations,
at 492 and 517 meV. Both of these may couple to optical phonons with energies
between 57 and 72 meV. A more complete model for electron phonon coupling in
NiRh2O4 would consider the separate $d$-$d$ excitations and their coupling to
multiple phonons, but the resolution of our measurement is not sufficient to
constrain such a model.

##  IV Discussion and Conclusion

We have characterized the site-specific local electronic structure of NiRh2O4
using resonant inelastic x-ray scattering and x-ray absorption spectroscopy.
We have compared two possible scenarios for the tetragonal splitting within a
single ion model and showed that $\Delta t_{2}>0$ is more likely than $\Delta
t_{2}<0$, and estimated that $\Delta t_{2}=70$ meV and $\lambda=13$ meV. These
parameters are the most relevant to modeling the magnetism in NiRh2O4 [32,
35]. The crystal field splittings are in agreement with DFT calculations that
determined $\Delta t_{2}=100$ meV and $\lambda=10$ meV from NMTO downfolding.
The single ion model also required a 50% reduction of the $F_{dd}$ Slater
parameters, suggesting a significant degree of covalency in NiRh2O4.

The O $K$ edge XAS data suggests a significant degree of hybridization between
O $p$ states and empty metal $d$ states, which can only be explained by a
metal-metal charge transfer between adjacent Ni and Rh sites. The Rh $L$ edge
XAS confirmed that the Rh ions maintain the nominal $3+$ oxidation state,
despite the strong hybridization and tendency toward charge disproportionation
in related systems. By extending the RIXS model to include metal-metal charge
transfer between Ni and Rh sites as parameterized by an effective hopping
$t\\!=\\!30$ meV, we captured Rh-Ni two-site excitations observed at 250 and
750 meV in the RIXS spectrum. Such an explicit demonstration of the failure of
a single ion model and requirement for metal-metal charge transfer in NiRh2O4
highlights the importance of metal-metal hybridization in mixed 3d-4d/5d
compounds in general. This hybridization can affect both the magnetic degrees
of freedom and exchange interactions so should be an essential consideration
in the design of novel magnetic states in materials such as spinels, double
perovskites [70, 71], or $A_{2}$Mo3O8 compounds [28, 72].

A detailed analysis of the RIXS lineshape also revealed that lattice
vibrations influence the magnetism in NiRh2O4 through strong electron-phonon
coupling. Inelastic neutron scattering reveals multiple optical phonons
overlapping in energy with the intra-$t_{2}$ excitation, suggesting a
hybridized orbital-lattice excitation between 60 and 70 meV. This effect was
observed in the RIXS spectra as a phonon-dressed crystal field excitation at
higher energy. Our results provide quantitative constraints on the key
parameters for modeling the single-ion ground state and low-lying excitations,
as well the long-range superexchange mechanism in NiRh2O4. We also demonstrate
the importance of additional degrees of freedom, namely covalency and phonons,
which can alter magnetic ground states [73]. Our results also demonstrate the
use of RIXS for probing spin-orbit entangled states in $3d$ transition metal
compounds [74] and for probing hybridized states in mixed 3d-4d/5d compounds
[75].

In addition to the predicted phenomena associated with spin-1 frustrated
diamond lattice antiferromagnets, we propose that NiRh2O4 may host a variety
of novel magneto-elastic and magneto-optical effects due to its rich spectrum
of low-energy spin-orbital-lattice excitations [40, 65, 76, 77]. Future
studies on NiRh2O4 could use pressure, magnetic fields or chemical
substitution to explore the predicted phase diagram and observe quantum
critical phenomena [31, 30, 32]. These studies would greatly benefit from the
synthesis of single crystals or thin films. Should such crystals become
available, another potentially fruitful route would be ultrafast optics to
study the spin-orbital and crystal field excitations with the lattice out-of-
equilibrium [78, 79, 80]. By resonantly exciting the lattice, it may be
possible to stabilize an excitonic condensate of $J=1$ moments [32].

## Acknowledgements

We thank Arun Paramekanti, Pat Clancy, and Hlynur Gretarsson for useful
discussions. We are grateful to Matthew Stone for assistance in collecting
data at ORNL. We thank Fazel Tafti and Faranak Bahrami for providing the
Ag3LiRh2O6 reference sample for XAS. Work at Brown University was supported by
the U.S. Department of Energy, Office of Science, Office of Basic Energy
Sciences, under Award Number DE-SC002165. TMM and JC acknowledge support from
the Institute for Quantum Matter, an Energy Frontier Research Center funded by
the U.S. Department of Energy, Office of Science, Office of Basic Energy
Sciences, under Award DE-SC0019331.The work of L.G. and M.M. at Georgia Tech
(inelastic neutron scattering measurements) was supported by NSF-DMR-1750186.
This research used beamline 2-ID of the National Synchrotron Light Source II,
a U.S. Department of Energy (DOE) Office of Science User Facility operated for
the DOE Office of Science by Brookhaven National Laboratory under Contract No.
DE-SC0012704. Use of the Advanced Photon Source at Argonne National Laboratory
was supported by the U. S. Department of Energy, Office of Science, Office of
Basic Energy Sciences, under Contract No. DE-AC02-06CH11357.A portion of this
research used resources at the Spallation Neutron Source, a DOE Office of
Science User Facility operated by the Oak Ridge National Laboratory.

## References

  * Lacroix  _et al._ [2011] C. Lacroix, P. Mendels, and F. Mila,  _Introduction to Frustrated Magnetism: Materials, Experiments, Theory_, Vol. 164 (Springer Science & Business Media, 2011).
  * Khomskii [2014] D. I. Khomskii,  _Transition Metal Compounds_ (Cambridge University Press, 2014).
  * Basov  _et al._ [2017] D. Basov, R. Averitt, and D. Hsieh, Nature materials 16, 1077 (2017).
  * Takayama  _et al._ [2021] T. Takayama, J. Chaloupka, A. Smerald, G. Khaliullin, and H. Takagi, Journal of the Physical Society of Japan 90, 062001 (2021).
  * Streltsov and Khomskii [2020] S. V. Streltsov and D. I. Khomskii, Phys. Rev. X 10, 031043 (2020).
  * Tsurkan  _et al._ [2021] V. Tsurkan, H.-A. Krug von Nidda, J. Deisenhofer, P. Lunkenheimer, and A. Loidl, Physics Reports 926, 1 (2021).
  * Roth [1964] W. Roth, Journal of Physics and Chemistry of Solids 25, 1 (1964).
  * Tristan  _et al._ [2005] N. Tristan, J. Hemberger, A. Krimmel, H.-A. Krug von Nidda, V. Tsurkan, and A. Loidl, Phys. Rev. B 72, 174404 (2005).
  * Krimmel  _et al._ [2009] A. Krimmel, H. Mutka, M. M. Koza, V. Tsurkan, and A. Loidl, Phys. Rev. B 79, 134406 (2009).
  * Fritsch  _et al._ [2004] V. Fritsch, J. Hemberger, N. Büttgen, E.-W. Scheidt, H.-A. Krug von Nidda, A. Loidl, and V. Tsurkan, Phys. Rev. Lett. 92, 116401 (2004).
  * Gao  _et al._ [2017] S. Gao, O. Zaharko, V. Tsurkan, Y. Su, J. S. White, G. S. Tucker, B. Roessli, F. Bourdarot, R. Sibille, D. Chernyshov, _et al._ , Nature physics 13, 157 (2017).
  * Guratinder  _et al._ [2022] K. Guratinder, V. Tsurkan, L. Prodan, L. Keller, J. P. Embs, F. Juranyi, M. Medarde, C. Rüegg, and O. Zaharko, Phys. Rev. B 105, 174423 (2022).
  * Ge  _et al._ [2017] L. Ge, J. Flynn, J. A. M. Paddison, M. B. Stone, S. Calder, M. A. Subramanian, A. P. Ramirez, and M. Mourigal, Phys. Rev. B 96, 064413 (2017).
  * Chamorro  _et al._ [2018] J. R. Chamorro, L. Ge, J. Flynn, M. A. Subramanian, M. Mourigal, and T. M. McQueen, Phys. Rev. Materials 2, 034404 (2018).
  * Bordelon  _et al._ [2021a] M. M. Bordelon, C. Liu, L. Posthuma, E. Kenney, M. J. Graf, N. P. Butch, A. Banerjee, S. Calder, L. Balents, and S. D. Wilson, Phys. Rev. B 103, 014420 (2021a).
  * Bordelon  _et al._ [2021b] M. M. Bordelon, J. D. Bocarsly, L. Posthuma, A. Banerjee, Q. Zhang, and S. D. Wilson, Phys. Rev. B 103, 024430 (2021b).
  * Chen  _et al._ [2009] G. Chen, L. Balents, and A. P. Schnyder, Phys. Rev. Lett. 102, 096406 (2009).
  * Savary  _et al._ [2011] L. Savary, E. Gull, S. Trebst, J. Alicea, D. Bergman, and L. Balents, Phys. Rev. B 84, 064438 (2011).
  * Nair  _et al._ [2014] H. S. Nair, Z. Fu, J. Voigt, Y. Su, and T. Brückel, Phys. Rev. B 89, 174431 (2014).
  * MacDougall  _et al._ [2016] G. J. MacDougall, A. A. Aczel, Y. Su, W. Schweika, E. Faulhaber, A. Schneidewind, A. D. Christianson, J. L. Zarestky, H. D. Zhou, D. Mandrus, and S. E. Nagler, Phys. Rev. B 94, 184422 (2016).
  * Plumb  _et al._ [2016] K. W. Plumb, J. R. Morey, J. A. Rodriguez-Rivera, H. Wu, A. A. Podlesnyak, T. M. McQueen, and C. L. Broholm, Phys. Rev. X 6, 041055 (2016).
  * Biffin  _et al._ [2017] A. Biffin, C. Rüegg, J. Embs, T. Guidi, D. Cheptiakov, A. Loidl, V. Tsurkan, and R. Coldea, Phys. Rev. Lett. 118, 067205 (2017).
  * Tsurkan  _et al._ [2017] V. Tsurkan, L. Prodan, V. Felea, I. Filippova, V. Kravtsov, A. Günther, S. Widmann, H.-A. Krug von Nidda, J. Deisenhofer, and A. Loidl, Phys. Rev. B 96, 054417 (2017).
  * Naka  _et al._ [2020] T. Naka, J. Valenta, J. Kaštil, M. Míšek, J. Prchal, V. Sechovský, H. Abe, T. Nakane, M. Nakayama, and T. Uchikoshi, Materials Research Express 7, 056105 (2020).
  * Cho  _et al._ [2021] H. Cho, C. H. Kim, Y. Lee, K. Komatsu, B.-G. Cho, D.-Y. Cho, T. Kim, C. Kim, Y. Kim, T. Y. Koo, Y. Noda, H. Kagi, D. I. Khomskii, D. Seoung, and J.-G. Park, Phys. Rev. B 103, L081101 (2021).
  * Li and Chen [2018] F.-Y. Li and G. Chen, Phys. Rev. B 98, 045109 (2018).
  * Liu  _et al._ [2020] J. Q. Liu, F.-Y. Li, G. Chen, and Z. Wang, Phys. Rev. Research 2, 033260 (2020).
  * Morey  _et al._ [2019] J. R. Morey, A. Scheie, J. P. Sheckelton, C. M. Brown, and T. M. McQueen, Phys. Rev. Materials 3, 014410 (2019).
  * Wang  _et al._ [2015] C. Wang, A. Nahum, and T. Senthil, Phys. Rev. B 91, 195131 (2015).
  * Buessen  _et al._ [2018] F. L. Buessen, M. Hering, J. Reuther, and S. Trebst, Phys. Rev. Lett. 120, 057201 (2018).
  * Chen [2017] G. Chen, Phys. Rev. B 96, 020412(R) (2017).
  * Li and Chen [2019] F.-Y. Li and G. Chen, Phys. Rev. B 100, 045103 (2019).
  * Horiuti and Miyahara [1964] S. Horiuti and S. Miyahara, Journal of the Physical Society of Japan 19, 423 (1964).
  * Blasse and Schipper [1963] G. Blasse and D. Schipper, Physics Letters 5, 300 (1963).
  * Das  _et al._ [2019] S. Das, D. Nafday, T. Saha-Dasgupta, and A. Paramekanti, Phys. Rev. B 100, 140408(R) (2019).
  * Wang  _et al._ [2019] Y. Wang, G. Fabbris, M. Dean, and G. Kotliar, Computer Physics Communications 243, 151 (2019).
  * Note [1] See Supplemental Material.
  * Cowan [1981] R. D. Cowan,  _The Theory of Atomic Structure and Spectra_ (Univ of California Press, 1981).
  * Gütlich  _et al._ [1984] P. Gütlich, K. M. Hasselbach, H. Rummel, and H. Spiering, The Journal of Chemical Physics 81, 1396 (1984).
  * Kocsis  _et al._ [2013] V. Kocsis, S. Bordács, D. Varjas, K. Penc, A. Abouelsayed, C. A. Kuntscher, K. Ohgushi, Y. Tokura, and I. Kézsmárki, Phys. Rev. B 87, 064416 (2013).
  * de Groot [1994] F. de Groot, Journal of Electron Spectroscopy and Related Phenomena 67, 529 (1994).
  * Kurmaev  _et al._ [2008] E. Z. Kurmaev, R. G. Wilks, A. Moewes, L. D. Finkelstein, S. N. Shamin, and J. Kuneš, Phys. Rev. B 77, 165127 (2008).
  * Wahila  _et al._ [2019] M. J. Wahila, Z. W. Lebens-Higgins, A. J. Jackson, D. O. Scanlon, T.-L. Lee, J. Zhang, K. H. L. Zhang, and L. F. J. Piper, Phys. Rev. B 100, 085126 (2019).
  * Frati  _et al._ [2020] F. Frati, M. O. Hunault, and F. M. de Groot, Chemical reviews 120, 4056 (2020).
  * Lenglet [2000] M. Lenglet, Materials Research Bulletin 35, 531 (2000).
  * Mayer and Angelov [1981] I. Mayer and S. Angelov, physica status solidi (b) 106, 467 (1981).
  * Zhu  _et al._ [2015] M. Zhu, D. Do, C. R. Dela Cruz, Z. Dun, J.-G. Cheng, H. Goto, Y. Uwatoko, T. Zou, H. D. Zhou, S. D. Mahanti, and X. Ke, Phys. Rev. B 92, 094419 (2015).
  * Katukuri  _et al._ [2020] V. M. Katukuri, P. Babkevich, O. Mustonen, H. C. Walker, B. Fåk, S. Vasala, M. Karppinen, H. M. Rønnow, and O. V. Yazyev, Phys. Rev. Lett. 124, 077202 (2020).
  * Bahrami  _et al._ [2022] F. Bahrami, X. Hu, Y. Du, O. I. Lebedev, C. Wang, H. Luetkens, G. Fabbris, M. J. Graf, D. Haskel, Y. Ran, and F. Tafti, Science Advances 8, eabl5671 (2022).
  * Burnus  _et al._ [2008] T. Burnus, Z. Hu, H. Wu, J. C. Cezar, S. Niitaka, H. Takagi, C. F. Chang, N. B. Brookes, H.-J. Lin, L. Y. Jang, A. Tanaka, K. S. Liang, C. T. Chen, and L. H. Tjeng, Phys. Rev. B 77, 205111 (2008).
  * Clancy  _et al._ [2014] J. P. Clancy, A. Lupascu, H. Gretarsson, Z. Islam, Y. F. Hu, D. Casa, C. S. Nelson, S. C. LaMarra, G. Cao, and Y.-J. Kim, Phys. Rev. B 89, 054409 (2014).
  * Chikara  _et al._ [2017] S. Chikara, G. Fabbris, J. Terzic, G. Cao, D. Khomskii, and D. Haskel, Phys. Rev. B 95, 060407(R) (2017).
  * Wu and Ellis [1995] Y. Wu and D. E. Ellis, Journal of Physics: Condensed Matter 7, 3973 (1995).
  * Meyers  _et al._ [2014] D. Meyers, S. Middey, J.-G. Cheng, S. Mukherjee, B. Gray, Y. Cao, J.-S. Zhou, J. Goodenough, Y. Choi, D. Haskel, _et al._ , Nature communications 5, 1 (2014).
  * Kobayashi  _et al._ [2022] M. Kobayashi, M. Seki, M. Suzuki, R. Okano, M. Kitamura, K. Horiba, H. Kumigashira, A. Fujimori, M. Tanaka, and H. Tabata, Phys. Rev. B 105, 205103 (2022).
  * Takubo  _et al._ [2005] K. Takubo, S. Hirata, J.-Y. Son, J. W. Quilty, T. Mizokawa, N. Matsumoto, and S. Nagata, Phys. Rev. Lett. 95, 246401 (2005).
  * Fontijn  _et al._ [1999] W. F. J. Fontijn, P. J. van der Zaag, L. F. Feiner, R. Metselaar, and M. A. C. Devillers, Journal of Applied Physics 85, 5100 (1999).
  * Kim  _et al._ [2002] K. J. Kim, H. S. Lee, M. H. Lee, and S. H. Lee, Journal of Applied Physics 91, 9974 (2002).
  * Agui  _et al._ [2009] A. Agui, T. Uozumi, M. Mizumaki, and T. Käämbre, Phys. Rev. B 79, 092402 (2009).
  * Chiuzbăian  _et al._ [2017] S. G. Chiuzbăian, S. Brignolo, C. F. Hague, R. Delaunay, M. Guarise, A. Nicolaou, Z. Yang, H. Zhou, and J.-M. Mariot, The Journal of Physical Chemistry C 121, 22369 (2017).
  * Liu  _et al._ [2018] B. Liu, E. N. Glass, R.-P. Wang, Y.-T. Cui, Y. Harada, D.-J. Huang, S. Schuppler, C. L. Hill, and F. M. F. de Groot, Phys. Chem. Chem. Phys. 20, 4554 (2018).
  * Green  _et al._ [2016] R. J. Green, M. W. Haverkort, and G. A. Sawatzky, Phys. Rev. B 94, 195127 (2016).
  * Guedes  _et al._ [2019] E. B. Guedes, F. Abud, H. P. Martins, M. Abbate, R. F. Jardim, and R. J. O. Mossanek, Phys. Rev. B 100, 075132 (2019).
  * Haverkort  _et al._ [2012] M. W. Haverkort, M. Zwierzycki, and O. K. Andersen, Phys. Rev. B 85, 165113 (2012).
  * Kocsis  _et al._ [2018] V. Kocsis, S. Bordács, J. Deisenhofer, L. F. Kiss, K. Ohgushi, Y. Kaneko, Y. Tokura, and I. Kézsmárki, Phys. Rev. B 97, 125140 (2018).
  * Krimmel  _et al._ [2005] A. Krimmel, M. Mücksch, V. Tsurkan, M. M. Koza, H. Mutka, and A. Loidl, Phys. Rev. Lett. 94, 237402 (2005).
  * Ptak  _et al._ [2013] M. Ptak, M. Maczka, A. Gągor, A. Pikul, L. Macalik, and J. Hanuza, Journal of Solid State Chemistry 201, 270 (2013).
  * Wang  _et al._ [2003] Z. Wang, S. Saxena, P. Lazor, and H. O’Neill, Journal of Physics and Chemistry of Solids 64, 425 (2003).
  * Geondzhian  _et al._ [2020] A. Geondzhian, A. Sambri, G. M. De Luca, R. Di Capua, E. Di Gennaro, D. Betto, M. Rossi, Y. Y. Peng, R. Fumagalli, N. B. Brookes, L. Braicovich, K. Gilmore, G. Ghiringhelli, and M. Salluzzo, Phys. Rev. Lett. 125, 126401 (2020).
  * Saha-Dasgupta [2020] T. Saha-Dasgupta, Materials Research Express 7, 014003 (2020).
  * Lee  _et al._ [2018] M.-C. Lee, S. Lee, C. J. Won, K. D. Lee, N. Hur, J.-L. Chen, D.-Y. Cho, and T. W. Noh, Phys. Rev. B 97, 125123 (2018).
  * Park  _et al._ [2021] K. Park, G. L. Pascut, G. Khanal, M. O. Yokosuk, X. Xu, B. Gao, M. J. Gutmann, A. P. Litvinchuk, V. Kiryukhin, S. W. Cheong, D. Vanderbilt, K. Haule, and J. L. Musfeldt, Phys. Rev. B 104, 195143 (2021).
  * Xu  _et al._ [2016] L. Xu, N. A. Bogdanov, A. Princep, P. Fulde, J. Van Den Brink, and L. Hozoi, npj Quantum Materials 1, 1 (2016).
  * Huang  _et al._ [2022] H. Huang, A. Singh, C. Wu, J. Xie, J. Okamoto, A. Belik, E. Kurmaev, A. Fujimori, C. Chen, S. Streltsov, _et al._ , npj Quantum Materials 7, 1 (2022).
  * Jin  _et al._ [2022] W. Jin, S. H. Chun, J. Kim, D. Casa, J. P. C. Ruff, C. J. Won, K. D. Lee, N. Hur, and Y.-J. Kim, Phys. Rev. B 105, 054419 (2022).
  * Smerald and Jackeli [2019] A. Smerald and G. Jackeli, Phys. Rev. Lett. 122, 227202 (2019).
  * Amelin  _et al._ [2020] K. Amelin, Y. Alexanian, U. Nagel, T. Rõõm, J. Robert, J. Debray, V. Simonet, C. Decorse, Z. Wang, R. Ballou, E. Constable, and S. de Brion, Phys. Rev. B 102, 134428 (2020).
  * Marciniak  _et al._ [2021] A. Marciniak, S. Marcantoni, F. Giusti, F. Glerean, G. Sparapassi, T. Nova, A. Cartella, S. Latini, F. Valiera, A. Rubio, _et al._ , Nature Physics 17, 368 (2021).
  * Afanasiev  _et al._ [2021] D. Afanasiev, J. Hortensius, B. Ivanov, A. Sasani, E. Bousquet, Y. Blanter, R. Mikhaylovskiy, A. Kimel, and A. Caviglia, Nature materials 20, 607 (2021).
  * Ergeçen  _et al._ [2022] E. Ergeçen, B. Ilyas, D. Mao, H. C. Po, M. B. Yilmaz, J. Kim, J.-G. Park, T. Senthil, and N. Gedik, Nature Communications 13, 1 (2022).