A Supplemental Material

# Chiral spin-orbital liquids with nodal lines

## Abstract

Strongly correlated materials with strong spin-orbit coupling hold promise for realizing topological phases with fractionalized excitations. Here we propose a chiral spin-orbital liquid as a stable phase of a realistic model for heavy-element double perovskites. This spin liquid state has Majorana fermion excitations with a gapless spectrum characterized by nodal lines along the edges of the Brillouin zone. We show that the nodal lines are topological defects of a non-Abelian Berry connection and that the system exhibits dispersing surface states. We discuss some experimental signatures of this state and compare them with properties of the spin liquid candidate BaYMoO.

###### pacs:
71.70.Ej, 75.10.Jm, 75.10.Kt

Quantum spin liquids (QSLs) are Mott insulators in which quantum fluctuations prevent long-range magnetic order down to zero temperature Balents (2010). They have received both experimental and theoretical attention due to predictions of unusual phenomena such as spin-gapped phases with topological order or gapless phases without spontaneous breaking of continuous symmetries Wen (2004). In recent years the evidence for QSLs in nature has started to look more auspicious thanks to the discovery of new compounds that realize the Heisenberg model on frustrated lattices Han et al. (2012). While frustration is a desirable ingredient, seminal work by Kitaev Kitaev (2006) has demonstrated that bond-dependent exchange interactions may provide another route towards QSL ground states. The key idea is that a spin- model on the (bipartite) honeycomb lattice with judiciously chosen anisotropic interactions can be rewritten in terms of free Majorana fermions hopping in the background of a static Z gauge field. The result is a QSL with exotic fractional excitations. The same idea has been applied to construct other exactly solvable models, including cases of higher spins Yao, Zhang, and Kivelson (2009); Mandal and Surendran (2009); Wu, Arovas, and Hung (2009); Hermanns and Trebst (2014); O’Brien, Hermanns, and Trebst (2016).

From a broader perspective, Kitaev’s model is an instance of a quantum compass model Dorier, Becca, and Mila (2005); Lacroix, Mendels, and Mila (2011); Nussinov and van den Brink (2015). Although Kitaev-type exactly solvable models are artificial, the kind of anisotropic interactions they presuppose arises naturally in Mott insulators with orbital degeneracy and strong spin-orbit coupling Jackeli and Khaliullin (2009); Chaloupka, Jackeli, and Khaliullin (2010). There is recent evidence that bond-dependent interactions are dominant in NaIrO Hwan Chun et al. (2015). While this compound is in a zigzag-ordered phase at low temperatures, the prospect of finding QSLs in compass models suggests inspecting other families of heavy-element transition metal oxides Pesin and Balents (2010); Wan et al. (2011); Witczak-Krempa et al. (2014).

All the conditions leading to quantum compass models can be found in Mott-insulating rock-salt-ordered double perovskites Chen, Pereira, and Balents (2010). Given the chemical formula ABBO, particularly interesting properties are found in compounds where the B magnetic ions have a 4d or 5d configuration. These ions are arranged in a face-centered-cubic (fcc) lattice, which, unlike the honeycomb lattice, is geometrically frustrated. The magnetic properties within this family are diverse Wiebe et al. (2002); Stitzer, Smith, and zur Loye (2002); Yamamura, Wakeshima, and Hinatsu (2006); Erickson et al. (2007), but the material that stands out is BaYMoO de Vries, Mclaughlin, and Bos (2010); Aharen et al. (2010); Carlo et al. (2011); de Vries et al. (2013). Despite a Curie-Weiss temperature K de Vries, Mclaughlin, and Bos (2010), several experiments point to the absence of long-range order down to K. Moreover, there is no sign of structural transitions, implying that the lattice remains cubic at low temperatures. Thus, the degeneracy of the orbitals is only partially lifted by the spin-orbit coupling, leading to a low-lying quadruplet de Vries, Mclaughlin, and Bos (2010). The effective model contains bond-dependent interactions between nearest-neighbor spins and is closely related to -matrix generalizations of Kitaev’s model Yao, Zhang, and Kivelson (2009); Wu, Arovas, and Hung (2009). Remarkably, the analysis in Chen, Pereira, and Balents (2010) revealed that, when antiferromagnetic exchange is dominant, ordered phases become unstable against quantum fluctuations, making this an interesting point to look for QSLs.

In this Letter we investigate a QSL in a realistic model for double perovskites with strong spin-orbit coupling. Using a representation of operators in terms of six Majorana fermions, we start by showing that a hidden SU(2) symmetry of the Hamiltonain becomes an explicit SO(3) symmetry for three of these fermions, whereas the other three exhibit a compass-model-type Z symmetry. As the model is not exactly solvable, we proceed with a mean-field approach and propose a spin liquid ansatz that preserves the SO(3) and Z symmetries. The ansatz breaks inversion and time reversal symmetry, thus describing a chiral spin liquid. Most interestingly, we find that the excitation spectrum is gapless along nodal lines which are vortices of a Berry connection in momentum space. This feature makes this chiral spin liquid a strongly correlated analogue of line-node semimetals and superconductors discussed in the context of topological phases of weakly interacting electrons Burkov, Hook, and Balents (2011); Chiu and Schnyder (2014); Yang, Pan, and Zhang (2014); Bian et al. (2016) and photonic crystals Lu et al. (2013). Going beyond the mean-field level, we use variational Monte Carlo (VMC) techniques Gros (1989); Ceperley, Chester, and Kalos (1977); Edegger et al. (2005) to show that our spin liquid state yields a remarkably low energy and should be regarded as a competitive candidate for the ground state of the spin-orbital model. Finally, we argue that the vanishing density of states at low energies predicted by our theory can account for some unusual properties observed in BaYMoO.

The spin-orbital model for cubic double perovskites with d electronic configuration is given by Chen, Pereira, and Balents (2010)

 H=J∑α,⟨ij⟩∈α(Si,α⋅Sj,α−14ni,αnj,α)−λ∑ili⋅Si. (1)

Here is the antiferromagnetic exchange between nearest-neighbor spins and is the atomic spin-orbit coupling. The index labels both planes (XY, YZ or XZ) and orbitals (, or ) Tokura and Nagaosa (2000). The operators and describe the number and the spin of electrons occupying the orbital on site , with the constraint , and is the effective orbital angular momentum of the states Ballhausen (1962). The total spin on site is .

In the regime , spin and orbital operators can be projected into the low-energy subspace of total angular momentum Chen, Pereira, and Balents (2010). The projected Hamiltonian , where is the projector, contains multipolar interactions in terms of . Our first step is to introduce operators and at each site as

 s=12(−Γ23,Γ13,Γ12),τ=12(Γ4,−Γ45,Γ5). (2)

The notation refers to five Dirac matrices given explicitly by

 Γ1 = σz⊗σy,Γ2=σz⊗σx,Γ3=σy⊗\mathbbm1, Γ4 = σx⊗\mathbbm1,Γ5=−Γ1Γ2Γ3Γ4, (3)

where , , are Pauli matrices, and 10 matrices Murakami, Nagaosa, and Zhang (2004); Yao, Zhang, and Kivelson (2009). The components of and satisfy the SU(2) algebra , , and . The relation between the basis of and the basis , with , is .

In the new representation the projected Hamiltonian assumes a relatively simple form:

 ~H=J9∑⟨i,j⟩∈α(%si⋅sj−14)(1−2ταi)(1−2ταj), (4)

where are given by , . A few comments are in order. First, Eq. (4) has the familiar form of a Kugel-Khomskii model Kugel and Khomski (1982); Feiner, Oleś, and Zaanen (1997). However, here the Kugel-Khomskii coupling involves pseudospins and pseudo-orbitals defined in the subspace, where the original spins and orbitals are highly entangled. Second, the Hamiltonian commutes with . This is a manifestation of the hidden global SU(2) symmetry discussed in Chen, Pereira, and Balents (2010). This continuous symmetry is unexpected, given that spin-orbit coupling breaks the conservation of , but appears in related models for orbitals Harris et al. (2003) and at special points of the Kitaev-Heisenberg model Chaloupka and Khaliullin (2015). Finally, the pseudo-orbital coupling has the form of a 120 compass model Nussinov and van den Brink (2015). There is a Z symmetry generated by followed by a C rotation of the planes.

In analogy with the spin liquid in Kitaev’s model Kitaev (2006), we now introduce a Majorana parton representation for the generators of SU(4) (i.e. the basis of operators). We write and operators as Tsvelik (2007); Coleman, Miranda, and Tsvelik (1994); Wang and Vishwanath (2009); Biswas et al. (2011); Herfurth, Streib, and Kopietz (2013); Hermanns, O’Brien, and Trebst (2015)

 saj=−i4ϵabcηbjηcj,τaj=−i4ϵabcθbjθcj. (5)

The components of and are Majorana fermions that obey , where . As the signs of the fermions can be changed (, ) without affecting the physical operators, this representation bears a Z redundancy. To eliminate the extra states, one can impose the local constraint Wang and Vishwanath (2009)

 Dj≡iη1jη2jη3jθ1jθ2jθ3j,Dj=1∀j. (6)

With this constraint we have , and Hamiltonian (4) becomes quartic in Majorana fermions:

 ~H = Missing or unrecognized delimiter for \left (7) +(η1iη2iη3j+η2iη3iη1j+η3iη1iη2j)¯θαj+(i↔j) +¯θαi¯θαjηi⋅ηj−θαiθαjθ2iθ2j],

where is the number of sites and and are defined by , , and , .

Hamiltonian (7) is invariant under global SO(3) rotations of the vector. The couplings involving the components of have only a discrete symmetry, namely the octahedral point group symmetry O of the lattice. The latter contains the Z that rotates and by 120 in the plane. In addition, is invariant under time reversal , where denotes complex conjugation. In terms of Majorana fermions, .

Next, we perform a mean-field decoupling of Hamiltonian (7). This is equivalent to neglecting fluctuations of the Z gauge field and yields qualitatively correct results as long as the system is in a QSL phase with deconfined Majorana fermions Tsvelik (2007). Our choice of mean field ansatz is guided by the condition of preserving the SO(3) and Z symmetries. This restricts the set of nonzero parameters . We obtain

 ~HMF = −NJ6+J36∑⟨j,l⟩∈α[i(2ujl+¯wαjl)ηj⋅ηl (8) +3iujl¯θαj¯θαl−iwαjlθ2jθ2l−ivjlθαjθαl +3u2jl+3¯wαjlujl−wαjlvjl],

where , , , and play the role of imaginary hopping amplitudes. Note that the symmetry implies decoupling of and fermions at the level of bilinear terms; yet, and remain coupled.

Seeking a translationally invariant ansatz, we set the order parameters to have uniform magnitude: , , , , with to be determined by self-consistent equations, whereas the ’s are chosen to be on each bond. Since e.g. , the choice of is equivalent to a choice of bond orientation and determines the gauge-invariant flux through elementary plaquettes. Noticing that the fcc lattice can be viewed as a network of edge-sharing tetrahedra, we obtain a symmetric ansatz by requiring that the Z fluxes, e.g. , be the same on all faces of a given tetrahedron, with sites on every triangle oriented counterclockwise with respect to an outward normal vector. This leads to the four-sublattice ansatz illustrated in Fig. 1.

Let us discuss the symmetry of our ansatz. First, we note that the Z gauge flux through triangles is odd under time reversal and is related to the spin chirality order parameter Baskaran (1989); Wen, Wilczek, and Zee (1989). The state also breaks inversion ; this can be seen from Fig. 1(b), since a mirror-plane reflection exchanges tetrahedra with opposite chiralities. Thus, our ansatz describes a chiral spin liquid with spontaneous breaking of and . However, is still a symmetry. Similarly, a projective symmetry group analysis Wen (2004) shows that broken rotational symmetries can be combined with the broken time reversal to restore an O point group symmetry, ensuring the orbital degeneracy assumed at the outset (see Supplemental Material Sup ()).

Having fixed the ansatz, we can calculate the resulting spectrum of the Majorana fermions. For simplicity, first we focus on the mean-field Hamiltonian for , i.e., the flavors which are decoupled in Eq. (8). In this case the Hamiltonian can be written in the form

 ~HMF=∑k∈12BZγ†kH(k)γk=|t|∑k∈12BZγ†k(hk⋅Σ)γk, (9)

where is the corresponding hopping amplitude in Eq. (8), is a spinor with components labeled by sublattice index, , , and the sum is restricted to half Brillouin zone since Coleman, Miranda, and Tsvelik (1994). As the components of obey , the spectrum is given simply by

 ε±(k)=±|t||hk|. (10)

The dispersion relation is illustrated in Fig. 2. There are two doubly degenerate bands Sup (). Since , the Hamiltonian has a chiral symmetry Wen and Zee (1989); Koshino, Morimoto, and Sato (2014) and the spectrum is symmetric between positive- and negative-energy states. The defining feature of the band structure is the band touching along the edges of the Brillouin zone. These are nodal lines parametrized, e.g., by . Expanding , with , we obtain the effective Hamiltonian on a plane perpendicular to a line node: . The latter is formally equivalent to the Hamiltonian for graphene and yields linear dispersion at low energies with -dependent velocity . These nodal lines can be characterized as topological defects of an SU(2) Berry connection Wilczek and Zee (1984) in reciprocal space (see Supplemental Material Sup ()). The three nodal lines related by C symmetry cross at . Expanding , we find that R is a quadratic band touching point Sun et al. (2009); Herbut and Janssen (2014) with anisotropic dispersion .

Another feature of topologically nontrivial states of matter is the presence of protected surface states. We identify the surface states by calculating the spectrum for in a slab geometry with open boundary conditions in the (111) direction (Fig. 3). There appear two pairs of doubly degenerate bands separated from the continuum, with dispersion terminating at the projections of the nodal lines. Remarkably, the positive-energy surface states are spatially separated from the negative-energy ones as their wave functions are localized at opposite surfaces (which surface depends on the sign of the hopping parameter). This is a direct manifestation of the breaking of inversion symmetry.

We calculate the ground state energy at mean-field level by solving the self-consistent equations that determine the order parameters in Eq. (8). For this purpose we had to diagonalize the Hamiltonian for coupled and fermions in Eq. (8) and found that the spectrum again displays nodal lines Sup (). We obtain . A better estimate of can be obtained by implementing a Gutzwiller projection according to Eq. (6) using VMC Gros (1989); Wang and Vishwanath (2009). Considering a restrictive form of the wave function which neglects variations in the population of the fermionic flavors (see Supplemental Material Sup ()), we obtain . This energy is already comparable to that of the best variational state identified in Ref. Chen, Pereira, and Balents (2010), namely a valence bond solid with . We expect the spin liquid to be stable since small fluctuations of the Z gauge field only induce weak short-range interactions Wen (2004), which are irrelevant in the renormalization group sense for topological semimetals with point or line band touching in three dimensions Herbut and Janssen (2014).

The low-temperature thermodynamic properties are governed by the density of states of the Majorana fermions, which is due to the quadratic band touching point. It follows that the QSL has heat capacity , magnetic susceptibility , and thermal conductivity for . Another important property is the correlation function . We find that vanishes when connects sites on the same sublattice. For vectors connecting different sublattices along (100) directions in the form , where and , the correlation decays at large distances as . This power-law decay coincides with the result for a Dirac point in two dimensions Wang and Vishwanath (2009).

Finally, we address the comparison with available experimental results for the spin liquid candidate BaYMoO. Aharen et al. Aharen et al. (2010) observed that both the heat capacity and the magnetic susceptibility vanishes at low temperatures and have attributed this behavior to a gapped collective spin singlet. de Vries et al. de Vries, Mclaughlin, and Bos (2010) proposed a picture of a valence bond glass, but noted that the muon spin relaxation is comparable to that of QSLs de Vries et al. (2013). Here we propose that an alternative explanation for the vanishing heat capacity and susceptibility at low temperatures is the vanishing density of states of our gapless spin-orbital liquid with nodal lines. A comprehensive study of the properties of this QSL in comparison with experimental results will be presented elsewhere Nat ().

To summarize, we have studied a realistic model for double perovskites in the regime of strong spin-orbit coupling. We proposed a new spin liquid ansatz that gives rise to nodal lines in the spectrum of Majorana fermions. We argued that some experimental results for BaYMoO can be interpreted in terms of the vanishing density of states predicted by our theory. We hope this work will stimulate the search for strongly correlated materials hosting fractional excitations with nontrivial momentum-space topology Hermanns, O’Brien, and Trebst (2015); Maciejko and Fiete (2015).

###### Acknowledgements.
This work was supported by Brazilian agencies FAPESP (W.M.H.N., E.C.A.) and CNPq (E.M., R.G.P.).

## Appendix A Supplemental Material

Supplemental Material for “Chiral spin-orbital liquids with nodal lines”

W. M. H. Natori E. C. Andrade E. Miranda R. G. Pereira

## Appendix B 1. Symmetry of the spin liquid ansatz

Time reversal acts on peudospins and pseudo-orbitals as

 T−1sjT=sj,T−1τx,zjT=τx,z,T−1τyjT=−τyj. (11)

In the Majorana fermion representation for and this transformation can be implemented by . This is equivalent to complex conjugation combined with the operation and . Thus, if we focus on the decoupled flavors , we can take time reversal to be represented simply by complex conjugation.

Let denote a Majorana fermion on site belonging to sublattice . The operators in momentum space are defined by

 γjℓ = √8N∑k∈12BZ[γkℓeik⋅Rj+γ†kℓe−ik⋅Rj], (12) γkℓ = √2N∑j∈ℓγjℓe−ik⋅Rj, (13)

and are normalized such that . For each flavor of Majorana fermion we combine the four sublattice modes into a single “spinor” .

In momentum space, time reversal takes . Up to a hopping amplitude (determined by self-consistent equations, see next section), the mean-field Hamiltonian for a decoupled flavor is of the form with

 H(k)=i⎛⎜ ⎜ ⎜ ⎜⎝0f(kx,ky)f(ky,kz)f(kx,kz)−f(kx,ky)0−f(kx,kz)f(ky,kz)−f(ky,kz)f(kx,kz)0−f(kx,ky)−f(kx,kz)−f(ky,kz)f(kx,ky)0⎞⎟ ⎟ ⎟ ⎟⎠, (14)

where . Notice the factor of . It follows that

 T−1H(k)T=H∗(−k)=−H(k). (15)

We define inversion as the reflection by the mirror plane that exchanges A and C sublattices (plane perpendicular to ). In momentum space, . In addition, we have the action in the internal (sublattice) space given by the matrix (with determinant -1)

 P=⎛⎜ ⎜ ⎜⎝0010010010000001⎞⎟ ⎟ ⎟⎠. (16)

We also define the Z gauge transformation that changes the sign of fermions on the B sublattice:

 GP=⎛⎜ ⎜ ⎜⎝10000−10000100001⎞⎟ ⎟ ⎟⎠. (17)

It is easy to check that inversion anticommutes with the mean-field Hamiltonian:

 (PGP)−1H(Pk)PGP=−H(k). (18)

It follows that the combined transformation is a symmetry of the Hamiltonian:

 (PTGP)−1H(Pk)PTGP=H(k). (19)

The C rotation about a axis that leaves an A site invariant is represented by

 C3=⎛⎜ ⎜ ⎜⎝1000000101000010⎞⎟ ⎟ ⎟⎠, (20)

and the rotation in momentum space takes . In this case a gauge transformation is not required; we obtain immediately that

 C−13H(R3k)C3=H(k). (21)

The C rotation along the axis that exchanges AB, C D is represented by

 C2=⎛⎜ ⎜ ⎜⎝0100100000010010⎞⎟ ⎟ ⎟⎠, (22)

and in momentum space . We need to combine the C rotation with the gauge transformation

 G2=⎛⎜ ⎜ ⎜⎝10000−1000010000−1⎞⎟ ⎟ ⎟⎠. (23)

We then have

 (C2G2)−1H(R2k)(C2G2)=H(k). (24)

Translation by , which we denote by , has the same effect of exchanging sublattices as the above C rotation. Thus, conjugation by , together with in real space, is also a symmetry of the Hamiltonian (and likewise for the equivalent translations in and planes).

Now consider a C rotation along the axis going through an A site, which exchanges C and D sublattices. This can be represented in sublattice space by

 ^C4=⎛⎜ ⎜ ⎜⎝10000100000−10010⎞⎟ ⎟ ⎟⎠. (25)

In momentum space, . Combining with the gauge transformation:

 ^G4=⎛⎜ ⎜ ⎜⎝10000−10000−100001⎞⎟ ⎟ ⎟⎠, (26)

we find

 Extra open brace or missing close brace (27)

Thus, like and , the C rotation inverts the chirality of the ansatz. It is then easy to see that is a symmetry of the Hamiltonian.

The C rotation can be used to construct a symmetry transformation that accounts for the twofold degeneracy of the Majorana fermion bands.

It is also interesting to consider the C rotation that exchanges A and B sublattices, given by

 ^C′4=⎛⎜ ⎜ ⎜⎝0−100100000100001⎞⎟ ⎟ ⎟⎠. (28)

If we define this to be a rotation around axis in the opposite direction than the one in Eq. (25), the transformation in momentum space is , i.e., . It is easy to check that the composition commutes with the Hamiltonian, , and obeys . Thus, we can block diagonalize by sectors labeled by the eigenvalue of the matrix :

 U†MH(k)UM=(~σ⋅hk00−~σ⋅hk), (29)

where and is the unitary matrix that diagonalizes . It is then clear that the spectrum of is twofold degenerate with eigenvalues . Two degenerate states can be distinguished by the eigenvalue of the Hermitean matrix (which is analogous to the chirality of Weyl fermions in the massless Dirac equation).

In summary, the chiral spin-orbital liquid ansatz lowers the symmetry of the Hamiltonian from OZ (where Z is time reversal) to O (where the new group contains combinations of broken point group symmetries with the broken time reversal).

## Appendix C 2. Berry connection

The nodal lines can be characterized as topological defects of a Berry connection in reciprocal space. In our case, the Berry connection has to be non-Abelian due to the double degeneracy of the bands. Away from the nodal lines, we define the SU(2) connection

 Aamn(k)=i⟨ψm(k)|∂kaψn(k)⟩, (30)

where and , , are degenerate eigenstates of (say with energy ) chosen so as to obey and to diagonalize . The generalized Berry phase is the Wilson loop

 U=Pexp[−i∮dkaAa(k)], (31)

where denotes path ordering. The calculation of is simplified if we consider a path around the line node parametrized by , with . For infinitesimal radius , we obtain

 Ax=−12ϵsinασy+O(ϵ0),Ay=12ϵcosασy+O(ϵ0), (32)

which is precisely the singular dependence of a vortex line. We then find , equivalent to a Berry phase.

## Appendix D 3. Solving the mean-field Hamiltonian

In this section we outline the steps required to diagonalize the mean-field Hamiltonian and calculate the ground state energy.

Using the mode expansion Eq. (12), we can rewrite the various hopping terms for Majorana fermions in terms of operators in reciprocal space. For instance,

 i∑⟨j,l⟩∈XYγjAγlB=2i∑k∈12BZh1k(γ†kAγkB−γ†kBγkA), (33)

where is the first component of . The mean-field Hamiltonian becomes

 ~HMF = J18∑k∈12BZ[(2u+¯w)3∑a=1(ηak)†H1(k)ηak−w(θ2k)†H1(k)θ2k+(Θk)†H2(k)Θk] (34) −NJ6+J36∑⟨i,j⟩∈α(3u2ij+3¯wαijuij−wαijvij),

where is the matrix in Eq. (14), is an eight-component spinor that combines and fermions and is an matrix to be specified below.

First consider the fermions , whose spectrum is determined by . Let be the unitary transformation that diagonalizes :

 U†kH1(k)Uk=Λ1(k), (35)

where is a diagonal matrix. The operators that annihilate fermions in eigenstates of are

 ~γkλ = ∑ℓ=A,…,D(U†)λℓγkℓ, (36) γkℓ = 4∑λ=1(U)ℓλ~γkλ, (37)

where is the band index. The mean-field ground state is the state in which all single-fermion states with negative energy are occupied. This leads to the self-consistent equation for expectation values, e.g.

 ⟨γj,Aγj+δxy,B⟩ = 8N∑k∈12BZ[eik⋅δxy⟨γ†kAγkB⟩+e−ik⋅δxy⟨γkAγ†kB⟩] (38) = iIm⎡⎢ ⎢⎣16N∑k∈12% BZ∑λ(U†)λAUBλeik⋅δxy⟨~γ†kλ~γkλ⟩⎤⎥ ⎥⎦ = i2π3Im⎧⎨⎩∑λ(ϵλ<0)∫12BZd3k(U†)λAUBλ⎫⎬⎭,

where in the last line the sum is over bands with negative energy and we took the thermodynamic limit to replace (corresponding to states in the Brillouin zone of the cubic sublattice).

Since determines the spectrum of and fermions, the self-consistent equations for and are the same up to an overall minus sign, depending on the relative sign of the hopping amplitudes and in Eq. (34) We then have the constraint , but must analyze two possibilities, namely and . Without loss of generality (by choosing one of the two degenerate ground states with opposite chiralities), we can set . Numerical evaluation of the integral in Eq. (38) then yields .

The relation between and determines the matrix for . For , we obtain

 H2(k)=uhk⋅Σ′, (39)

where , with

 Kx=12⎛⎜ ⎜ ⎜ ⎜⎝020√320√300√300√3000⎞⎟ ⎟ ⎟ ⎟⎠,Ky=i2⎛⎜ ⎜ ⎜ ⎜⎝0−20√320−√300√300−√3000⎞⎟ ⎟ ⎟ ⎟⎠,Kz=12⎛⎜ ⎜ ⎜⎝cccc1−1−33⎞⎟ ⎟ ⎟⎠. (40)

The components of the matrix vector satisfy the SU(2) algebra. We then obtain the spectrum of and use it to solve the self-consistent equations for and analogous to Eq. (38). In this case of we find a self-consistent solution with and . Having fixed the order parameters, we obtain the mean-field ground states energy .

For we obtain

 H(2)k=uhk⋅Σ′′, (41)

where , with the matrix vector defined in the main text, and the matrices given by , . In this case we find a self-consistent solution with and . The ground state energy is , slightly lower than the result for . This is the value quoted in the main text. We note that the small difference between the two energies may change beyond the mean-field level. However, we have verified that both solutions give rise to a spectrum with nodal lines along MR directions, qualitatively similar to the spectrum for and fermions. Therefore, the properties derived from the low-energy density of states are generic.

## Appendix E 4. Variational Monte Carlo

To check the viability of the proposed chiral spin-orbital liquid beyond the mean-field level, we now enforce the local constraint exactly by considering a Gutzwiller projection of the mean-field wave function (Gros, 1989) by means of a variational Monte Carlo calculation (Ceperley, Chester, and Kalos, 1977).

We begin by rewriting the Majorana fermions in terms of three Dirac fermions, closely following the representation used in Ref. (Wang and Vishwanath, 2009)

 c†i=12(η2i+