Dynamics of a j=3/2 quantum spin liquid

TL;DR

This paper investigates a $j=3/2$ quantum spin liquid in double perovskites, proposing a chiral spin-orbital liquid with Majorana fermions, and characterizes its thermodynamic and spectroscopic properties using theoretical models.

Contribution

It introduces a novel chiral spin-orbital liquid model for $j=3/2$ systems and analyzes its experimental signatures, including RIXS cross sections, in relation to Kitaev-like physics.

Findings

Calculated heat capacity and neutron scattering spectra match experimental data for Ba$_{2}$YMoO$_{6}$.

Proposed RIXS operators can selectively probe pseudospin and pseudo-orbital excitations.

Identified signatures of Majorana fermions in spectroscopic measurements.

Abstract

We study a spin-orbital model for 4$d^{1}$ or 5$d^{1}$ Mott insulators in ordered double perovskites with strong spin-orbit coupling. This model is conveniently written in terms of pseudospin and pseudo-orbital operators representing multipoles of the effective $j=3/2$ angular momentum. Similarities between this model and the effective theories of Kitaev materials motivate the proposal of a chiral spin-orbital liquid with Majorana fermion excitations. The thermodynamic and spectroscopic properties of this quantum spin liquid are characterized using parton mean-field theory. The heat capacity, spin-lattice relaxation rate, and dynamic structure factor for inelastic neutron scattering are calculated and compared with the experimental data for the spin liquid candidate Ba$_{2}$YMoO$_{6}$. Moreover, based on a symmetry analysis, we discuss the operators involved in resonant inelastic X-ray scattering (RIXS) amplitudes for double perovskite compounds. In general, the RIXS cross sections allow one to selectively probe pseudospin and pseudo-orbital degrees of freedom. For the chiral spin-orbital liquid in particular, these cross sections provide information about the spectrum for different flavors of Majorana fermions.

Tables

Table 1. Table 1: Operators describing active multipoles within a cubic Γ 8 subscript Γ 8 \Gamma_{8} quartet. Bars over functions of 𝐉 𝐉 \mathbf{J} indicate the symmetrization with respect to all the possible permutations of the indices, e.g., J x ​ J y 2 ¯ = J x ​ J y 2 + J y ​ J x ​ J y + J y 2 ​ J x ¯ subscript 𝐽 𝑥 superscript subscript 𝐽 𝑦 2 subscript 𝐽 𝑥 superscript subscript 𝐽 𝑦 2 subscript 𝐽 𝑦 subscript 𝐽 𝑥 subscript 𝐽 𝑦 superscript subscript 𝐽 𝑦 2 subscript 𝐽 𝑥 \overline{J_{x}J_{y}^{2}}=J_{x}J_{y}^{2}+J_{y}J_{x}J_{y}+J_{y}^{2}J_{x} . Adapted from Refs. Santini  et al. ( 2009 ); Shiina  et al. ( 1998 ) .

Moment	Symmetry	J multipoles	$(\text{s},𝝉)$ multipoles
Dipoles	${\Gamma }_4$	$J^x$	$s^x\left(1+4{\tau }^{yz}\right)$
		$J^y$	$-s^y\left(1+4{\tau }^{xz}\right)$
		$J^z$	$s^z\left(1+4{\tau }^{xy}\right)$
Quadrupoles	${\Gamma }_3$		$6{\tau }^z$
		$O_{x^2-y^2}={\left(J^x\right)}^2-{\left(J^y\right)}^2\equiv {\hat{O}}_2^2$	$2\sqrt{3}{\tau }^x$
	${\Gamma }_5$	$O_{xy}=\frac{1}{2}\overline{J^x{J}^y}\equiv {\hat{O}}_2^{-2}$	$2\sqrt{3}{s}^z{\tau }^y$
		$O_{yz}=\frac{1}{2}\overline{J^y{J}^z}\equiv {\hat{O}}_2^{-1}$	$2\sqrt{3}{s}^x{\tau }^y$
		$O_{xz}=\frac{1}{2}\overline{J^x{J}^z}\equiv {\hat{O}}_2^1$	$-2\sqrt{3}{s}^y{\tau }^y$
Octupoles	${\Gamma }_2$	$T_{xyz}=\frac{\sqrt{15}}{6}\overline{J^x{J}^y{J}^z}$	$\frac{3\sqrt{5}}{2}{\tau }^y$
	${\Gamma }_4$	$T_x^{\alpha }={\left(J^x\right)}^3-\frac{1}{2}\left(\overline{J^x{\left(J^y\right)}^2}+\overline{{\left(J^z\right)}^2{J}^x}\right)$	$3s^x\left(1-{\tau }^{yz}\right)$
		$T_y^{\alpha }={\left(J^y\right)}^3-\frac{1}{2}\left(\overline{J^y{\left(J^z\right)}^2}+\overline{{\left(J^x\right)}^2{J}^y}\right)$	$-3s^y\left(1-{\tau }^{xz}\right)$
		$T_z^{\alpha }={\left(J^z\right)}^3-\frac{1}{2}\left(\overline{J^z{\left(J^x\right)}^2}+\overline{{\left(J^y\right)}^2{J}^z}\right)$	$3s^z\left(1-{\tau }^{xy}\right)$
	${\Gamma }_5$	$T_x^{\beta }=\frac{\sqrt{15}}{6}\left[\overline{J^x{\left(J^y\right)}^2}-\overline{{\left(J^z\right)}^2{J}^x}\right]$	$3\sqrt{5}{s}^x{\overline{\tau }}^{yz}$
		$T_y^{\beta }=\frac{\sqrt{15}}{6}\left[\overline{J^y{\left(J^z\right)}^2}-\overline{{\left(J^x\right)}^2{J}^y}\right]$	$-3\sqrt{5}{s}^y{\overline{\tau }}^{xz}$
		$T_z^{\beta }=\frac{\sqrt{15}}{6}\left[\overline{J^z{\left(J^x\right)}^2}-\overline{{\left(J^y\right)}^2{J}^z}\right]$	$3\sqrt{5}{s}^z{\overline{\tau }}^{xy}$

Table 2. Table 2: Coefficients of the scattering operators in Eq. 124 as a function of the angle parameter φ 𝜑 \varphi . The columns with the cubic limit values are obtained by taking φ = arcsin ⁡ ( 1 / 3 ) 𝜑 1 3 \varphi=\arcsin(1/\sqrt{3}) .

	$L_3$ edge	Cubic	$L_2$ edge	Cubic
$a_{\mu ,U}$	$\frac{1}{9}\left({\cos }^2\phi -\sqrt{2}\sin 2\phi +3{\sin }^2\phi -1\right)$	0	$\frac{1}{9}\left({\cos }^2\phi +\sqrt{2}\sin 2\phi -2\right)$	0
$a_{\mu ,Q_2}$	$\frac{\sqrt{3}}{18}\left({\cos }^2\phi -\sqrt{2}\sin 2\phi +2-6{\sin }^2\phi \right)$	0	$\frac{\sqrt{3}}{18}\left({\cos }^2\phi +\sqrt{2}\sin 2\phi +4\right)$	$\frac{\sqrt{3}}{3}$
$a_{\mu ,Q_3}$	$2\sqrt{2}\left(\cos \phi -\sqrt{2}\sin \phi \right)$	0	$2\left(\sqrt{2}\cos \phi +\sin \phi \right)$	$\frac{\sqrt{3}}{3}$
$a_{\mu ,T_x}$	$-\frac{2\sqrt{2}}{3}\left(\cos \phi -\sqrt{2}\sin \phi \right)$	0	$-\frac{4}{3}\left(\sqrt{2}\cos \phi +\sin \phi \right)$	$-\frac{4\sqrt{3}}{3}$
$a_{\mu ,T_y}$	$-\frac{2\sqrt{2}}{3}\left(\cos \phi -\sqrt{2}\sin \phi \right)$	0	$-\frac{4}{3}\left(\sqrt{2}\cos \phi +\sin \phi \right)$	$-\frac{4\sqrt{3}}{3}$
$a_{\mu ,T_z}$	$\frac{4\sqrt{2}}{3}\left(\cos \phi -\sqrt{2}\sin \phi \right)$	0	$-\frac{4}{3}\left(\sqrt{2}\cos \phi +\sin \phi \right)$	$-\frac{4\sqrt{3}}{3}$
$a_{\mu ,P_x}$	$\frac{1}{3}\left(2+\frac{3\sqrt{2}}{2}\sin 2\phi \right)$	$\frac{4}{3}$	$-\frac{2}{3}$	$-\frac{2}{3}$
$a_{\mu ,P_y}$	$\frac{1}{3}\left(2+\frac{3\sqrt{2}}{2}\sin 2\phi \right)$	$\frac{4}{3}$	$-\frac{2}{3}$	$-\frac{2}{3}$
$a_{\mu ,P_z}$	$\frac{1}{3}\left(4{\cos }^2\phi +\sqrt{2}\sin 2\phi \right)$	$\frac{4}{3}$	$-\frac{1}{3}\left({\cos }^2\phi +\sqrt{2}\sin 2\phi \right)$	$-\frac{2}{3}$
$b_{\mu ,P_x}$	$\frac{2\sqrt{2}}{3}\left(\cos \phi -\sqrt{2}\sin \phi \right)$	0	$\frac{4}{3}\left(\sqrt{2}\cos \phi +\sin \phi \right)$	$\frac{4\sqrt{3}}{3}$
$b_{\mu ,P_y}$	$-\frac{2\sqrt{2}}{3}\left(\cos \phi -\sqrt{2}\sin \phi \right)$	0	$-\frac{4}{3}\left(\sqrt{2}\cos \phi +\sin \phi \right)$	$-\frac{4\sqrt{3}}{3}$
$b_{\mu ,P_z}$	0	0	$0$	0
$c_{\mu ,P_x}$	$\frac{1}{3}\left(4-3\sqrt{2}\sin 2\phi \right)$	0	$-\frac{4}{3}$	$-\frac{4}{3}$
$c_{\mu ,P_y}$	$\frac{1}{3}\left(4-3\sqrt{2}\sin 2\phi \right)$	0	$-\frac{4}{3}$	$-\frac{4}{3}$
$c_{\mu ,P_z}$	$-\frac{2}{3}\left(\sqrt{2}\sin 2\phi -4{\sin }^2\phi \right)$	0	$\frac{1}{3}\left(5+\cos 2\phi +2\sqrt{2}\sin 2\phi \right)$	$\frac{8}{3}$

Equations

$$\text{l}=-{\mathcal{P}}_{t_{2g}}\text{L}{\mathcal{P}}_{t_{2g}},$$

$$d_{0,\sigma }$$

$$d_{\pm 1,\sigma }$$

$$H_{\text{ion}}=-\lambda \text{l}\cdot \text{S},$$

$$A_{\sigma }$$

$$B_{\sigma }$$

$$C_{\sigma }$$

$$|B_{\uparrow }⟩$$

$$|B_{\downarrow }⟩$$

$$|C_{\uparrow }⟩$$

$$|C_{\downarrow }⟩$$

$$s=$$

$$\tau =$$

$$|B_{\uparrow }⟩$$

$$|B_{\downarrow }⟩$$

$$|C_{\uparrow }⟩$$

$$|C_{\downarrow }⟩$$

$$J^z=s^z+4s^z{\tau }^z,$$

$$s^z|s,\tau ⟩$$

$${\tau }^z|s,\tau ⟩$$

$$s^a$$

$${\tau }^a$$

$$s^a{\tau }^b$$

$$[s^a,s^b]$$

$$[{\tau }^a,{\tau }^b]$$

$$[s^a,{\tau }^b]$$

$${\tau }^{xy}$$

$${\tau }^{yz(zx)}$$

$${\overline{\tau }}^{xy}$$

$${\overline{\tau }}^{yz(zx)}$$

$$T^{-1}\mathbf{\tau }T=({\tau }^x,-{\tau }^y,{\tau }^z).$$

$$H=$$

$$+U\sum _{i,a}{n}_{i,\alpha \beta ,\uparrow }{n}_{i,\alpha \beta ,\downarrow }-J_H\sum _{i,\sigma }\sum _{\gamma \delta <\alpha \beta }{n}_{i,\alpha \beta ,\sigma }{n}_{i,\gamma \delta ,\sigma }$$

$$+(U-2J_H)\sum _{i,\sigma ,{\sigma }^{\prime }}\sum _{\gamma \delta <\alpha \beta }{n}_{i,\alpha \beta ,\sigma }{n}_{i,\gamma \delta ,{\sigma }^{\prime }}$$

$$-J_H\sum _{i,\gamma \delta <\alpha \beta }\left(d_{i,\alpha \beta ,\uparrow }^{†}{d}_{i,\alpha \beta ,\downarrow }{d}_{i,\gamma \delta ,\downarrow }^{†}{d}_{i,\gamma \delta ,\uparrow }+\text{h.c.}\right)$$

$$-J_H\sum _{i,\gamma \delta <\alpha \beta }\left(d_{i,\alpha \beta ,\uparrow }^{†}{d}_{i,\gamma \delta ,\downarrow }{d}_{i,\alpha \beta ,\downarrow }^{†}{d}_{i,\gamma \delta ,\uparrow }+\text{h.c.}\right).$$

$$H_{\text{so}}=$$

$$-J^{\prime }\sum _{⟨ij{⟩}_{\gamma }}{\text{S}}_i\cdot {\text{S}}_j{P}_{ij}^{(\gamma )}+\frac{3}{2}{J}^{\prime }\sum _{⟨ij{⟩}_{\gamma }}{n}_{i,\alpha \beta }{n}_{j,\alpha \beta }$$

$$+V\sum _{⟨ij{⟩}_{\gamma }}{n}_{i,\alpha \beta }{n}_{j,\alpha \beta },$$

$$J$$

$$J^{\prime }$$

$$V$$

$$H_{\text{eff}}={\mathcal{P}}_{3/2}{H}_{\text{so}}{\mathcal{P}}_{3/2}.$$

Full text

I Introduction

Quantum spin liquids (QSLs) are highly entangled phases of matter arising in strongly interacting spin systems Savary and Balents (2017). Their intrinsic nonlocal character makes them elusive, since standard experimental techniques probe two-point correlation functions. Experiments performed on QSL candidates must then combine the outcomes of different techniques with a careful theoretical analysis Savary and Balents (2017); Lacroix et al. (2011); Avella and Mancini (2014). The difficulties to experimentally verify these quantum states of matter highlight the importance of studying effective Hamiltonians which stabilize them as ground states. If one could calculate the response functions for these Hamiltonians, general properties of QSLs could be investigated accurately, thus guiding the design and interpretation of experiments.

The Kitaev model on the honeycomb lattice plays an important role in this context Kitaev (2006). This spin Hamiltonian displays Ising-like interactions along different quantization axes depending on the bond directions. This causes an exchange frustration that drives the system to a Majorana QSL ground state Baskaran et al. (2007). Thanks to its integrability, several thermodynamic Nasu et al. (2014, 2015); Nasu and Motome (2016); Yoshitake et al. (2016) and spectroscopic Knolle et al. (2014, 2015); Smith et al. (2015, 2016); Halász et al. (2016) responses of the Kitaev model have been calculated exactly.

Remarkably, the seminal work of Jackeli and Khaliullin Jackeli and Khaliullin (2009) showed that the Kitaev model is a good starting point to describe the magnetism of certain $4d^{5}$ or 5$d^{5}$ Mott insulators. To derive the Kitaev Hamiltonian, they considered the interplay of $t_{2g}$ orbital directionality, hole virtual transfer through intermediate oxygen orbitals, electronic correlation and strong spin-orbit coupling (SOC) Khaliullin (2005). The work in Ref. Jackeli and Khaliullin (2009) made specific proposals for candidate materials that could exhibit a QSL ground state, leading to a manifold of experimental studies, exemplified by Refs. Singh et al. (2012); Plumb et al. (2014); Takayama et al. (2015); Banerjee et al. (2016). Unfortunately, none of the compounds studied so far harbors a Majorana QSL, showing instead different types of magnetic order at low temperatures. The magnetic order in these materials can be explained by the effects of competing exchange interactions which have to be added to the Kitaev model Chaloupka et al. (2010); Rau et al. (2014); Winter et al. (2016). The effective Hamiltonians generated by the Jackeli-Khaliullin mechanism are examples of quantum compass models, which are known to host unusual magnetism Nussinov and van den Brink (2015). The wealth of theoretical proposals and experiments has led much of the research on QSLs to turn to compounds that combine strong correlations and SOC Witczak-Krempa et al. (2014); Schaffer et al. (2016); Winter et al. (2017).

Mott insulators in ordered double perovskites based on heavy $d$ ions satisfy the conditions leading to quantum compass models. Ordered double perovskites are oxides of general stoichiometry A2BB’O6, where A corresponds to an alkaline-earth or lanthanide, and B, B’ are transition metal ions (Fig. 1a). Chen et al. Chen et al. (2010) put forward a spin-orbital model for compounds in which B’ is the only magnetically active ion in a 4$d^{1}$ or 5$d^{1}$ electronic configuration. In materials that retain cubic symmetry, the spin and orbital angular momenta of the electron in the $t_{2g}$ orbital combine to form an effective $j=3/2$ magnetic moment. The effective spin Hamiltonian in this case contains bond-dependent anisotropic interactions between $j=3/2$ moments distributed on an fcc lattice (Fig. 1b). However, in contrast to the Kitaev model, the interactions involve higher multipoles of the angular momentum.

The combination of geometric frustration in the fcc lattice and multipolar interactions induced by SOC can favor exotic phases such as valence bond solids or QSLs Chen et al. (2010). In coherence with these predictions, experimental results show that the double perovskite Ba2YMoO6 does not present any structural transition or magnetic order down to 2K (much lower than its Curie-Weiss temperature) de Vries et al. (2010); Aharen et al. (2010); Carlo et al. (2011); de Vries et al. (2013). Motivated by these observations, a chiral spin-orbital liquid has been proposed as a possible ground state of the double perovskite model in a particular regime Natori et al. (2016). This QSL is similar to some three-dimensional versions of the Kitaev model Hermanns et al. (2015); O’Brien et al. (2016), as it exhibits Majorana fermion excitations with a gapless nodal-line spectrum instead of a Fermi surface. Another theoretical proposal to explain the properties of Ba2YMoO6 is the disordered dimer-singlet phase Romhányi et al. (2017). The latter shares with the chiral spin-orbital liquid the property of pseudo-gapped low-energy excitations, which are however due to a random distribution of dimerized bonds.

In this work, we calculate various response functions of the $j=3/2$ chiral spin-orbital liquid proposed in Ref. Natori et al. (2016) within the mean-field approximation of free Majorana fermions. We calculate the magnetic specific heat, spin-lattice relaxation rate, and inelastic neutron scattering (INS) cross section and compare the theoretical results with the experimental data for Ba2YMoO6 de Vries et al. (2010); Aharen et al. (2010); Carlo et al. (2011); de Vries et al. (2013). We also investigate the expected resonant inelastic X-ray scattering (RIXS) cross sections of the chiral spin-orbital liquid. RIXS has grown in importance in recent years Haverkort (2010); Ament et al. (2011), due to its ability of probing orders that are hidden from neutron experiments Savary and Senthil (2015) and of distinguishing different types of excitations by tuning the polarization and energy of the photons. In fact, recent studies of RIXS cross sections for the Kitaev honeycomb model showed that they can probe gapless Majorana fermions and gapped visons separately Halász et al. (2016). This result is in sharp contrast to the dynamic structure factor measured in neutron scattering, which probes spin excitations that in the Kitaev model always excite a gapped vison Knolle et al. (2014, 2015); Smith et al. (2015, 2016). Therefore, RIXS can give more information about Majorana fermions than the usual experiments. This is particularly interesting for the chiral spin-orbital liquid Natori et al. (2016), whose spectrum contains non-degenerate flavors of Majorana fermions associated with different spin-orbital excitations.

The paper is organized as follows. In Section II, we describe the electronic structure of the Mo5+ ion, relevant for the magnetism in Ba2YMoO6, and derive the microscopic Hamiltonian for ordered double perovskites. Section III discusses the parton mean-field theory for the chiral spin-orbital liquid state. In Section IV, we present our results for specific heat, spin-lattice relaxation rate and INS cross section, providing comparison with available experimental data. Section V contains our results for RIXS scattering operators for ordered double perovskites, based on a symmetry analysis of the $L$ absorption edge. These results apply in general to 4$d^{1}$ and 5$d^{1}$ based compounds. An important outcome of this analysis is a proposal of how to directly probe pseudospin and pseudo-orbital degrees of freedom of the $j=3/2$ multiplet. We apply these results in particular to calculate the RIXS cross sections of the chiral spin-orbital liquid. Finally, we summarize our results and suggest future developments for theory and experiments in Sec. VI. Technical details of the calculations and complementary results are left to the Appendices.

II Model and Symmetry

II.1 $t_{2g}$ orbitals in $d^{1}$ configuration

We start by discussing the orbital physics of singly occupied $t_{2g}$ orbitals. Double perovskites with stoichiometry A2BB’O6 are structurally formed by corner-sharing BO6 and B’O6 octahedra, arranged as shown in Fig. 1a. The projection of the angular momentum $\mathbf{L}$ ($L=2$ for $d$ orbitals) onto the $t_{2g}$ triplet defines a $l=1$ effective angular momentum $\mathbf{l}$ Fazekas (1999); Khomskii (2014):

[TABLE]

in which $\mathcal{P}_{t_{2g}}$ is the projection operator. Let $d_{\alpha\beta,\sigma}$ be the annihilation operator for an electron in the $\alpha\beta$ orbital $(\alpha\beta\in\{xy,yz,xz\})$ with spin $\sigma\in\{\uparrow,\downarrow\}$, and $d_{m_{l},\sigma}$ the corresponding operators for eigenstates of $l^{z}$, with eigenvalue $m_{l}\in\{-1,0,1\}$. The relation between these operators is Chen et al. (2010)

[TABLE]

Equations (2) and (3) provide a complete basis to describe the physics of $d^{1}$ strongly correlated systems. Spin-orbital models for double perovskites considering all states spanned by this basis were studied in Ref. Svoboda et al. (2017). These general models interpolate between the weak and strong SOC limits. Here we focus on the limit in which the SOC is strong enough to justify a projection of the Hamiltonian onto a low-energy subspace. The ionic spin-orbit Hamiltonian is written as

[TABLE]

in which $\mathbf{S}$ is the electronic spin and $\lambda>0$ is the SOC constant. The effect of $H_{\text{ion}}$ is to split the $t_{2g}$ levels into one $j=1/2$ and one $j=3/2$ manifold ($\textbf{J}=\textbf{l}+\textbf{S}$), the latter being energetically favored by a gap of $3\lambda/2$.

It is convenient to organize the six eigenstates of $H_{\text{ion}}$ into three Kramers’ doublets. We define the corresponding annihilation operators by Kim et al. (2012); Chaloupka and Khaliullin (2016)

[TABLE]

where $\sigma=\uparrow,\downarrow=\pm 1/2$ distinguishes between Kramers-degenerate states. Note that in this notation the index $\sigma$ in $A_{\sigma}$ and $B_{\sigma}$ is not directly connected with the actual spin eigenvalue in the $d$ operators on the right-hand side of Eqs. (5).

Figure 2 shows the level splitting and the electronic density profiles of the $A$, $B$ and $C$ states. The $A$ states are associated with the higher-energy $j=1/2$ subspace. In the presence of cubic symmetry (which is the case for Ba2YMoO6 Chen et al. (2010)), the $B$ and $C$ states are degenerate and form the $j=3/2$ multiplet. However, since they have different electronic distributions, their degeneracy would be lifted by a tetragonal lattice distortion (see Appendix A). In terms of eigenstates $|j,m_{j}\rangle$ of $\mathbf{J}^{2}$ and $J^{z}$, we can identify the states created by $B_{\sigma}^{\dagger}$ and $C_{\sigma}^{\dagger}$ as

[TABLE]

Alternatively, we can use two pseudospins $1/2$ to label the four states in the $j=3/2$ subspace Natori et al. (2016). The first pseudospin is denoted $s$ and is associated with the Kramers degeneracy:

[TABLE]

The second pseudospin, hereafter called pseudo-orbital $\tau$, is defined by

[TABLE]

In the notation of $|s,\tau\rangle$, with $s,\tau=\pm 1/2$, we write

[TABLE]

This definition is such that the $z$ component of the total angular momentum is given by

[TABLE]

where the operators $s^{z}$ and $\tau^{z}$ are defined by

[TABLE]

More generally, if we define the vector of annihilation operators $\xi\equiv(C_{\uparrow},C_{\downarrow},B_{\uparrow},B_{\downarrow})^{t}$, we have a basis of operators in the space of a singly occupied $j=3/2$ level:

[TABLE]

where $\sigma^{a}$, with $a\in\{x,y,z\}$, are Pauli matrices and $\mathbb{I}$ is the $2\times 2$ identity matrix. The pseudospin and pseudo-orbital operators obey the SU(2) algebra:

[TABLE]

Similar spin and orbital operators appear in Kugel-Khomskii models for $e_{g}$ orbitals without SOC Feiner et al. (1997); Khaliullin and Oudovenko (1997); Li et al. (1998); Joshi et al. (1999); Wang and Vishwanath (2009); Corboz et al. (2012); Lajkó and Penc (2013). A crucial difference is that here $\mathbf{s}$ and $\boldsymbol{\tau}$ act on $j=3/2$ states, in which spin and orbital degrees of freedom are entangled as described by Eqs. (5).

Table 1 shows how the multipoles of J are written in terms of the components of s and $\boldsymbol{\tau}$. Here we introduce the linear combinations of $\tau^{x}$ and $\tau^{z}$:

[TABLE]

According to Table 1, s is a linear combination of dipole and octupole moments of $\mathbf{J}$ in the $\Gamma_{4}$ representation. Similarly, $\tau^{x}$ and $\tau^{z}$ correspond to quadrupoles in the $\Gamma_{3}$ representation. The component $\tau^{y}$ appears separately as a one-dimensional representation $\Gamma_{2}$. We conclude that all the components of s are odd under conjugation by the time-reversal operator $T$. As for the pseudo-orbital $\boldsymbol{\tau}$, the $\tau^{x}$ and $\tau^{z}$ components are even while $\tau^{y}$ is odd under time reversal. More explicitly,

[TABLE]

II.2 Interacting spin model

In this subsection, we reproduce the derivation of the effective spin model for $d^{1}$ double perovskites following Ref. Romhányi et al. (2017). We present this derivation here for completeness and to mention some important aspects in the interpretation of the model parameters.

We start from the multi-orbital Hubbard model

[TABLE]

Here $i$ labels the lattice sites, $\alpha\beta$ labels the $t_{2g}$ orbitals, $\sigma$ is the electronic spin, $U$ is the Coulomb interaction, and $J_{H}$ is Hund’s coupling. We use the ordering convention $xy<yz<zx$. The hopping processes are restricted to nearest-neighbor sites, as in Ref. Chen et al. (2010), such that $\gamma$ labels the axis perpendicular to the $\alpha\beta$ plane of the $\langle i,j\rangle$ bond (see Fig. 1b).

Using the single-occupancy constraint $\underset{\alpha\beta}{\sum}n_{i,\alpha\beta}=1$ and applying perturbation theory in the regime $t\ll U,J_{H}$, we obtain the spin-orbital model Romhányi et al. (2017)

[TABLE]

in which $P_{ij}^{(\gamma)}=n_{i,\alpha\beta}\bar{n}_{j,\alpha\beta}+\bar{n}_{i,\alpha\beta}n_{j,\alpha\beta}$ with $\bar{n}_{i,\alpha\beta}=n_{i,\beta\gamma}+n_{i,\gamma\alpha}=1-n_{i,\alpha\beta}$. The coupling constants $J$, $J^{\prime}$ and $V$ are given by

[TABLE]

where $K=4t^{2}/U$, $r_{1}=1/(1-3\eta)$, $r_{2}=1/(1-\eta)$ and $r_{3}=1/(1+2\eta)$, with $\eta=J_{H}/U$.

Equation (22) can be compared with the model derived in Ref. Chen et al. (2010). The first line corresponds to the antiferromagnetic exchange interaction, with a correction in the sign of the spin-independent term. The second line is formally the ferromagnetic Hamiltonian; the difference lies in the interpretation of the parameter $J^{\prime}$. We find that $J^{\prime}$ is related with the ratio $J_{H}/U$ of the transition metal, instead of the ratio at the oxygen site. The third line is similar to the electric quadrupole interaction discussed in Ref. Chen et al. (2010), differing by the absence of a term proportional to $(n_{i,\beta\gamma}-n_{i,\gamma\alpha})(n_{j,\beta\gamma}-n_{j,\gamma\alpha})$. In summary, our minimal model also contains antiferromagnetic, quadrupole and ferromagnetic interactions. However, the explanation of the model parameters comes from a different mechanism.

Hereafter we focus on the limit of vanishing Hund’s coupling $\eta=0$, in which $J=K$ and $J^{\prime}=V=0$. This corresponds to the regime in which quantum fluctuations are maximized and may favor a QSL ground state Chen et al. (2010). The final step is the projection of $H_{\text{so}}$ onto the $j=3/2$ manifold in the limit $\lambda\gg K$:

[TABLE]

Using the pseudospins and pseudo-orbitals discussed in Subsection II.1, we find for $\eta=0$ Natori et al. (2016); Romhányi et al. (2017):

[TABLE]

The bond-dependent exchange processes are represented in Fig. 1b. In analogy with the effective models for iridates Jackeli and Khaliullin (2009); Chaloupka et al. (2010); Rau et al. (2014), the anisotropy arises from the directionality of the $t_{2g}$ orbitals Khaliullin (2005). In this notation, the hidden SU(2) symmetry of the effective model discussed in Ref. Chen et al. (2010) becomes transparent. More explicitly, if we define the total pseudospin operator

[TABLE]

then $[H_{\text{eff}},\textbf{s}_{\text{tot}}]=0$. This continuous symmetry, unexpected for general spin-orbit coupled systems, enhances quantum fluctuations and favors unconventional magnetic states Chen et al. (2010). In addition, a $Z_{3}$ symmetry corresponding to a $2\pi/3$ rotation of the $(\tau^{z},\tau^{x})$ vector and analogous to the symmetry of quantum compass models Nussinov and van den Brink (2015) is made evident by the $\tau^{\alpha\beta}$ pseudo-orbital operators. These symmetry properties of $H_{\text{eff}}$ play an important role in the ansatz for the chiral spin-orbital liquid to be discussed in Sec. III. The expression for the more general projected Hamiltonian with $\eta\neq 0$ is given in Appendix B.

III Mean-Field Theory of the Chiral Spin-Orbital Liquid

In the following we describe the parton mean-field theory that gives rise to the chiral spin-orbital liquid studied in Ref. Natori et al. (2016). The motivation for considering a Majorana fermion parton construction arises mainly from the similarities between the double perovskite model and the Kitaev model, as they both contain bond-dependent anisotropic exchange interactions. The main point of this section is to discuss the spectrum of Majorana fermion excitations, which will be important to interpret the response functions discussed in Secs. IV and V.

The operators s and $\boldsymbol{\tau}$ obeying the algebra in Eqs. (16) through (18) can be represented by Majorana fermions in the following way Tsvelik (1992); Coleman et al. (1994); Biswas et al. (2011); Herfurth et al. (2013); Wang and Vishwanath (2009)

[TABLE]

where $a=x,y,z=1,2,3$ for the Majorana fermion flavors. The six Majorana fermions $\zeta^{a}\in\{\eta^{a},\theta^{a}\}$ obey $(\zeta^{a})^{\dagger}=\zeta^{a}$ and $\left\{\zeta^{a},\zeta^{b}\right\}=2\delta^{ab}$. This representation has a $Z_{2}$ gauge structure because the sign of the fermions can be changed ($\eta^{a}\rightarrow-\eta^{a}$ and $\theta^{a}\rightarrow-\theta^{a}$) without modifying the local physical operators. Since the Hilbert space is enlarged, one needs to impose a local constraint at each site $j$ to identify the physical states:

[TABLE]

The local constraint also implies that

[TABLE]

Using Eq. (29), we rewrite the Hamiltonian in Eq. (27) as

[TABLE]

Here we have introduced the fermions $\theta^{\alpha\beta}$ and $\bar{\theta}^{\alpha\beta}$ as linear combinations of $\theta^{1}$ and $\theta^{3}$ in analogy with Eqs. (19):

[TABLE]

Let us analyze some symmetries of the Hamiltonian. In general, we can define a six-component column vector of Majorana fermions $\zeta=(\eta^{1},\dots,\theta^{3})^{t}$ that transform as

[TABLE]

where $R$ is an SO(6) matrix. Although the Hamiltonian in Eq. (32) is not invariant under global SO(6) transformations, it is invariant under a subset that includes global rotations of the form $R=R_{\eta}\oplus\mathit{I}_{\theta}$, where $R_{\eta}$ corresponds to an SO(3) rotation of the vector $\boldsymbol{\eta}=(\eta^{1},\eta^{2},\eta^{3})$, and $\mathit{I}_{\theta}$ is the identity matrix in the $\theta$ sector. This symmetry is nothing but the global SU(2) invariance of Hamiltonian (27) expressed in terms of Majorana fermions. Moreover, we can identify the $Z_{3}$ symmetry as being generated by the transformation $R=\mathit{I}_{\eta}\oplus M_{\theta}$, where $M_{\theta}$ is the $2\pi/3$ rotation matrix acting on the two-component vector $(\theta^{1},\theta^{3})$ leaving $\theta^{2}$ invariant.

The action of time reversal $T$ on Eq. (32) follows from the symmetry properties of s and $\boldsymbol{\tau}$ discussed in Subsection II.1. In the representation of Eq. (29), $T$ can be defined as complex conjugation supplemented by $T^{-1}\theta^{2}T=-\theta^{2}$, while leaving the other flavors invariant. With this rule, the Hamiltonian in Eq. (32) is explicitly time-reversal invariant.

We construct a mean-field theory with the expectation values of bond operators $\langle\zeta_{i}^{a}\zeta_{j}^{b}\rangle$ as order parameters. These parameters are chosen in a way that preserves as many symmetries of the Hamiltonian (32) as possible. Since the fcc lattice contains triangular plaquettes, a Majorana QSL necessarily breaks time reversal and reflection symmetries Fiete et al. (2012). As a result, the mean-field theory can preserve at most the SO(3), $Z_{3}$ and some point-group symmetries. To preserve the SO(3) symmetry, the state must remain invariant under any global rotation of s. Consequently, all order parameters of the type $\langle\eta_{i}^{a}\eta_{j}^{b}\rangle$ with $a\neq b$ vanish, and $\langle\eta_{i}^{a}\eta_{j}^{a}\rangle=\langle\eta_{i}^{b}\eta_{j}^{b}\rangle$, for $a,b=1,2,3$. Similarly, $\langle\eta_{i}^{a}\theta_{j}^{b}\rangle$ vanishes for any pair $(a,b)$. The $Z_{3}$ symmetry rotates the $\tau^{\alpha\beta}$ operators among themselves. Thus, requiring $Z_{3}$ invariance implies that $\langle\theta^{\alpha\beta}\theta^{2}\rangle$ must also be zero. Applying these restrictions, we perform the mean-field decoupling of Eq. (32) and obtain the mean-field Hamiltonian

[TABLE]

where $iu_{ij}=\langle\eta_{i}^{a}\eta_{j}^{a}\rangle$, $iv_{ij}=\langle\theta_{i}^{2}\theta_{j}^{2}\rangle$, $iw_{ij}^{\alpha\beta}=\langle\theta_{i}^{\alpha\beta}\theta_{j}^{\alpha\beta}\rangle$ and $i\bar{w}_{ij}^{\alpha\beta}=\langle\bar{\theta}_{i}^{\alpha\beta}\bar{\theta}_{j}^{\alpha\beta}\rangle$. Notice that all quadratic terms are diagonal in the flavor index except those involving $\theta^{1}$ and $\theta^{3}$. Moreover, the decoupled terms for $\eta^{a}$ and $\theta^{2}$ fermions differ only by the corresponding mean-field amplitudes.

Our study of Eq. (38) is also restricted to translationally invariant ansätze. In this regard, one must impose the magnitude of each order parameter to be uniform:

[TABLE]

where $\phi_{ij}=\pm 1$ are $Z_{2}$ link variables. The anticommutation relations of the Majorana fermions ensure that $\phi_{ij}=-\phi_{ji}$, which gives an orientation to the links between the sites. To orient the links, it is convenient to subdivide the fcc lattice into four cubic sublattices, labeled by $X=1,2,3,4$. As can be seen in Fig. 3b, there is a correspondence between the sublattice sites and the vertices of elementary tetrahedra. The orientation of $\phi_{ij}$ between sublattices can be represented on a plane as shown in Fig. 3a. With two possible values for each $\phi_{ij}$, there are in total $2^{6}=64$ different “hopping” configurations, which can be grouped into eight non-gauge-equivalent ansätze.

While the order parameters are not gauge invariant, physically distinct ansätze can be labeled by the gauge-invariant $Z_{2}$ fluxes through the elementary plaquettes. Choosing three nearest-neighbor sites $(i,j,k)$ in a fixed orientation, one can define the flux $\chi_{ijk}\equiv-i\phi_{ij}\phi_{jk}\phi_{ki}$. The latter is closely related to the scalar spin chirality of $S=1/2$ QSLs Baskaran (1989); Wen et al. (1989); Wen (2002, 2004). Using Eq. (29), we can write

[TABLE]

The operator in Eq. (43) is odd under time reversal and is analogous to the spin chirality order parameter. Using our mean-field decoupling, we obtain

[TABLE]

The physical state is determined by the $Z_{2}$ flux configuration on all plaquettes. We should note that the fluxes through the faces of any tetrahedron are not all independent. If the sites on any given face are oriented clockwise with respect to an outward normal vector, the four fluxes obey the relation $\prod_{r=1}^{4}\chi_{r}=1$, where $r$ labels the faces of the tetrahedron.

Time reversal plays an important role in choosing the mean-field theory, since it relates pairs of non-equivalent gauge configurations (see Fig. 3a). In terms of $Z_{2}$ fluxes, $T$ inverts $\chi_{ijk}$ of every elementary plaquette of the lattice. Although not related by gauge transformations, two gauge choices related by $T$ lead to degenerate mean-field ground states. Still guided by symmetry principles, we study here the most symmetric ansätze, which are characterized by the same $Z_{2}$ flux through all faces of a tetrahedron. The imposition of translation invariance implies that two tetrahedra sharing an edge have opposite $Z_{2}$ fluxes (see Fig. 3b). In other words, our ansatz is a staggered-flux Majorana QSL, where the staggering is between nearest-neighbor tetrahedra.

We define the parity transformation $P$ as a reflection by a symmetry plane of the fcc lattice. As can be seen in Fig. 3b, $P$ inverts all flux orientations in the mean-field ansatz. Since the Majorana QSL breaks both $P$ and $T$ symmetries, it is classified as a chiral spin(-orbital) liquid Wen et al. (1989). However, notice that the antiunitary operator $PT$ is still a symmetry, with $(PT)^{2}=+1$. Since we are dealing with system of fractionalized quasiparticles, crystalline symmetries must be studied by means of a projective symmetry group (PSG) analysis Wen (2004), which was discussed in Ref. Natori et al. (2016). Due to the breaking of $P$ and $T$, the point group symmetry of Hamiltonian (32) is reduced from $O_{h}\times Z_{2}$, with $Z_{2}$ corresponding to time reversal, to a group isomorphic to $O_{h}$.

After fixing the ansatz, we solve the mean-field Hamiltonian using the Fourier mode expansion

[TABLE]

where $X=1,2,3,4$ is the sublattice index and $N$ is the total number of sites in the fcc lattice. The positions of the sites in sublattice $X$ are given by

[TABLE]

with $\boldsymbol{\delta}_{1}=(0,0,0)$, $\boldsymbol{\delta}_{2}=(1/2,1/2,0)$, $\boldsymbol{\delta}_{3}=(0,1/2,1/2)$, and $\boldsymbol{\delta}_{4}=(1/2,0,1/2)$ in units where the lattice parameter is set to 1. The operators $\zeta_{\textbf{k}X}^{a}$ obey $(\zeta_{\textbf{k}X}^{a})^{\dagger}=\zeta_{-\textbf{k}X}^{a}$ and $\{\zeta_{\textbf{k}X}^{a},\zeta_{\textbf{k}^{\prime}X^{\prime}}^{a}\}=\delta_{\textbf{k},-\textbf{k}^{\prime}}\delta_{X,X^{\prime}}$. Thus, $\zeta_{\textbf{k}X}^{a}$ can then be treated as complex fermions with well-defined occupation numbers if we split the first Brillouin zone of the cubic lattice into two halves, which can be mapped into each other by inversion. Only one of these halves is taken into account and will be called $\frac{1}{2}\text{BZ}$. It is worth pointing out that the PSG analysis shows that the mean-field ansatz is invariant under translations on the fcc lattice Natori et al. (2016). This can be understood intuitively by noting that translations by $\boldsymbol{\delta}_{X}$ exchange the sublattices but do not change the signs of the gauge-invariant fluxes represented in Fig. 3b.

The mean-field Hamiltonian in Eq. (38) can be rewritten in the form

[TABLE]

where $\zeta_{\mathbf{k}}=(\zeta_{\mathbf{k},1},\zeta_{\mathbf{k},2},\zeta_{\mathbf{k},3},\zeta_{\mathbf{k},4})^{t}$ for $\zeta\in\{\eta^{a},\theta^{2}\}$ are four-component spinors, and $\Theta_{\textbf{k}}=(\theta_{\textbf{k}1}^{1},\theta_{\textbf{k}2}^{1},...,\theta_{\textbf{k}4}^{3})^{t}$ is an eight-component spinor. To find the ground state of Eq. (47), first we study the $4\times 4$ matrix $\mathcal{H}_{1}(\textbf{k})$, given by

[TABLE]

with

[TABLE]

and

[TABLE]

Let $U_{\textbf{k}}$ be the unitary matrix that diagonalizes $\mathcal{H}_{1}(\mathbf{k})$:

[TABLE]

where $\Lambda_{1}(\textbf{k})$ is diagonal. Since the matrices in Eq. (48) obey the Clifford algebra $\{\Sigma^{a},\Sigma^{b}\}=2\delta^{ab}$, the eigenvalues of $\mathcal{H}_{1}(\mathbf{k})$ are simply $\pm|\mathbf{h}(\mathbf{k})|$ and are doubly degenerate. This is a Kramers-type degeneracy that can be explained by point group symmetries, as discussed in Ref. Natori et al. (2016). The mean-field Hamiltonian is diagonal in the basis of operators $\tilde{\zeta}_{\textbf{k}\lambda}$ given by

[TABLE]

with $\lambda=1,\dots,4$ being the band index. The order parameters $u$ and $v$ are determined by self-consistent equations:

[TABLE]

where the sum over $\mathbf{k}$ is restricted to $\mathbf{k}\in\frac{1}{2}\text{BZ}$. At zero temperature, we can replace the average occupation of the single-particle states by

[TABLE]

where $\Theta(x)$ is the Heaviside step function and

[TABLE]

are the dispersion relations of the $\eta$ and $\theta^{2}$ fermions, respectively, with $C_{\lambda}=-1$ for $\lambda=1,2$ and $C_{\lambda}=+1$ for $\lambda=3,4$.

The expressions for $u$ and $v$ coincide at zero temperature, except possibly for a sign depending on the relative sign between the parameters $2u+\bar{w}$ and $w$. Without loss of generality, we fix $u>0$ (which corresponds to fixing the sign of the $T$-symmetry-breaking order parameter). As discussed in Ref. Natori et al. (2016), the two cases $v=u$ or $v=-u$ give rise to two different ansätze, with different expressions for the $8\times 8$ matrix $\mathcal{H}_{2}(\mathbf{k})$. In the remainder of this work, we will deal with the case $v=u$ com . In this case, self-consistency of the mean-field equations implies $2u+\bar{w}>0$ and $w>0$ .

Having fixed $\text{sgn}(uv)>0$, we find that the $8\times 8$ matrix $\mathcal{H}_{2}(\textbf{k})$ in Eq. (47) is given by

[TABLE]

where

[TABLE]

We denote by $V_{\mathbf{k}}$ the matrix that diagonalizes $\mathcal{H}_{2}(\mathbf{k})$:

[TABLE]

The order parameters $w$ and $\bar{w}$ can be calculated similarly to Eq. (54), using the components of $V_{\mathbf{k}}$ instead of $U_{\mathbf{k}}$.

Figure 4 shows the dispersion relation for the different flavors of Majorana fermions. The bands are particle-hole symmetric and doubly degenerate for all flavors. The dispersion relations of all bands are qualitatively similar, differing mainly by their bandwidths. The most remarkable feature is that the band structure displays nodal lines along the edges of the Brillouin zone, a consequence of the vanishing of $\mathbf{h}(\mathbf{k})$ when two components of $\mathbf{k}$ are equal to $\pi$ Natori et al. (2016). The energy increases linearly with the distance in momentum space from a generic point on a nodal line. The exception is the vertex point $R=(\pi,\pi,\pi)$, where the nodal lines cross and the dispersion becomes approximately quadratic but anisotropic. For $\mathbf{k}=(\pi,\pi,\pi)+\mathbf{q}$, with $|\mathbf{q}|\ll 1$, we obtain for all bands

[TABLE]

which is of the form $\epsilon_{\mathbf{k}\lambda}=q^{2}f_{\lambda}(\Omega)$, with $f_{\lambda}(\Omega)$ a function of the spherical angle coordinates of $\mathbf{q}$.

The single-particle states in the neighborhood of the $R$ point dominate the low-energy physics due to the quadratic dispersion. To see this, we can compute the corresponding contribution to the density of states

[TABLE]

Thus, we find $\rho_{\text{point}}(E)\propto\sqrt{E}$, a vanishing density of states characteristic of a pseudogap. The same analysis for the density of states around generic points on the nodal line parallel to the $k_{z}$ axis yields

[TABLE]

where $v_{\lambda}(k_{z})$ is the effective velocity of the linear dispersion around the nodal line and we cut off the integral at $|k_{z}-\pi|=\epsilon>0$ to exclude the contribution from the $R$ point [since $v_{\lambda}(k_{z}\to\pi)\to 0$]. Thus, the contribution from the nodal lines to the density of states is $\rho_{\text{line}}(E)\propto E$. This is the same result as for a Dirac point in two dimensions. The comparison of Eqs. (65) and (66) suggests that the low-temperature thermodynamics of the chiral spin-orbital liquid should be governed by the quadratic band touching point.

IV Specific heat, spin-lattice relaxation rate and dynamical spin structure

factor

In the absence of a “smoking-gun” signature of QSLs Savary and Balents (2017), a proper characterization of such states must combine the response to different perturbations. In this section, we calculate the response of our chiral spin-orbital liquid to three well-established probes: specific heat, nuclear magnetic resonance, and inelastic neutron scattering.

IV.1 Specific Heat

The specific heat $C_{V}$ of Ba2YMoO6 was measured by de Vrie et al. de Vries et al. (2010) and Aharen et al. Aharen et al. (2010). In both experiments, the magnetic contribution was obtained by subtracting off the data for the isostructural nonmagnetic compound Ba2YNbO6 from the total specific heat of Ba2YMoO6. The measurements agree about the presence of a broad peak around 50 K. However, the reported values of $C_{V}$ at the maximum are different: $7.5$ J/mol$\cdot$K in Ref. de Vries et al. (2010) versus $2.5$ J/mol$\cdot$K in Ref. Aharen et al. (2010). By integrating $C_{V}$ out to $T\approx 200$ K, de Vries et al. de Vries et al. (2010) found that the entropy released is close to $k_{B}\ln 4$, as expected for a $j=3/2$ system. In Ref. de Vries et al. (2013), the low-temperature behavior of $C_{V}$ was interpreted as evidence for a pseudogap in magnetic excitations. On the other hand, Aharen et al. Aharen et al. (2010) noted that the entropy lost below 50 K is lower than $k_{B}\ln 2$ and found an abrupt drop in the magnetic specific heat above 60 K. While it would be desirable to clarify the disagreement between these experiments, here we will focus on the common observation of a broad peak in $C_{V}$ and use this information to set the energy scale in our spin-orbital model.

To calculate $C_{V}$, we follow the method of Ref. Herfurth et al. (2013), which studied a Majorana QSL on a $S=1/2$ Heisenberg model. The mean-field theory described in Sec. III can be extended to $T>0$ by replacing the average occupation of single-particle states by the Fermi-Dirac distribution:

[TABLE]

where $\beta=1/(k_{B}T)$. The order parameters calculated using Eq. (67) define a temperature-dependent mean-field Hamiltonian $H_{\text{MF}}(T)$. We fix these parameters by minimizing the free energy,

[TABLE]

and solving the self-consistent equations numerically. The absolute values of the order parameters decrease with increasing temperature as shown in Fig. 5a.

We analyze the free energy for small values of the order parameters in Appendix C and show that they vanish at the critical temperature $k_{B}T_{c}=J/12$. The parameters $u$ and $\bar{w}$ vanish as $(T_{c}-T)^{1/2}$, as expected for primary order parameters at the mean-field level. Note that this is a well-defined second-order phase transition because a nonzero value of $u$ implies spontaneous breaking of time reversal symmetry [see Eq. (44)]. On the other hand, $v$ and $w$ behave as secondary order parameters Stokes and Hatch (1991); Hatch and Stokes (2001) and vanish as $(T_{c}-T)^{3/2}$ (see Fig. 5a).

At low temperatures $T\ll J$, we can approximate the order parameters by their zero-temperature values. The main effect of thermal fluctuations in this regime is to change the occupation of the states in a band with fixed bandwidth. Using the density of states in Eq. (65), we find

[TABLE]

Thus, at sufficiently low temperatures we obtain the power-law behavior $C_{V}\propto T^{3/2}$.

The whole temperature dependence of the specific heat is shown in Fig. 5b. Starting from the low-temperature limit, we see that the $T^{3/2}$ behavior turns into a small plateau at $k_{B}T\approx 0.02J$. Above this temperature, there is a regime where $C_{V}$ increases approximately linear with $T$, followed by a sharp drop at $k_{B}T_{c}=J/12$ (which is a discontinuity at the mean-field level). Our theoretical result shows qualitative agreement with the experimental data obtained in Ref. Aharen et al. (2010). To make some quantitative predictions, we use the experimental data from Ref. Aharen et al. (2010) to estimate $T_{c}\approx 70$ K. This fixes the exchange coupling constant at $J\approx 72\text{ meV}$.

The lost entropy per site calculated within the parton mean-field theory is approximately $1.98k_{B}$, significantly higher than the expected for a $j=3/2$ system ($k_{B}\ln 4\approx 1.39k_{B}$). We expect the mean-field result to overestimate the entropy since this approximation violates the local constraint in Eq. (30). As a result, the number of microstates in this approach is higher than the actual number of physical states. At zero temperature, this problem was circumvented by using the Variational Monte Carlo (VMC) method to find a better estimate of the ground state energy Natori et al. (2016). To our knowledge, the only calculations of thermodynamic quantities in controlled approximations of QSLs at finite temperatures were done recently for the Kitaev model Nasu et al. (2014, 2015); Nasu and Motome (2016); Yoshitake et al. (2016). The numerical methods benefit from the exact solvability of the Kitaev model, a feature not available in our case.

Another problem with the mean-field approach used in this section is that it implicitly assumes that the Z2 gauge configuration in the mean-field ansatz remains frozen at finite temperatures. Without this assumption, we would not be able to diagonalize a free-fermion Hamiltonian and find the dispersion relations used in Eq. (67). Remarkably, studies of thermodynamics of the Kitaev model Nasu et al. (2014, 2015); Nasu and Motome (2016); Yoshitake et al. (2016) found that thermal fluctuations of the $Z_{2}$ gauge field are activated at temperatures much lower than the bandwidth of the Majorana fermions at zero temperature. The proliferation of thermally excited visons is detected as an additional peak in the specific heat. In the case of three-dimensional QSLs Nasu et al. (2014), the lower-temperature peak in $C_{V}$ is a true singularity and signals a topological phase transition predicted by $Z_{2}$ gauge theory Polyakov (1978); Senthil and Fisher (2000).

Nevertheless, we still argue in favor of using the broad peak at higher temperature to determine the energy scale of the exchange interactions. We note that in controlled numerical calculations for Kitaev models the corresponding peak in $C_{V}$ is well described by the approximations of either fixing a uniform configuration or treating the $Z_{2}$ gauge field as a completely random variable Nasu et al. (2015). In the following we will use the estimate $J\approx 72\text{ meV}$ to analyze the energy scales that appear in INS and RIXS.

IV.2 Spin-lattice relaxation rate

Nuclear magnetic resonance is a technique that relies on nuclear spins to probe the local environment. In spin systems, the energy transfer between electrons and nuclei is mediated by the hyperfine coupling

[TABLE]

where $\mathbf{I}_{i}$ is the nuclear spin at site $i$ and $\textbf{B}_{\text{hf}}(i)$ is the hyperfine effective field due to neighboring electrons. In the experiment of Ref. Aharen et al. (2010), the excited nuclear spin was the $I=1/2$ 89Y, which couples to the $j=3/2$ magnetic moment of Mo electronic spins. We can then write

[TABLE]

where $\boldsymbol{\delta}$ is the relative position of the atoms of 89Y and their neighboring Mo atoms, and $A_{0}$ is the constant hyperfine coupling for first-neighbor $\boldsymbol{\delta}$. The spin-lattice relaxation rate $1/T_{1}$ is given by

[TABLE]

where $\omega$ is the resonance frequency,

[TABLE]

is the hyperfine interaction form factor, and $\chi_{+-}^{\prime\prime}(\textbf{q},\omega)$ is the spectral function given by

[TABLE]

with $|n\rangle$ being an exact eigenstate of the spin Hamiltonian with energy $E_{n}$, and $Z=\sum_{n}e^{-\beta E_{n}}$ being the partition function. Here, $J_{j}^{-}=J_{j}^{x}-iJ_{j}^{y}$ is the angular momentum lowering operator at site $j$ and $J_{\mathbf{q}}^{-}$ is its Fourier transform.

We calculate $1/T_{1}$ for the chiral spin-orbital liquid using the parton mean-field theory. The main idea is to write $J_{j}^{-}$ in terms of $\mathbf{s}$ and $\boldsymbol{\tau}$ according to Table 1 and relate the spectral function $\chi^{\prime\prime}_{+-}(\mathbf{q},\omega)$ to finite-temperature correlations of free Majorana fermions. In this approach, we employ the order parameters calculated self-consistently at finite temperature as described in Subsection IV.1.

To gain some insight into the low-temperature behavior of $1/T_{1}$, we find it instructive to first analyze the contribution of the $\eta$ fermions to the total spectral function, since in this case we can derive some closed-form expressions. Using the procedure outlined in Appendix D, we find that the $\eta$-fermion contribution in the experimentally relevant regime $\omega\ll k_{B}T$ is given by

[TABLE]

where

[TABLE]

At low temperatures $k_{B}T\ll J$, the spin-lattice relaxation rate is dominated by excitations with small momentum transfer near the quadratic band touching point. We write $\textbf{k}=(\pi,\pi,\pi)+\textbf{q}$ and $\textbf{k}^{\prime}=(\pi,\pi,\pi)+\textbf{q}^{\prime}$, with $|\mathbf{q}|,|\mathbf{q}^{\prime}|\ll 1$. In this case, the energies can be approximated by Eq. (64) and the vector $\textbf{h}(\textbf{k})$ by

[TABLE]

where $\Omega$ is the solid angle in spherical coordinates. We can also approximate $A(\textbf{k}-\textbf{k}^{\prime})\approx 3$ and

[TABLE]

Gathering all these approximations, we verify that $(1/T_{1})_{\eta}$ scales as $T^{2}$ for $T\to 0$, as could be anticipated from the low-energy density of states in Eq. (65). A similar calculation assuming momenta near the Fermi lines leads to $(1/T_{1})_{\eta}\propto T^{3}$. While this result refers to the contribution from $\eta$ fermions, it also reflects the qualitative behavior of the total $1/T_{1}$ since the dispersion relation of the $\theta$ fermions is qualitatively similar.

We have calculated the total spectral function $\chi_{+-}^{\prime\prime}(\textbf{q},\omega)$ numerically, including the contribution from $\theta$ fermions and at arbitrary temperatures, as explained in Appendix D. The result for the spin-lattice relaxation rate is shown in Fig. 6. At low temperatures, the behavior is dominated by the $R$ point and is described by the power law $1/T_{1}\propto T^{2}$ discussed above. An abrupt increase of $1/T_{1}$ is verified near the critical temperature $T_{c}$, followed by a constant behavior at higher temperatures. This result should be compared with figure 15(b) of Ref. Aharen et al. (2010). While the suppression of $1/T_{1}$ at low temperatures was interpreted as evidence for a gapped QSL, the experimental result is also qualitatively consistent with a pseudogap in the low-energy density of states. This makes the chiral spin-orbital liquid state a valid alternative to explain the spin-lattice relaxation rate of Ba2YMoO6.

IV.3 Inelastic Neutron Scattering

Neutron scattering is the standard probe to study magnetic ordering and excitations in condensed matter. At zero temperature, the magnetic scattering cross section for polarized neutrons is proportional to one component of the dynamical structure factor

[TABLE]

where $|g\rangle$ is the ground state, $\mathbf{q}$ and $\omega>0$ are the momentum and energy transferred by the neutron, and $|n\rangle$ is an excited state of the many-body Hamiltonian.

Here we will calculate the dynamical structure factor for the chiral spin-orbital liquid. It follows from $PT$ and point group symmetries that $S^{ab}(\mathbf{q},\omega)\propto\delta^{ab}$ and $S^{aa}(\mathbf{q},\omega)=S^{aa}(-\mathbf{q},\omega)$. Writing the operator J in terms of s and $\boldsymbol{\tau}$ and then Majorana fermions, we obtain

[TABLE]

The form factor for $S^{zz}(\mathbf{q},\omega)$ is

[TABLE]

where $\mathbf{G}$ is a vector of the cubic reciprocal lattice chosen such that $\mathbf{q}-\mathbf{k}+\mathbf{G}$ is contained in the first Brillouin zone. The components $S^{xx}(\mathbf{q},\omega)$ and $S^{yy}(\mathbf{q},\omega)$ can be obtained from Eq. (81) by cyclic permutation of all indices $a=1,2,3=x,y,z$. Within the mean-field theory, the excited states are restricted to two-particle excitations. The form factors can be calculated using the matrix elements of $U_{\mathbf{k}}$ and $V_{\mathbf{k}}$ defined in Eqs. (51) and (63).

Figure 7b shows the dynamical structure factor $S^{xx}(\mathbf{q},\omega)$ along the high-symmetry lines of the Brillouin zone of the fcc lattice (see Fig. 7a). As expected for QSLs, the spectral weight is distributed over a continuum of fractionalized excitations. The maximum intensity is found at the $L$ point, corresponding to momentum transfer $\mathbf{q}=(\pi,\pi,\pi)$. The energy scale at the maximum is of the order of the bandwidth of the Majorana fermions $\eta^{a}$ and $\theta^{1,3}$ shown in Fig. 4b. Using $J\approx 72$ meV estimated from the specific heat, we find that the peak in the dynamical structure factor appears at $\omega\approx 0.25J\approx 18\text{ meV}$.

We now compare our theoretical results with the neutron scattering experiments reported by Carlo et al. Carlo et al. (2011) done in polycrystalline samples. To make the comparison, we average the dynamical structure factor over momenta with absolute value $Q$ in the range 1.5Å${}^{-1}<Q<1.8$Å*-1*. This range includes the point equivalent to $L$ called $L^{\prime}=a^{-1}(3\pi,3\pi,-\pi)$, at which $Q=1.63$Å*-1* if we use the lattice spacing $a=8.389$Å de Vries et al. (2010). Our result in Fig. 7c shows a single broad maximum at $\omega\approx 18$ meV. By contrast, the experimental result shows a three-peak structure, with a pronounced magnetic peak at $\omega\approx 28$ meV and two smaller ones at 11 and 17 meV. As noted by the authors of Ref. Carlo et al. (2011) , the energy scale of the broad peak is a factor of 2 larger than the one inferred from the spin-lattice relaxation rate.

Our result for the chiral spin-orbital liquid at the mean-field level does not predict such a three-peak structure. Our model does contain multiple energy scales associated with the nondegenerate Majorana fermion bands shown in Fig. 4b, but the bandwidths of the flavors $\eta^{a},\theta^{1,3}$, which appear in the form factor Eq. (81), are rather close to each other. We also recall that our calculations were done neglecting fluctuations of the $Z_{2}$ gauge field and interactions between Majorana fermions. The inclusion of these effects in a bosonic spin liquid on the kagome lattice Punk et al. (2014) led to broadening and shift of the spectral weight of $S(\textbf{q},\omega)$ when compared with the mean-field theory. We expect a similar broadening in our case if gauge fluctuations are taken into account.

V RIXS cross sections

RIXS is a photon-in photon-out spectroscopic technique that probes excitations in solid state systems by measuring the energy, momentum and polarization changes of the scattered photon Ament et al. (2011). It is a resonant technique because the x-ray is tuned to coincide with the atomic transition between a core and a valence level of a given atom. The resonance turns an otherwise negligible second-order perturbation into the dominant contribution to the scattering amplitude. Moreover, the transitions involved in the absorption and emission processes are more complex than the ones generated by the probes listed in Sec. IV, allowing for the experimental study of a manifold of elementary excitations.

In this section, we evaluate and analyze the RIXS scattering operators for cubic double perovskites and calculate the RIXS cross sections for the chiral spin-orbital liquid. In subsection V.1 we describe the RIXS processes as well as the approximations used in our calculation, and present a symmetry analysis of the scattering operators using the method described in Refs. Ament and Khaliullin (2010); Haverkort (2010); Savary and Senthil (2015). We stress that here the symmetry arguments are applied to $j=3/2$ operators in the strong SOC limit and $O_{h}$ point-group symmetry. This is in contrast to Ref. Savary and Senthil (2015), which focused on SU(2)-invariant spin-$1/2$ systems with negligible SOC. In Subsection V.2, we determine the operators that appear specifically in the scattering amplitudes for the $L$ edge. The RIXS cross sections for the chiral spin-orbital liquid are then calculated and analyzed.

V.1 Derivation and symmetry analysis of RIXS scattering operators

Consider a general $N$ electron system, described by a many-body Hamiltonian $H_{0}$. The total Hamiltonian describing the system is $H=H_{0}+H^{\prime}$, where $H^{\prime}$ describes the interaction between electrons and photons

[TABLE]

Concerning the electrons, $e$ is the charge, $m$ is the mass, $\textbf{r}_{i}$, and $\textbf{p}_{i}$ and $\boldsymbol{\sigma}_{i}$ are, respectively, the position, momentum and spin of the $i$-th electron. The photon is represented by the electromagnetic vector potential $\textbf{A}(\mathbf{r})$. In second quantization, $\textbf{A}(\mathbf{r})$ is written as

[TABLE]

where $\mathcal{V}$ is the volume, $\epsilon_{0}$ is the vacuum permittivity, and $a_{\textbf{k},\boldsymbol{\varepsilon}}^{\dagger}$ is the creation operator for a photon with wave vector k, frequency $\omega_{\mathbf{k}}$, and polarization vector $\boldsymbol{\varepsilon}$.

Our aim is to evaluate the x-ray scattering cross sections after treating the photons as perturbations. Let the initial electron-photon state be $|G\rangle$ and a set of final states be $\{|F\rangle\}$. Using Fermi’s golden rule to second order, we obtain the x-ray cross section

[TABLE]

in which $E_{\nu}$ and $1/\gamma_{\nu}$ are the energy and the lifetime of the intermediate state $|\nu\rangle$, respectively. We assume that the initial state corresponds to a direct product of a many-body electronic ground state $|g\rangle$ and an incident photon state: $|G\rangle=|g\rangle\otimes|\mathbf{k},\omega_{\mathbf{k}},\boldsymbol{\varepsilon}\rangle$. Similarly, the final state is a direct product of an excited electronic state $|n\rangle$ with energy $E_{n}$ and an emitted photon labeled by $|F\rangle=|n\rangle\otimes|\mathbf{k}^{\prime},\omega_{\mathbf{k}^{\prime}},\boldsymbol{\varepsilon}^{\prime}\rangle$. We also deal with the case where $\omega_{\textbf{k}}$ is tuned to the energy difference between an atomic core level and a valence shell state. The photon is totally absorbed and the $|\nu\rangle$ state contains an atomic core hole and an additional electron in the valence or conduction band. If the photon energy $\hbar\omega_{\mathbf{k}}$ is tuned so that $|E_{G}-E_{\nu}|\ll\gamma_{\nu}$, the system is said to be in resonance and the importance of second-order processes is maximized.

Four standard approximations will be used to evaluate the second-order terms in Eq. (84). First, we neglect the so-called “magnetic” contribution ($\propto\boldsymbol{\sigma}\cdot\nabla\times\textbf{A}$ ) of $H^{\prime}$. Second, we use the dipole approximation for the scattering amplitude and take $e^{i\mathbf{k}\cdot\mathbf{r}_{i}}\approx e^{i\mathbf{k}\cdot\mathbf{R}_{i}}$, where $\mathbf{R}_{i}$ represents the lattice point to which the $i$-th electron is bound. Third, we consider that the highly unstable core hole in the $|\nu\rangle$ state decays before it can hop to a different ion. Finally, we consider only direct RIXS processes, i.e., we neglect effects of the core-hole Coulomb potential on the valence electrons. Within this fast collision approximation Ament et al. (2011), RIXS probes only single-site operators. The cross section then depends only on $\mathbf{q}=\mathbf{k}-\mathbf{k}^{\prime}$ and $\omega=\omega_{\mathbf{k}}-\omega_{\mathbf{k}^{\prime}}$, which are, respectively, the momentum and energy transferred to the sample. Equation (84) can then be recast in the form

[TABLE]

where $\hat{\mathcal{O}}_{\mathbf{q}}$ is the so-called scattering operator in momentum space. The latter is obtained from the Fourier transform $\hat{\mathcal{O}}_{\mathbf{q}}=\underset{{i}}{\sum}e^{i\mathbf{q}\cdot\mathbf{R}_{i}}\hat{\mathcal{O}}_{i}$, where

[TABLE]

Here, the dipole operator

[TABLE]

acts on the electronic states bound to position $\mathbf{R}_{i}$.

Equations (86) and (87) show that the RIXS cross section (85) depends on the initial and final polarizations $\boldsymbol{\varepsilon}$ and $\boldsymbol{\varepsilon}^{\prime}$ and on the matrix elements of the electron position operator $\langle\nu|\mathbf{r}|g\rangle$ and $\langle n|\mathbf{r}|\nu\rangle$. The photon polarizations can be controlled in experiments (at least in principle). However, the matrix elements depend on details of the intermediate states for a particular compound. The general claim one can make is that, provided the final states are low-energy excitations, the scattering operators can be rewritten in terms of charge, spin and orbital degrees of freedom of the valence electrons. For magnetic insulators, RIXS operators correspond to a combination of spin and orbital angular momentum. This feature makes RIXS an attractive technique to investigate magnetic insulators with strong SOC, in which spins and orbitals cannot be treated as separate degrees of freedom.

Since $\hat{\mathcal{O}}_{i}$ is in general a complicated operator, it is desirable to start our RIXS analysis by determining: (i) which polarization vectors $\boldsymbol{\varepsilon}$ and $\boldsymbol{\varepsilon}^{\prime}$ we should choose to acquire the signatures of a given state; and (ii) which spin operators couple with these polarizations. The two issues can be tackled at once by an elementary symmetry analysis of Eq. (86). As the absorption and emission processes occur at the same ion, the operator $\hat{\mathcal{O}}_{i}$ should be invariant under operations of the point group symmetry of the site $\mathbf{R}_{i}$. In general, one starts by decomposing the scattering operator into irreducible representations of the point group, $\Gamma=\Gamma_{1}\oplus...\oplus\Gamma_{n}$. A basis for these representations is then constructed in terms of the polarization factors $\varepsilon^{\Gamma_{j}}$ and (pseudo)spins $\mathcal{J}^{\Gamma_{j}}$, in the form Ament and Khaliullin (2010); Haverkort (2010); Savary and Senthil (2015)

[TABLE]

where $\text{mul}(\Gamma_{j})$ is the multiplicity of the irreducible representation $\Gamma_{j}$, the dot represents a symmetric contraction of all indices, and $\kappa_{\Gamma_{j},l_{j}}$ are material specific coefficients.

The bases of the irreducible representations of the octahedral group in terms of multipoles of $j=3/2$ moments are known Santini et al. (2009); Shiina et al. (1998) and are reproduced in Table 1, together with their representation in terms of s and $\boldsymbol{\tau}$ pseudospins. It is also easy to verify that the following polarization factors form the bases $\varepsilon^{\Gamma_{j},l_{j}}$:

[TABLE]

Here, the vector P corresponds to the $\Gamma_{4}$ representation, T to $\Gamma_{5}$, $Q_{2}$ and $Q_{3}$ to $\Gamma_{3}$, and $U$ to the scalar representation. Combining the operators in the same irreducible representation in Table 1 with the polarization factors in Eqs. (89) according to Eq. (88), we find the general form of all transition operators [except for the scalar representation, which couples with the Casimir operator $\textbf{J}^{2}=j(j+1)=\text{const.}$]. We then see that RIXS can in principle directly probe pseudospin and pseudo-orbital excitations.

V.2 L-edge RIXS cross section of Mo5+

We now focus on the $L$-edge RIXS operators for $4d^{1}$ and $5d^{1}$ orbital systems retaining cubic symmetry, whose mechanism is illustrated in Fig. 8. At the $L$ edge Ament et al. (2011), 2$p$ core electrons are excited to the $B$ and $C$ states of Eq. (5). To describe the core-hole states, we first note that they are similar to the $t_{2g}$ valence states, since they result from the combination of spin-$1/2$ states with orbital angular momentum $L=1$ in the presence of SOC. The core-hole Hamiltonian is

[TABLE]

where $\Lambda>0$ is the SOC constant for the $2p$ states. Like in the $d^{1}$ valence electron, there is a lifting of the six-fold degeneracy into a doublet and a quadruplet. However, now the doublet has lower energy (see Fig. 8). We refer to the excited hole in the $j=1/2$ ($j=3/2$) multiplet as the resonant $L_{2}$ ($L_{3}$) edge.

Most of the dipole matrix elements of Eq. (87) vanish by symmetry. The remaining terms are written in second quantization as Kim and Khaliullin (2017)

[TABLE]

in which we have dropped a multiplicative factor $\langle 4d_{yz}|y|2p_{z}\rangle$. We can simplify Eq. (86) by writing $\gamma_{\nu}\approx\gamma_{\mu}=\text{const.}$, with $\mu=2,3$, for all intermediate states in the $L_{\mu}$ edge. Here, $\gamma_{\mu}$ is the average decay rate of the intermediate core-hole states. This approximation, together with the ones discussed in Section V.1, leads to

[TABLE]

where $\mathcal{P}_{\mu}$ is the projection operator of the intermediate states in the $L_{\mu}$ edge.

We can derive expressions for the scattering operator in terms of pseudospins $\mathbf{s}$ and $\mathbf{\tau}$ by taking the projection in the $j=3/2$ subspace and using the single-occupancy constraint $\sum_{\sigma}(B_{\sigma}^{\dagger}B_{\sigma}+C_{\sigma}^{\dagger}C_{\sigma})=1$. In Appendix E, we provide a general expression for Eq. (92) including the effects of a tetragonal distortion that lifts the degeneracy between $B$ and $C$ states. Here, we restrict the discussion to the cubic limit. For the $L_{2}$ edge, we find

[TABLE]

where

[TABLE]

On the other hand, the scattering operator for the $L_{3}$ edge involves only the pseudospin $\mathbf{s}$:

[TABLE]

After calculating the scattering operators $\hat{\mathcal{O}}_{i}^{L_{\mu}}$, the RIXS cross section in Eq. (85) can be calculated like the INS dynamical structure factor discussed in Section IV.3. The results for some representative operators are shown in Fig. 9. Once again, we find that the spectral weight is distributed over a broad continuum. A common feature for all these results is a maximum of intensity for transferred momentum at the $L$ point, $\textbf{q}=(\pi,\pi,\pi)$.

Interestingly, the $\theta^{2}$ fermion is excited in the cross section of the $L_{2}$ edge through the operators $\tau^{x}$, $\tau^{z}$ and $\textbf{s}\tau^{y}$, in sharp contrast with the dynamical structure factor for INS. Due to the reduced bandwidth of the $\theta^{2}$ fermions, the spectrum probed by RIXS (with the proper polarization) displays a narrower energy range when compared to the one measured by INS. This feature is readily verified when comparing Figs. 9b, 9c and 9d with Fig. 7b.

Let us turn to the $L_{3}$ edge, which detects pseudospin excitations directly. We can simplify the result by choosing $\boldsymbol{\varepsilon}$ and $\boldsymbol{\varepsilon}^{\prime}$ such that $P_{x}=P_{y}=0$, but $P_{z}\neq 0$. The cross section in this case is given by

[TABLE]

where $s_{\textbf{q}}^{z}$ is the Fourier transform of $s_{j}^{z}$. At the special point $\textbf{q}=0$, the form factor involves the conserved quantity $\textbf{s}_{\textbf{q}=0}=\textbf{s}_{\text{tot}}$ [see Eq. (28)], which commutes with the spin Hamiltonian. Since the ground state is a singlet of the pseudospin SU(2) symmetry, we have $\textbf{s}_{\textbf{q}=0}|g\rangle=0$. Thus, it follows from Eq. (96) that

[TABLE]

for any transferred energy $\omega$. This feature is clearly seen in Fig. 9a, and should be contrasted with the dynamical structure factor $S(\mathbf{q}=0,\omega)\neq 0$ for INS in Fig. 7b. This result is explicitly confirmed by the computation of the form factor in Eq. (96). At the mean field level, the excited state $\left|n\right\rangle$ is a two-particle excitation, in which the particles are characterized by well-defined momenta k and $\textbf{k}^{\prime}$. We can write $\left|n\right\rangle=\left|n(\mathbf{k},\mathbf{k}^{\prime})\right\rangle$, in which the vector $\mathbf{k}^{\prime}$ can take the values $\pm\mathbf{k}\pm\mathbf{q}$ according to the type of two-particle excitation under consideration. As shown in Appendix D (see Eq. (123)), the form factor in this case is

[TABLE]

This form factor clearly vanishes for $\mathbf{q}=0$. Therefore, this RIXS cross section could be used to detect the hidden SU(2) symmetry of the spin-orbital model for double perovskites.

We note that the dynamic structure factor for the operator $\textbf{s}\tau^{y}$ calculated at mean-field level also vanishes at $\mathbf{q}=0$ (see Fig. 9b). The reason is that the Majorana representation $s^{a}\tau^{y}=-i\eta^{a}\theta^{2}/4$, involves only $\theta^{2}$ and $\eta$ fermions, whose mean-field Hamiltonian is diagonalized by the same unitary transformation $U_{\mathbf{k}}$ given by Eq. (51). As a result, the form factor associated with $\textbf{s}\tau^{y}$ is also given by Eq. (98). However, since $\sum_{j}\textbf{s}_{j}\tau_{j}^{y}$ does not commute with the Hamiltonian, the vanishing of the spectral weight at $\mathbf{q}=0$ in this case is an artifact of the mean-field approximation.

We make here a final remark on the usefulness of RIXS to probe our QSL. Our discussion was restricted to one-site operators, but this technique can, in principle, probe operators involving two or three sites. Extending the symmetry arguments presented here, we predict that the chiral operator $\textbf{s}_{i}\cdot(\textbf{s}_{j}\times\textbf{s}_{k})$ could be probed and would couple with the polarization factor $\boldsymbol{\varepsilon}^{\prime\ast}\cdot\boldsymbol{\varepsilon}$ (see Table I of Ref. Savary and Senthil (2015)). Therefore, RIXS could in principle detect the $P$ and $T$ symmetry breaking of the chiral spin-orbital liquid in the elastic limit. Once again, we emphasize that our results were obtained within the parton mean-field theory. The role of $Z_{2}$ gauge fluctuations in the RIXS response deserves a separate and detailed study.

VI Conclusions

This paper presented a theoretical study of thermodynamic and spectroscopic properties of a $j=3/2$ Majorana chiral spin-orbital liquid. In the process, we fully developed a pseudospin representation of the $d^{1}$ orbital physics in Mott insulators with strong SOC. These results can guide the theoretical modeling, as well as the design and interpretation of experiments in compounds with similar local physics. Interestingly, the thermodynamic properties of the chiral spin-orbital liquid agree qualitatively with the available experimental results for the material Ba2YMoO6 de Vries et al. (2010, 2013); Carlo et al. (2011); Aharen et al. (2010). In particular, we find a sharp drop in the spin-lattice relaxation rate at low temperatures, even though the chiral spin-orbital liquid is a gapless phase. On the other hand, the inelastic neutron scattering cross section measured in Ref. Carlo et al. (2011) was not reproduced, since we found a single broad peak instead of the three-peak structure observed in polycrystalline samples. Adding effects beyond mean-field theory may explain this difference and will be left for future work.

As the main result of this paper, we showed that RIXS can selectively probe pseudospin and pseudo-orbital operators, and thus provide a direct way to detect quadrupolar and octupolar orders and excitations. Our results give some guidance to interpret RIXS spectra in 4$d^{1}$ and 5$d^{1}$ based compounds. In particular, we showed that the hidden SU(2) symmetry of the double perovskite model without Hund’s coupling can be demonstrated by probing pseudospin s excitations and observing the suppression of the spectral weight for momentum transfer at the $\Gamma$ point.

Finally, we note that the analysis of RIXS scattering operators studied here is also useful for other ordered double perovskites Chen et al. (2010); Romhányi et al. (2017); Svoboda et al. (2017). For magnetically ordered systems, the excitation spectrum can be fitted using a microscopic model [such as Eq. (106) in Appendix B] and representing spin-orbital excitations in terms of magnons within a spin-wave theory Li et al. (1998); Joshi et al. (1999). For instance, the onset of quadrupolar order in some osmium-based compounds observed in a recent study Lu et al. (2017) can be investigated in more detail using RIXS. Two recent RIXS measurements of spin waves, one in a compound preserving cubic structure Betto et al. (2017) and another in a $j=3/2$ compound Jeong et al. (2017), indicate that the theory developed in this paper can be tested in the near future.

We thank E. Andrade, F. A. Garcia, G. Jackeli, and E. Miranda for helpful discussions. This work was supported by Brazilian agencies FAPESP (W.M.H.N.) and CNPq (R.G.P.).

Appendix A Orbital physics with distortion

Throughout the main text, we kept our discussion of the orbital physics restricted to case of the cubic symmetry. In this appendix, we discuss the effects of tetragonal distortions on the $A,B,C$ states. The Hamiltonian (4) is redefined by

[TABLE]

where $\delta$ is the energy scale associated with the distortion. Notice that the time-reversal symmetry of $H_{\text{ion}}$ is preserved, which means that the eigenstates can still be organized into three Kramers pairs. In analogy with Eq. (5), we define

[TABLE]

in which the angle $\varphi$ is defined by

[TABLE]

The corresponding energies are given by

[TABLE]

showing how distortion lifts the degeneracy of the cubic limit. Notice that the s and $\boldsymbol{\tau}$ operators can still be used to describe the physics of the quadruplet formed by the orbitals $B$ and $C$.

Appendix B Effective Hamiltonian for nonzero Hund’s coupling

In the limit of strong SOC, we can project the Hamiltonian in Eq. (22) in the $j=3/2$ manifold as written in Eq. (26). Here we present the more general effective Hamiltonian for $\eta\neq 0$. We introduce the pseudo-orbital-dependent operators:

[TABLE]

in which $\bar{\tau}^{\alpha\beta}=\frac{1}{\sqrt{3}}(\tau^{\beta\gamma}-\tau^{\gamma\alpha})$. The Hamiltonian is given by

[TABLE]

The coupling constants $J$, $J^{\prime}$ and $V$ are defined by Eqs. (23), (24) and (25).

Appendix C Free energy near the critical point

In this appendix, we find an approximate expression for Eq. (68) near the finite-temperature critical point where the order parameters of the parton mean-field theory vanish. Expanding (68) up to the fourth-order in $\beta\epsilon_{\textbf{k}\lambda}$, we find

[TABLE]

where $\mathcal{K}=\beta J$. We reorganize $\Phi$ in the form

[TABLE]

where $\Phi_{2}$ contains the terms that are quadratic in the order parameters and $\Phi_{4}$ contains the quartic terms. The quadratic term can be written in matrix form $\Phi_{2}=\textbf{t}^{T}M\textbf{t}$, where $\textbf{t}^{T}=(u,v,w,\bar{w})$. Diagonalizing $M$, we find the set of eigenvalues $a_{n}$, $n=1,\dots,4$, given by

[TABLE]

The eigenvalues $a_{2}$ and $a_{4}$ vanish, respectively, at the temperatures $k_{B}T_{p}=J/36$ and $k_{B}T_{c}=J/12$. The critical temperature where the numerically calculated specific heat in Fig. 5 drops to zero corresponds to the higher value $T=T_{c}$.

Appendix D Computation of correlation functions

In this appendix, we outline the calculation of finite-temperature spectral functions such as the one in Eq. (74).

We start by considering the correlation function

[TABLE]

where $\hat{O}_{l}$ is a local operator acting on the $j=3/2$ subspace associated with site $l$, $\hat{O}_{l}(\tau)=e^{H_{\text{eff}}\tau}\hat{O}_{l}e^{-H_{\text{eff}}\tau}$ is the operator evolved in imaginary time, $T_{\tau}$ denotes time ordering, and $\langle\cdot\rangle=\text{Tr}(\rho\,\cdot)$ denotes the thermal average with density matrix $\rho=e^{-\beta H_{\text{eff}}}/Z$.

Quite generally, the local operator $\hat{O}_{l}$ can be written as a combination of Majorana fermion bilinears, $\zeta^{a}_{l}\zeta^{b}_{l}$, with $\zeta^{a}\in\{\eta^{a},\theta^{a}\}$. Let us illustrate the procedure by taking

[TABLE]

Within the mean-field approximation, the correlation function can be written as

[TABLE]

where

[TABLE]

is the noninteracting fermion Green’s function. If $\mathbf{R}_{l}$ belongs to the $X$ sublattice, $X=1,\dots,4$, and $\mathbf{R}_{m}$ to the $Y$ sublattice, we can write for $0<\tau<\beta$ [using momentum conservation and Eq. (52)]

[TABLE]

The thermal average yields

[TABLE]

Taking the Fourier transform of Eq. (113), we obtain

[TABLE]

where $\omega_{n}=(2n+1)\pi/\beta$, $n\in\mathbb{Z}$, are fermionic Matsubara frequencies.

Similarly, we obtain the Fourier transform of the correlation in Eq. (111):

[TABLE]

where $\omega_{m}=2\pi m/\beta$, $m\in\mathbb{Z}$, are bosonic Matsubara frequencies. The form factors are given by

[TABLE]

where $\mathbf{G}=2\pi(n_{x},n_{y},n_{z})$ with $n_{a}\in\mathbb{Z}$ are reciprocal lattice vectors chosen such that the momenta $\pm\mathbf{k}\pm\mathbf{q}+\mathbf{G}$ in each form factor lies in $\frac{1}{2}$BZ.

After an analytical continuation $i\omega_{m}\to\omega+i0^{+}$, we can take the imaginary part of the retarded correlation function $\chi^{\prime\prime}(\mathbf{q},\omega)$ in a standard way. In the regime $\beta\omega\ll 1$, we can approximate the factors of Fermi-Dirac distributions using $n_{F}(\epsilon+\omega)-n_{F}(\epsilon)\approx\omega dn_{F}/d\epsilon$. We then obtain the expression for the contribution from the $\eta$ fermions to the spin-lattice relaxation rate

[TABLE]

Eq. (118) can be further simplified since $|\epsilon_{\textbf{k}\lambda_{i}}^{(\eta)}|=|\epsilon_{\textbf{k}\lambda_{j}}^{(\eta)}|$ for $\lambda_{i},\lambda_{j}=1,...,4$ (see Eq. (57)). For $\mathcal{F}_{\lambda_{1}\lambda_{2}}^{(1)}(\textbf{k},\textbf{q})$, the sum over eigenstates yields

[TABLE]

where

[TABLE]

The expressions obtained for other sums differ from (119) only by the combination of vectors $\mathbf{k}$ and $\mathbf{q}$. Eq. (120) is the form factor $\mathscr{F}^{\eta}$ stated in Eq. (76). Notice that the reciprocal lattice vector $\mathbf{G}$ does not appear in this final expression.

From $\chi(\textbf{q},\omega)$ we can also recover the RIXS dynamical structure factor of $s^{z}$. Taking the zero temperature limit, we find

[TABLE]

in which

[TABLE]

Once again, summing over the eigenstates, we find for $\omega=0$

[TABLE]

with similar expressions for other summations. The expression in brackets is just the form factor written in Eq. (98). As stated in the main text, it is clear that this form factor will vanish when $\mathbf{q}=0$.

The procedure outlined in this appendix can be generalized for the $\theta$ fermions as well. In particular, if $\hat{O}_{l}=-i\eta_{l}^{a}\theta_{l}^{2}$, the corresponding form factor of a RIXS experiment will be exactly the one given in Eq. (123). For operators involving the fermions $\theta^{1}$ and $\theta^{3}$, it is not possible to find exact expressions to the form factors, since there is no closed form to the matrix $V_{\mathbf{k}}$ (see Eq. 63). The response functions must then be computed numerically.

Appendix E RIXS scattering operators

In this appendix, we present the RIXS scattering operators discussed in Section V considering an arbitrary tetragonal distortion. In general, we write

[TABLE]

where $\mu=2,3$ for the $L_{2,3}$ edge and we use the polarization factors given by Eq. (89). The above equation corresponds to the projection of the operators listed in Ref. Kim and Khaliullin (2017) to the $B$ and $C$ states discussed in Appendix A. In Table 2 we show the explicit values of the coefficients in terms of the angle parameter $\varphi$ in Eq. (101). We also highlight the coefficients in the cubic limit, which were expressed in Eqs. (93) and (95).