Spin-orbit entangled j=1/2 moments in Ba$_2$CeIrO$_6$ -- a frustrated fcc quantum magnet

TL;DR

Ba$_2$CeIrO$_6$ is an almost ideal j=1/2 Mott insulator with frustrated fcc lattice magnetism, exhibiting antiferromagnetic order driven by Kitaev interactions and magneto-elastic effects, revealing complex spin-orbit physics.

Contribution

This work identifies Ba$_2$CeIrO$_6$ as a near-ideal model for j=1/2 moments and uncovers the role of Kitaev exchange and lattice distortions in its magnetic ordering.

Findings

Less than 1% deviation from ideal j=1/2 state.

Antiferromagnetic order at 14K despite high Curie-Weiss temperature.

Enhanced geometric frustration due to next-nearest neighbor exchange.

Abstract

We establish the double perovskite Ba$_2$CeIrO$_6$ as a nearly ideal model system for j=1/2 moments, with resonant inelastic x-ray scattering indicating a deviation of less than 1% from the ideally cubic j=1/2 state. The local j=1/2 moments form an fcc lattice and are found to order antiferromagnetically at $T_N$=14K, more than an order of magnitude below the Curie-Weiss temperature. Model calculations show that the geometric frustration of the fcc Heisenberg antiferromagnet is further enhanced by a next-nearest neighbor exchange, indicated by ab initio theory. Magnetic order is driven by a bond-directional Kitaev exchange and by local distortions via a strong magneto-elastic effect - both effects are typically not expected for j=1/2 compounds making Ba2CeIrO6 a riveting example for the rich physics of spin-orbit entangled Mott insulators.

Tables

Table 1. Table 1: Results of structure refinements with single-crystal X-ray data obtained at room temperature and 100 K. In space group F ​ m ​ 3 ¯ ​ m 𝐹 𝑚 ¯ 3 𝑚 Fm\bar{3}m Ba occupies an 8 ​ c 8 𝑐 8c site at (1/4,1/4,1/4), Ce a 4 ​ a 4 𝑎 4a site at (0,0,0), Ir a 4 ​ b 4 𝑏 4b site at (0,0,1/2), and O a 24 ​ e 24 𝑒 24e site at ( x 𝑥 x ,0,0). Thermal parameters are given in 10 − 5 superscript 10 5 10^{-5} Å 2 , and only the O parameters are anisotropic. The weighted R 𝑅 R values for the structure factor amount to 4.22 % and 3.96 % at 300 K and 100 K, respectively, and refining occupations reduces them to 2.57 % and 3.11 %, respectively. The occupation of Ba and O was fixed. For the occupation of Ce and Ir we find 99.6 99.6 99.6 % ( 99.1 99.1 99.1 %) and 94.6 94.6 94.6 % ( 95.0 95.0 95.0 %) at 300 K (100 K), respectively, with an uncertainty of roughly 1%.

$T$(K)	$U$(Ba)	$U$(Ce)	$U$(Ir)	$x$(O)	$U_{11}$(O)	$U_{22}$(O)
300	1299(8)	740(8)	475(6)	0.2592(3)	1230(140)	2390(110)
100	1087(13)	738(18)	585(12)	0.2590(7)	1100(300)	2200(200)

Equations

$$t_p=\frac{r_A+r_O}{\sqrt{2}(r_B+r_O)}$$

$$\chi (T)=N_A\frac{{\mu }_{\mathrm{eff}}^2}{3k_B(T-{\Theta }_{\mathrm{CW}})}+{\chi }_0,$$

$$H_{\mathrm{single}}=\lambda \vec{S}\cdot \vec{L}+{\Delta }_{\mathrm{CF}}{L}_z^2,$$

$$|0⟩=0.991|\frac{1}{2},\frac{1}{2}⟩-0.130|\frac{3}{2},\frac{1}{2}⟩$$

$$\mathcal{H}=J_1\sum _{⟨i,j⟩}{\vec{\mathcal{S}}}_i\cdot {\vec{\mathcal{S}}}_j+K\sum _{⟨i,j{⟩}_{\gamma }}{\mathcal{S}}_i^{\gamma }{\mathcal{S}}_j^{\gamma }+J_2\sum _{⟨⟨i,j⟩⟩}{\vec{\mathcal{S}}}_i\cdot {\vec{\mathcal{S}}}_j,$$

$$J_{1}\to\left\{\begin{array}[]{l}J_{1}+\delta J(\Delta_{\rm CF}),\hskip 23.9002ptxy\,\,{\rm plane},\\ J_{1}-\delta J(\Delta_{\rm CF}),\hskip 7.11317ptyz,\,xz\,\,{\rm planes},\end{array}\right.$$

$$\delta J(\Delta_{\rm CF})=\frac{\Delta_{\rm CF}}{\Delta^{\rm max}_{\rm CF}}\big{(}J(\Delta^{\rm max}_{\rm CF})-J(\Delta_{\rm CF}=0)\big{)}\,.$$

$$H_{\mathrm{site}}$$

$$|\frac{1}{2},\mathrm{up}⟩=\frac{1}{\sqrt{3}}\left(|yz,\downarrow ⟩+i|xz,\downarrow ⟩+|xy,\uparrow ⟩\right).$$

Full text

I Introduction

Spin-orbit entangled Mott insulators stand out in the growing family of quantum materials with strong spin-orbit coupling for their correlation-driven phenomena Witczak-Krempa et al. (2014). Of particular interest are materials with partially filled $4d$ and $5d$ orbitals, such as the iridates, in which the formation of local $j$ = $1/2$ moments is an iridescent source of rich physics Rau et al. (2016). The spin-orbit entangled wavefunction of these Kramers doublets gives rise to fundamentally different types of exchange interactions depending on the geometric arrangement of the elementary octahedral IrO6 building blocks Khaliullin (2005); Chen and Balents (2008); Jackeli and Khaliullin (2009). Corner-sharing octahedra yield isotropic Heisenberg exchange, which has been explored as a potential source of spin-orbit assisted superconductivity Kim et al. (2014a); de la Torre et al. (2015); Cao et al. (2016); Yan et al. (2015); Kim et al. (2016) in the context of Sr2IrO4 Kim et al. (2008, 2009), an isostructural analogue of the high-Tc parent compound La2CuO4. Edge-sharing octahedra, in contrast, give rise to Kitaev-type bond-directional exchange, which has initiated an intense search for spin-orbit driven frustrated quantum magnetism in so-called Kitaev materials Trebst (2017) such as the honeycomb iridates Na2IrO3, $\alpha$-Li2IrO3, and H3LiIr2O6 Singh and Gegenwart (2010); Singh et al. (2012); Kitagawa et al. (2018) and the related $\alpha$-RuCl3 Plumb et al. (2014). Possibly the most spectacular experimental result in this realm is the recent claim of a quantized thermal Hall effect in $\alpha$-RuCl3 Kasahara et al. (2018), a direct signature of the long sought-after Kitaev spin liquid Kitaev (2006).

In this paper, we first demonstrate experimentally that the double perovskite Ba2CeIrO6 is a nearly ideal realization of a $j$ = $1/2$ Mott insulator, forming a model system for frustrated quantum magnetism on the fcc lattice. Our x-ray diffraction results show a global cubic $Fm\bar{3}m$ structure, while resonant inelastic x-ray scattering (RIXS) reveals a small non-cubic distortion resulting in a ground-state wavefunction which overlaps by more than 99 % with the ideal cubic $j$ = 1/2 state. The magnetic susceptibility shows an antiferromagnetic ordering temperature $T_{N}$ = 14 K which is suppressed by more than an order of magnitude compared to the Curie-Weiss temperature $|\Theta_{\rm CW}|$, resulting in a large frustration parameter $f$ = $|\Theta_{\rm CW}|/T_{N}\gtrsim 13$. Employing a combination of density functional theory and microscopic model simulations, we address the minimal model for Ba2CeIrO6 and its phase diagram. The system shows a particularly high degree of frustration, since the geometric frustration of antiferromagnetic nearest-neighbor Heisenberg exchange on the fcc lattice is augmented by next-nearest-neighbor Heisenberg coupling, yielding a wide window of a quantum spin liquid ground state. However, an antiferromagnetic Kitaev-type bond-directional exchange is found to counteract this geometric frustration and turns out to be instrumental in stabilizing long-range magnetic order – in contrast to the common wisdom that Kitaev interactions in $j$ = $1/2$ compounds enhance frustration and induce spin liquid physics.

It is also common wisdom that the $j$ = $1/2$ wavefunction does not show orbital degeneracy and hence is not Jahn-Teller active. Commonly, this is interpreted as a protection of $j$ = $1/2$ physics against lattice distortions, and small deviations from cubic symmetry with the concomitant change of the wavefunction are typically neglected. We challenge this point of view and provide theoretical evidence for a strong magneto-elastic effect. Our theoretical analysis shows that even small deviations from the $j$ = $1/2$ wavefunction, associated with small lattice distortions, yield a massive bond-dependent variation of the nearest-neighbor exchange constants, as illustrated in Fig. 1c), effectively lifting the strong magnetic frustration. This dramatic magneto-elastic coupling is of general importance in the quest for exotic spin liquids based on $j$ = $1/2$ compounds.

II Synthesis and structure

Single crystals of Ba2CeIrO6 of about 1 mm3 size were grown by melt solution growth (see Appendix A). X-ray diffraction shows a well ordered double perovskite with Ce-Ir order as illustrated in Fig. 1a. The cation order can be explained by the notably different bond lengths of 2.20 Å for Ce-O and 2.04 Å for Ir-O. For $5d^{5}$ Ir4+, the formation of ideal $j$ = $1/2$ moments requires a cubic crystal field. Thus far, deviations from cubic symmetry were reported for all $5d^{5}$ iridate compounds Rau et al. (2016); Trebst (2017), as discussed in more detail in the section on RIXS below.

For Ba2CeIrO6, our powder diffraction peaks – measured using a Stoe Stadi MP 198 powder diffractometer – are very well described in the cubic space group $Fm\bar{3}m$ with a lattice constant of 8.47Å at 300 K. However, we find a clear broadening of Bragg peaks in particular for large diffraction angles $2\theta>80^{\circ}$. Such broadened Bragg peaks may explain a previous claim of tiny ($<\!0.2\,\%$) monoclinic distortions of the metric in polycrystalline Ba2CeIrO6 Wakeshima et al. (2000). Note that the issue of cubic or non-cubic symmetry is often discussed controversially in double perovskites, for instance for the closely related Ba2PrIrO6 Wakeshima et al. (2000); Kockelmann et al. (2006).

To resolve this issue, we collected single-crystal x-ray diffraction data. Our results strongly support a cubic structure of Ba2CeIrO6. We employed a Bruker X8 Apex diffractometer, a sample with octahedral shape ($\{111\}$ faces), and a distance to the center of 12.5 $\mu$m. At room temperature (100 K) 32921 (13274) Bragg reflection intensities were recorded, yielding 198 (205) independent reflections in space group $Fm\bar{3}m$. The single-crystal data do not yield any evidence for significant superstructure reflections with respect to $Fm\bar{3}m$, neither at room temperature nor at 100 K (see Appendix A). From this and the description of the powder diffraction pattern with the cubic lattice we must conclude that the average structure of Ba2CeIrO6 is cubic.

However, the atomic displacement parameters shown in table I provide evidence for local distortions since they are (i) larger than expected for a purely dynamical displacement, (ii) very similar at 300 K and 100 K, and (iii) similar for the heavy Ba ions and the lighter O ions. The large values observed for O perpendicular to its bond at room temperature reflect the general instability of a perovskite against tilting. But the small difference in the room-temperature and 100 K displacement values in general indicates some local distortions. Moreover, a normal dynamical effect cannot explain the fact that the atomic displacement parameter of the heavy Ba is of the same magnitude (a root mean square displacement of the order of 0.1 Å) as the one of the much lighter O. Fits of the data in space groups with the same translation lattice but reduced symmetry do not yield significant improvement. However, a split model in which the Ba ions are statistically distributed over sites slightly displaced by $\delta_{\rm Ba}$ against the cubic (0.25,0.25,0.25) position results in $\delta_{\rm Ba}$ = 0.13(1) Å and 0.14(2) Å at 300 K and 100 K, respectively. The statistical character may be related to the existence of about 5 % of vacancies on the Ir sites.

One example for structural distortions that were first sensed by enlarged atomic displacement parameters is KH2PO4, a prototype ferroelectric material that exhibits highly enlarged atomic displacement parameters above its ferroelectric transition of order-disorder character Nelmes . Another example is La1.85Sr0.15CuO4, in which enhanced atomic displacement factors are observed in samples which do not show long-range tilt order Braden et al. (2001). In Ba2CeIrO6, the presence of local distortions from cubic symmetry is supported by our RIXS data, see below. This can be reconciled with the observation of a global cubic structure in x-ray diffraction by assuming a negligible correlation length of the distortions. We conclude that Ba2CeIrO6 is cubic on average but exhibits small local distortions.

This result is in contrast to an earlier report on a monoclinic structure of Ba2CeIrO6 in Ref. Wakeshima et al. (2000) which was based on powder data and a tiny monoclinic distortion of the metric ($a/b$ = $1.0016$ Wakeshima et al. (2000)) that can result from the broadening of Bragg peaks. Combining the typical rotation of octahedra in the GdFeO3 structure type of a perovskite $AB$O3 with the doubling of the unit cell in the double perovskite results in a monoclinic distortion, $P2_{1}/c$. One may examine the possible instability of Ba2CeIrO6 by calculating the Goldschmidt tolerance factor for perovskites

[TABLE]

with the ionic radii $r_{i}$. For an ideal cubic perovskite, $t_{p}$ = 1. For the hypothetic perovskites BaIrO3 and BaCeO3 this yields $t_{p}$ = 1.06 and 0.94 not indicating sizable bond-length mismatch. The same analysis for distorted Sr2CeIrO6 yields $t_{p}$ = 0.90 and 0.80. Alternatively, one may consider the tolerance factor $t_{dp}$ for an ordered double perovskite $A_{2}BB^{\prime}$O6 with $r_{B}$ in Eq. (1) to be replaced by $(r_{B}+r_{B^{\prime}})/2$. A monoclinic structure is favored for $t_{dp}\leq 0.96$ while values close to 1 point towards a cubic structure Vasala and Karppinen (2015). For Ba2CeIrO6, one finds $t_{dp}$ = $0.991$, supporting a cubic structure.

III Magnetic susceptibility

To explore the magnetism of the local moments in Ba2CeIrO6 we measured the magnetization and the magnetic susceptibility $\chi(T)$. We used an assembly of 20 small single crystals in order to enhance the magnetic signal. The crystals were not aligned because an isotropic magnetic susceptibility is expected in the paramagnetic phase of the (global) cubic structure. As shown in the lower inset of Fig. 2, we observe a field-linear magnetization. The main panel of Fig. 2 shows $\chi(T)$. Its high-temperature behavior essentially follows a Curie-Weiss behavior from 300 K down to about $T_{N}$ = 14 K, where a distinct drop in $\chi(T)$ signals antiferromagnetic ordering. For the quantitative analysis we use

[TABLE]

where $N_{A}$ and $k_{B}$ denote Avogadro’s and Boltzmann’s constant, respectively, $\Theta_{\rm CW}$ is the Curie-Weiss temperature, and the constant $\chi_{0}$ = $\chi_{dia}+\chi_{vV}$ represents core diamagnetism $\chi_{dia}\!\simeq\!-1.7\cdot 10^{-4}$ emu/mol and van Vleck paramagnetism $\chi_{vV}\!>\!0$, which are expected to be of the same order of magnitude 111 $\chi_{dia}\simeq-1.7\cdot 10^{-4}$ emu/mol results from the tabulated values Bain and Berry (2008) for the ions in Ba2CeIrO6.. A fit based on Eq. (2) describes the data above $T_{N}$ very well, see red line in Fig. 2, and yields the parameters $\chi_{0}$ = $1.2\cdot 10^{-4}$ emu/mol, $\mu_{\rm eff}$ = $1.41$ $\mu_{B}$, and $\Theta_{\rm CW}$ = $-184$ K. Very similar values, $\mu_{\rm eff}$ = $1.3$ $\mu_{B}$ and $\Theta_{\rm CW}$ = $-177$ K, were reported in a previous study on polycrystalline Ba2CeIrO6 Wakeshima et al. (2000). To estimate the reliability of our result, we compare with a fit assuming $\chi_{0}$ = 0, which shows Curie-Weiss behavior above about 120 K (blue line in Fig. 2), $\mu_{\rm eff}$ = $1.69$ $\mu_{B}$ and $\Theta_{\rm CW}$ = $-263$ K, i.e., an even larger value of $|\Theta_{\rm CW}|$. Thus, both fits result in an effective magnetic moment that is moderately reduced from $\mu_{\rm eff}$ = $1.73\,\mu_{B}$ expected for $j$ = $1/2$ moments in an ideal cubic crystal field Moretti Sala et al. (2014) and indicate substantial frustration with a frustration parameter $f$ = $\left|\Theta_{\rm CW}\right|/T_{N}\!>\!13$.

IV RIXS

Ba2CeIrO6 indeed realizes nearly ideal local $j$ = $1/2$ moments, which can be inferred from our RIXS results. In cubic symmetry, a single $5d^{5}$ Ir4+ site with a $t_{2g}^{5}$ configuration is expected to show a local $j$ = 1/2 ground state and a $j$ = 3/2 excited state, the so-called spin-orbit exciton, at 1.5 $\lambda$ with $\lambda$ = 0.4-0.5 eV. The effect of a non-cubic crystal field is described by the single-site Hamiltonian

[TABLE]

which shows a crystal-field splitting of the $j$ = $3/2$ quartet and a mixing of $j$ = $1/2$ and $3/2$ wavefunctions in the ground state, $|0\rangle=\alpha\left|\frac{1}{2},\frac{1}{2}\right\rangle+\beta\left|\frac{3}{2},\frac{1}{2}\right\rangle$ in the $|j,j_{z}\rangle$ basis. With $\alpha$ = $(\sin\theta+\sqrt{2}\cos\theta)/\sqrt{3}$ and $\tan 2\theta$ = $\sqrt{8}/(1-2\Delta_{\rm CF}/\lambda)$ Jackeli and Khaliullin (2009) we can readily infer the ground state wavefunction by measuring $\Delta_{\rm CF}$.

To do so, we performed RIXS measurements at the Ir $L_{3}$ edge, the most sensitive probe for the corresponding intra-$t_{2g}$ excitations. For $\Delta_{\rm CF}/\lambda\ll 1$, the experimentally observed peak splitting amounts to $\Delta_{\rm exp}$ = $\frac{2}{3}\Delta_{\rm CF}$. Thus far, all experimental results on the spin-orbit exciton in iridates show a finite non-cubic crystal-field splitting Rau et al. (2016); Trebst (2017); Gretarsson et al. (2013); Kim et al. (2014b); Rossi et al. (2017); Liu et al. (2012); Sala et al. (2014). The smallest values $\Delta_{\rm exp}$ = 0.11-0.14 eV were reported for Rb2IrF6, Na2IrO3, and Sr2IrO4 Rossi et al. (2017); Gretarsson et al. (2013); Kim et al. (2014b). In Rb2IrF6, F-Ir-F bond angles vary from $87^{\circ}$ to $93^{\circ}$ Rossi et al. (2017), while Sr2IrO4 shows distorted IrO6 octahedra with Ir-O bond lengths of 1.98-2.06 Å and Ir-O-Ir bond angles of $157^{\circ}$ Crawford et al. (1994). Despite the substantial distortions, these compounds are widely accepted as realizations of the $j$ = 1/2 scenario. In contrast, strong deviations from the $j$ = 1/2 model are reported for Sr3CuIrO6 and CaIrO3 with $\Delta_{\rm exp}$ = 0.23 eV and 0.6 eV, respectively Liu et al. (2012); Sala et al. (2014).

For Ba2CeIrO6, we measured RIXS data on a polished (0 0 1) surface at the ID20 beamline at ESRF using an incident energy of 11.215 keV with an overall resolution of 25 meV Moretti Sala et al. (2013, 2018). The incident photons were $\pi$ polarized. Our data offer a textbook example of the spin-orbit exciton by showing two narrow RIXS peaks on a negligible background, see Fig. 3. Similar RIXS spectra with a slightly larger peak splitting were reported for Rb2IrF6 Rossi et al. (2017) and Ba3Ti2.7Ir0.3O9 Revelli et al. (2019), two compounds with well separated Ir4+ ions. In comparison, $5d^{5}$ iridates with stronger hopping such as Na2IrO3 and Sr2IrO4 show more complex RIXS features Gretarsson et al. (2013); Kim et al. (2014b) with, e.g., further peaks, broader line widths, and/or a continuum contribution. In Ba2CeIrO6, the peaks are located at about 0.61 eV and 0.71 eV, both at 10 K and at 300 K. The observation of two peaks signals non-cubic local distortions in agreement with our analysis of the x-ray diffraction data. A fit using two peaks with the Pearson VII line shape Wang and Zhou (2005) that mimics a convolution of an intrinsic Lorentzian line shape and a Gaussian profile with the experimental resolution yields a splitting $\Delta_{\rm exp}$ = $(100\pm 4)$ meV, the smallest splitting reported thus far in $L$ edge RIXS for the spin-orbit exciton in iridates Rossi et al. (2017); Gretarsson et al. (2013); Kim et al. (2014b); Liu et al. (2012); Sala et al. (2014); Revelli et al. (2019). The peak values of 0.61 eV and 0.71 eV allow for two different solutions of Eq. (3) with $\lambda$ = 0.43 eV and $\Delta_{\rm CF}$ = 0.17 eV or -0.15 eV, which correspond to elongation or compression, respectively. This results in a ground-state wavefunction

[TABLE]

in the $|j,j_{z}\rangle$ basis for elongation, while for compression the coefficients are 0.995 and 0.100, respectively. Note that both solutions deviate by less than 1 % from the ideal $j$ = 1/2 case.

To probe the intersite hopping interactions, we have measured the dispersion via RIXS for ${\bf q}$ along different high-symmetry directions. Data along $\Gamma$-$K$ and $\Gamma$-$L$ paths reveal a finite dispersion of up to 15-20 meV, while peak energies are nearly independent of ${\bf q}$ along $\Gamma$-$X$, see the lower panels of Fig. 3. The corresponding delocalization of the $j$ = $3/2$ excited state is a clear signature of microscopic hopping processes and intersite interactions that are closely related to the magnetic exchange interactions between localized $j$ = 1/2 moments Kim et al. (2012, 2014b). Roughly, this common microscopic origin is reflected in the common energy scale of 15-20 meV of the spin-orbit exciton dispersion and the Curie-Weiss temperature, which also is a measure of the size of magnetic exchange interactions.

V Microscopic Model

V.1 fcc lattice with cubic site symmetry

A symmetry analysis Cook et al. (2015); Aczel et al. (2016); Li et al. (2017) of exchange interactions on the undistorted fcc lattice shows that the most general nearest-neighbor spin Hamiltonian allows for Heisenberg coupling $J_{1}$, Kitaev coupling $K$, and symmetric off-diagonal exchange $\Gamma$. We estimate the coupling constants using density functional theory (GGA+U+SOC) for different magnetic configurations and $t/U$ perturbation theory for an effective tight-binding model (see Appendix B). Both approaches consistently yield an antiferromagnetic $J_{1}\!\approx\!5-7$ meV and two subdominant couplings $K\!\approx J_{2}\!\approx\!0.2\,J_{1}$, where $J_{2}$ denotes a next-nearest neighbor Heisenberg coupling. We find that $\Gamma/J_{1}\!\lesssim\!0.05$ is negligible. The corresponding Curie-Weiss temperature $\Theta_{\rm CW}$ = $-(3J_{1}+K+3J_{2}/2)\!\approx\!-200$ K to $-280$ K agrees with the experimental $\chi(T)$; see Fig. 2. Note that we find an antiferromagnetic Kitaev coupling, in contrast to the ferromagnetic ones inferred for the honeycomb-based iridates and $\alpha$-RuCl3 Winter et al. (2016). The ferromagnetic Kitaev coupling of the latter arises from Hund’s coupling in the virtually excited intermediate state with two holes on the same site, favoring parallel hole spins. For the honeycomb materials with a $90^{\circ}$ Ir-O-Ir exchange path, this translates into a ferromagnetic coupling of $j$ = 1/2 pseudo-spins. In Ba2CeIrO6, exchange proceeds via an Ir-O-O-Ir path with a different combination of orbitals in the virtual state. Again, Hund’s coupling favors parallel spins of the two holes, but for the relevant orbitals this translates to antiferromagnetic coupling of $j$ = 1/2 pseudo-spins (see Appendix C).

To study the competition of geometric and exchange frustration, we explore the minimal microscopic model

[TABLE]

where $\langle i,j\rangle_{\gamma}$ denotes nearest-neighbor pairs in the plane perpendicular to axis $\gamma$ (= $x$,$y$,$z$), $\langle\langle i,j\rangle\rangle$ runs over next-nearest-neighbor pairs, and the spin operators ${\vec{\cal S}}$ refer to $j$ = $1/2$ moments. We have calculated its rich phase diagram using a pseudofermion functional renormalization group (pf-FRG) approach Reuther and Wölfle (2010). This numerical scheme combines elements from $1/S$ expansion Baez and Reuther (2017) and $1/N$ expansion Buessen et al. (2018a); Roscher et al. (2018), allowing it to capture both magnetic order and spin-liquid ground states. There are four magnetically ordered phases, one of them showing incommensurate spiral order, see Fig. 4a). These phases can be readily understood in the classical limit of model (5) via a Luttinger-Tisza approach Luttinger and Tisza (1946); Luttinger (1951), with the classical phase boundaries also indicated in Fig. 4a). The quantum model additionally exhibits a spin-liquid phase with no magnetic order. Its origin is revealed by two points of special interest in the classical model, see white and black circles in Fig. 4a): (i) $J_{2}$ = $K$ = [math], the fcc nearest-neighbor Heisenberg antiferromagnet. It exhibits a degenerate manifold of coplanar spin spiral ground states Henley (1987). The corresponding set of ${\bf q}$ vectors is shown in Fig. 4b). (ii) $J_{2}$ = $J_{1}/2$, $K$ = [math], where three ordered phases meet in the classical model. This point features an even larger set of degenerate coplanar spin-spiral ground states, depicted by the surface of ${\bf q}$ vectors in Fig. 4c). The presence of a considerable (but still subextensive) manifold of (nearly) degenerate low-energy states appears to give rise to an extended spin liquid regime in the quantum model, centered around the classical high-degeneracy point 222 A similar scenario has recently been discussed in the context of the $J_{1}-J_{2}$ Heisenberg model on the diamond lattice Bergman et al. (2007); Buessen et al. (2018b). .

To further investigate the interplay of geometric and exchange frustration, we calculate Reuther et al. (2011) the dimensionless frustration parameter $f$ = $|\Theta_{\rm CW}|/T_{N}$, see Fig. 4d), using estimates of $\Theta_{\rm CW}$ and $T_{N}$ obtained from fits of the magnetic susceptibility numerically obtained by FRG calculations. The frustration parameter diverges in the spin liquid regime due to the absence of finite-temperature order. Furthermore, $f$ is particularly large along the phase boundary between the $(1\frac{1}{2}0)$ and $(\frac{1}{2}\frac{1}{2}\frac{1}{2})$ phases, where both $J_{2}$ and $K$ are substantial and antiferromagnetic. This boosts $|\Theta_{\rm CW}|$ while $T_{N}$ is small close to the phase boundary. Close to the spin-liquid regime for the parameter set estimated for Ba2CeIrO6 (cf. star in Fig. 4d)), we also find large values of $f$. However, moving away from the spin-liquid regime the frustration is quickly reduced with increasing strength of the Kitaev coupling. This is consistent with a previous classical Monte Carlo study Cook et al. (2015); Aczel et al. (2016), although such a classical analysis by itself is not reliable in the deep quantum limit of $j\!=\!1/2$. Our results show that the Kitaev coupling, in competition with the geometric frustration of the Heisenberg exchange, indeed induces magnetic order for the system at hand – in striking contrast to a number of $j$ = $1/2$ materials where the Kitaev coupling is primarily considered a source of frustration Rau et al. (2016); Trebst (2017).

V.2 Distortions

The strong frustration in Ba2CeIrO6 boosts the importance of magneto-elastic coupling. We find theoretically that even small local distortions severely affect the exchange couplings, although the ground state wavefunction remains close to the $j$ = $1/2$ limit, see Eq. (4). The precise character of the local distortions cannot be determined from our x-ray diffraction results, which show global cubic symmetry. A putative tetragonal distortion of strength $\Delta_{\rm CF}$ gives rise to a strong spatial anisotropy, which can be rationalized as follows. Focusing, e.g., on the dominant contribution to exchange within the $xy$ plane, we find $J_{1}^{xy}$ to depend quadratically on the occupation probability of the $xy$ orbital. Comparing cubic $\Delta_{\rm CF}$ = [math] with the distorted case $\Delta_{\rm CF}/\lambda\approx 0.4$ derived above, the $xy$ occupation is strongly enhanced from $1/3$ to $0.46$ and as a result the nearest-neighbor Heisenberg exchange $J_{1}^{xy}$ increases by about a factor of two, which corresponds to a dramatic magneto-elastic effect. In particular, $\Delta_{\rm CF}\!>\!0$ strengthens (weakens) $J_{1}$, $J_{2}$, and $K$ in the $xy$ plane ($yz$ and $xz$ planes), while $\Delta_{\rm CF}\!<\!0$ has the reverse effect. This strong spatial anisotropy of the couplings is sketched in Fig. 1c). Note that a change of the $xy$ occupation $\sin^{2}\theta$ has a much more pronounced effect on the exchange, $J_{1}\propto\sin^{4}\theta$, than on the coefficient $\alpha$ = $\sqrt{1/3}\sin\theta+\sqrt{2/3}\sqrt{1-\sin^{2}\theta}$ of the $|\frac{1}{2},\frac{1}{2}\rangle$ contribution to the ground state wavefunction. The comparably small change of $\alpha$ indicates a small deviation from a cubic charge distribution and a concomitant small energy cost for lattice distortions, while the larger change of $J_{1}$ and in particular its spatial anisotropy yield a significant gain of magnetic energy, particularly in the presence of frustration.

To analyze the effect of a global tetragonal lattice distortion, we have simulated a variant of the $J_{1}$-$J_{2}$-$K$ model, cf. Eq. (5), with anisotropic coupling strengths, enhancing/reducing the couplings as described above and illustrated in Fig. 1c). Specifically, we have modeled the dependence of the coupling parameters on the distortion $\Delta_{\rm CF}$ as

[TABLE]

and, analogously, for $J_{2}$ and $K$. For $\Delta_{\rm CF}$ = [math] this corresponds to the parameter set indicated by a star in Fig. 4d), i.e. $J_{2}$ = $K$ = $0.2J_{1}$, while for $\Delta_{\rm CF}$ = $\Delta^{\rm max}_{\rm CF}$ = $160$ meV we have the enhanced $xy$ couplings $J_{1}^{xy}$ = $2J_{1}$, $J_{2}^{xy}$ = $1.2J_{2}$, and $K^{xy}$ = $1.2K$ for simultaneously reduced parameters $J_{1}^{xz,yz}$ = $0.6J_{1}$, $J_{2}^{xz,yz}$ = $0.7J_{2}$, and $K^{xz,yz}$ = $0.8K$ in the $xz$ and $yz$ plane. The function $\delta J(\Delta_{\rm CF})$ is well-approximated as a linear interpolation

[TABLE]

Results for pf-FRG calculations are summarized in Fig. 5 showing the frustration parameter $f=|\Theta_{\rm CW}|/T_{\rm N}$ as a function of the distortion $\Delta_{\rm CF}$. As clearly visible, the frustration is strongly suppressed by the distortions, it quickly approaches a non-frustrated regime $f\lesssim 5$ for distortions of the order of $|\Delta_{\rm CF}|\approx 40$ meV, independent of the sign of the distortion. Additionally, the data indicate a potential change of magnetic order, for instance to (100) order, depending on the sign and the strength of $\Delta_{\rm CF}$. These results agree with the small $f\lesssim 2$ reported for the globally distorted monoclinic $5d^{5}$ double perovskites La2ZnIrO6 and La2MgIrO6 Aczel et al. (2016). However, this result for a global distortion is significantly smaller than the value of $f\!>\!13$ measured in Ba2CeIrO6. This suggests that the statistical distribution of local distortions along different tetragonal axes, in contrast to a global distortion, is important in order to recover the experimentally observed large frustration.

The existence of a weak (but unavoidable) magneto-elastic effect was recently discussed for Sr2IrO4 Liu and Khaliullin (2019). The strong effect of the magneto-elastic coupling in Ba2CeIrO6, however, is due to an additional mechanism arising from an interplay of distortions and magnetic frustration, which is not present in Sr2IrO4, but will be relevant, e.g., in tetragonal bilayer Sr3Ir2O7 and in the 3D honeycomb iridates.

VI Conclusions

The spin-orbit entangled $j$ = 1/2 wavefunction has proved to be a versatile source for novel states of quantum matter. Its experimental realization in the double perovskite Ba2CeIrO6 deviates less than 1 % from the ideal theoretical scenario for a cubic system and is one of the most pristine $j$ = 1/2 incarnations reported so far in the literature. Combining structural analysis, magnetic susceptibility measurements, and RIXS with ab initio and functional renormalization group calculations for the obtained microscopic model Hamiltonian we find that the collective magnetism of this fcc compound is governed by a competition of geometrical frustration, Kitaev-type bond-directional exchange, and magneto-elastic coupling. In striking contrast to the honeycomb-based Kitaev materials, the Kitaev exchange is antiferromagnetic and in fact stabilizes long-range magnetic order in proximity to a spin liquid phase. Importantly, the exchange couplings turn out to be highly sensitive to small deviations from cubic symmetry, giving rise to a dramatic magneto-elastic coupling. This should be contrasted with the common notion that $j$ = 1/2 moments are not Jahn-Teller active, as the orbital degeneracy is lifted by spin-orbit coupling. The strong magneto-elastic coupling resurrects the prominent role of lattice distortions on the low-energy properties of $j$ = 1/2 compounds.

Note added: After submission of our manuscript, Aczel et al. reported similar experimental results on polycrystalline samples of Ba2CeIrO6 Aczel et al. (2019). They find (100) magnetic order, in agreement with our calculations for negative $\Delta_{\rm CF}$. Furthermore, Khan et al. Khan et al. (2019) reported on the realization of cubic site symmetry in K2IrCl6 based on x-ray diffraction data; a spectroscopic proof of this claim is still missing.

Acknowledgements.

We acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project numbers 277146847 and 247310070 – CRC 1238 (projects A02, B01, B02, B03, C02, C03) and CRC 1143 (project A05), respectively. M.H. acknowledges partial funding by the Knut and Alice Wallenberg Foundation and the Swedish Research Council. A.P. acknowledges support from NSERC of Canada and the Canadian Institute for Advanced Research, and the support and hospitality of the University of Cologne. The numerical simulations were performed on the JURECA booster at FZ Jülich and the CHEOPS cluster at RRZK Cologne. DFT calculations were supported by the Russian Science Foundation via project 17-12-01207.

Appendix A Single-crystal growth and characterization

We have grown single crystals of Ba2CeIrO6 by melt solution growth using BaCl2 as melt solvent and BaCO3 (Merck, p.a.), IrO2 (Chempur, 99.9 %) and CeO2 (Auer Remy, 99.9 %) as educts for the crystals. Due to the moderate solubility of metal oxides in halide melts Elwell and Scheel (1975), a ratio flux/crystal of 15/1 was used to achieve sufficient dissolution of the oxides. The crucible was sealed with a lid to prevent evaporation of BaCl2 from the melt solution. Within three weeks single crystals of about 1 mm3 size were obtained. The black single crystals were separated from the flux by dissolving the flux in deionized water and analyzed by energy-dispersive x-ray scattering.

Our single-crystal x-ray diffraction data strongly support a cubic structure of Ba2CeIrO6, as explained in the main text. Figure A1 shows the distribution of the observed peak intensities $I$ divided by their error bars $\sigma(I)$ for all reflections that are not allowed in space group $Fm\bar{3}m$. The width of the distributions is approximately 1, which is an indication for meaningful statistical errors of the intensities. Both at room temperature and at 100 K, we find a Gaussian profile peaking at a value of $I/\sigma(I)$ close to zero or much smaller than 1, i.e., the intensities of forbidden peaks are negligible within the experimental error bars.

The excellent agreement of the RIXS spectra shown in Fig. 3 with the expectations for the spin-orbit exciton provides an unambiguous fingerprint of the Ir4+ valence state. RIXS spectra for $5d^{4}$ Ir5+ are distinctly different, as reported for, e.g., the double perovskites Sr2YIrO6 and Ba2YIrO6 Yuan et al. (2017); Nag et al. (2018). In Ba2CeIrO6, the Ir4+ valence means that also the Ce ions are tetravalent, as claimed before Wakeshima et al. (2000) based on the dependence of the lattice parameters on the ionic radius of the lanthanide ions $Ln$ in Ba2$Ln$IrO6.

The Mott-insulating character of Ba2CeIrO6 is demonstrated by the very low value of the optical conductivity in the mid-infrared range, $\sigma_{1}(\omega)\!\approx\!1$ ($\Omega$cm)-1, see Fig. A2. Using a Bruker IFS 66/v Fourier-transform infrared spectrometer, we measured the transmittance on a single crystal with a thickness of ($30\pm 5$) $\mu$m. On the low-frequency side, the accessible frequency range is cut off by strong phonon absorption suppressing the transmittance. The steep edge in $\sigma_{1}(\omega)$ at about 0.1 eV corresponds to the upper limit for single phonon absorption, while the weak features up to about 0.2 eV can be attributed to multi-phonon absorption. On the high-frequency side, the transmittance is suppressed by electron-hole excitations across the gap, giving rise to the increase of $\sigma_{1}(\omega)$ above about 0.2 eV.

Appendix B Derivation of microscopic model

*Ab initio calculations.– * In order to calculate Heisenberg-type and Kitaev-type exchange constants, we used the projector augmented-wave (PAW) method Blöchl (1994) as realized in the pseudopotential VASP code Kresse and Furthmüller (1996). The exchange-correlation potential was chosen in the form proposed by Perdew, Burke, and Ernzerhof Perdew et al. (1996). Electronic correlations and spin-orbit coupling were considered in the framework of the GGA+U+SOC formalism Dudarev et al. (1998) with $U\!-\!J_{H}$ = $2.0$ eV Kim et al. (2008), where $U$ and $J_{H}$ denote the on-site Coulomb repulsion and intra-atomic Hund’s exchange, respectively. The integration was performed on a $7\times 7\times 7$ mesh of the Brillouin zone. We calculated total energies of three magnetic configurations (FM, AFM type I and type II) to extract nearest and next-nearest neighbor exchange constants $J_{1}$ and $J_{2}$. The Kitaev coupling $K$ was computed via the difference of the total energies of two AFM type I configurations with spins pointing along $a$ and $c$, respectively.

*Perturbative approach.– * For the perturbative calculation of exchange couplings, we use one-hole and two-hole eigenstates and energies obtained by numerical diagonalization of the single-site Hamiltonian

[TABLE]

with $\lambda$ = $0.43$ eV deduced from our RIXS data, Hund’s coupling $J_{H}$ = $0.25$ eV Yuan et al. (2017) and $U$ = $2.25$ eV as a typical estimate for the Hubbard repulsion (which leads to $U\!-\!J_{H}$ = $2.0$ eV as in the ab initio calculations). For the inter-site Hamiltonian, we use $H_{\rm hop}$ = $\sum_{\langle i,j\rangle,\sigma}{\bf d}^{\dagger}_{i\sigma}\cdot{\rm\bf T}^{\phantom{\dagger}}_{ij}\cdot{\bf d}^{\phantom{\dagger}}_{j\sigma}$, with ${\bf d}^{\dagger}$ representing the creation operators for the three $t_{2g}$ orbitals, $(d^{\dagger}_{yz},d^{\dagger}_{xz},d^{\dagger}_{xy})$. For the hopping amplitudes $t_{\alpha-\beta}$ between a pair of orbitals ($\alpha,\beta$) on nearest-neighbor sites $(i,i+\hat{x}+\hat{y})$, we retain only the by far dominant one $t_{xy-xy}\!\approx\!-150$ meV, which is an order of magnitude larger than $t_{xz-yz}$ which again is larger than $t_{xz-xz}$ = $t_{yz-yz}$, while others vanish by symmetry. For second-neighbor pairs $(i,i+2\hat{x})$, we employ $t_{xy-xy}$ = $t_{xz-xz}\!\approx\!30$ meV, while other hopping amplitudes vanish by symmetry. For the cubic fcc lattice, the corresponding hopping amplitudes between all other nearest-neighbor or next-nearest neighbor pairs are determined by symmetry. To extract the two-site exchange Hamiltonian, we carry out second-order degenerate perturbation theory in $t/U$, evaluating matrix elements numerically using exact single-site eigenfunctions and energies. The resulting exchange Hamiltonian has dominant nearest-neighbor Heisenberg exchange interaction $J_{1}$ with subdominant Kitaev and second-neighbor Heisenberg terms as listed in the main text and a negligible $\Gamma\lesssim 0.05J_{1}$ exchange term.

Appendix C Antiferromagnetic character of Kitaev coupling

The $j$ = $1/2$ wave function is given by

[TABLE]

We address superexchange interactions between two sites A and B in the $xy$ plane. For comparison, we first consider edge-sharing geometry with $90^{\circ}$ Ir-O-Ir bonds Jackeli and Khaliullin (2009) as approximately realized in the honeycomb iridates. In this case, the hopping $t_{xy-xy}^{90^{\circ}}$ between $xy_{A}$ and $xy_{B}$ via an intermediate O ligand vanishes by symmetry. There are two finite hopping contributions, between $xz_{A}$ and $yz_{B}$ and between $yz_{A}$ and $xz_{B}$. The Heisenberg interaction vanishes due to destructive interference between these two. This destructive interference originates from the phase factor $i$ in the $j$ = $1/2$ wavefunction, see Eq. (A9). In contrast, Kitaev exchange remains finite and is ferromagnetic. The ferromagnetic character arises from the virtual intermediate state with two holes on the same site occupying the $xz$ and $yz$ orbitals. The energy of this intermediate state is lower for parallel spins in $xz$ and $yz$, which corresponds to parallel $j$ = $1/2$ pseudo-spins, see Eq. (A9).

The same ferromagnetic contribution to Kitaev exchange is present on the $fcc$ lattice as well where the Ir-O-Ir hopping of edge-sharing geometry has to be replaced by Ir-O-O-Ir hopping. However, the size of this ferromagnetic contribution to Kitaev exchange is negligibly small on the $fcc$ lattice due to the very small hopping between $xz$ and $yz$ orbitals. In contrast to the edge-sharing geometry discussed above, the hopping $xy_{A}\rightarrow xy_{B}$ via the two intermediate O ligands does not vanish on the $fcc$ lattice, it is in fact the dominant hopping term within the $xy$ plane. For parallel pseudo-spins $|\frac{1}{2},{\rm up}\rangle_{A}\,|\frac{1}{2},{\rm up}\rangle_{B}$, the configuration $|xy,\uparrow\rangle_{A}\,|yz,\downarrow\rangle_{B}$ in orbital/spin basis carries finite weight, see Eq. (A9). Hopping $xy_{A}\rightarrow xy_{B}$ then yields the doubly occupied virtual state $|xy,\uparrow\rangle_{B}\,|yz,\downarrow\rangle_{B}$, and a second hopping process $xy_{B}\rightarrow xy_{A}$ brings us back to the ground state without any spin flip. This exchange will lead to an Ising-like $S_{A}^{z}S_{B}^{z}$ term. Equivalent to the case discussed above, the lowest energy in the virtual doubly occupied state is realized for two parallel spins, $|xy,\uparrow\rangle_{B}\,|yz,\uparrow\rangle_{B}$, which for this combination of orbitals corresponds to antiparallel pseudo-spins, see Eq. (A9). As a result, pure $xy_{A}\rightarrow xy_{B}$ hopping within the $xy$ plane in combination with Hund’s coupling generates a small antiferromagnetic Kitaev coupling $K\propto J_{1}\cdot J_{H}/U$. By symmetry, it couples the $z$ ($x$, $y$) component of the moments for two nearest-neighbor sites within the $xy$ ($yz$, $zx$) plane.