Magnetic Order with Fractionalized Excitations in Pyrochlore Magnets with Strong Spin-Orbit Coupling

A recent inelastic neutron scattering experiment on the pyrochlore magnet Yb2Ti2O7 uncovers an unusual scattering continuum in the spin excitation spectrum despite the splayed ferromagnetic order in the ground state. While there exist well defined spin wave excitations at high magnetic fields, the one magnon modes and the two magnon continuum start to strongly overlap upon decreasing the field, and eventually they become the scattering continuum at zero field. Motivated by these observations, we investigate the possible emergence of a magnetically ordered ground state with fractionalized excitations in the spin model with the exchange parameters determined from two previous experiments. Using the fermionic parton mean field theory, we show that the magnetically ordered state with fractionalized excitations can arise as a stable mean field ground state in the presence of sufficiently strong quantum fluctuations. The spin excitation spectrum in such a ground state is computed and shown to have the scattering continuum. Upon increasing the field, the fractionalized magnetically ordered state is suppressed, and is eventually replaced by the conventional magnetically ordered phase at high fields, which is consistent with the experimental data. We discuss further implications of these results to the experiments and possible improvements on the theoretical analysis.

In Fig. S2, we have followed the choice of coordinates as in Ref. 31 , such that the fcc Bravais lattice points are located at the centers of tetrahedra, and the sublattices s = 0, 1, 2, and 3 are displaced by a/8 (1, 1, 1), a/8 (1, −1, −1), a/8 (−1, 1, −1), and a/8 (−1, −1, 1) from the tetrahedral centers respectively, where a is the lattice constant of the conventional cubic cell (which contains four fcc Bravais lattice points). The inversion center is chosen to be the sublattice s = 0 in the unit cell at the origin 0.

S2 Local Coordinates, Spin Hamiltonian, and g Tensor
As discussed in the main text, in the global coordinates, the most general nearest neighbor spin Hamiltonian allowed by the symmetries of pyrochlore lattice contains four independent exchange parameters J 1 , J 2 , J 3 , and J 4 31 . J 1 is the Heisenberg interaction (J), J 2 − J 1 the Kitaev interaction (K), J 3 the symmetric anisotropic exchange interaction (Γ), and J 4 the Dzyaloshinskii Moriya interaction (D). For instance, the interaction between the spins at sublattice 0 and 1 is given by, .
(S3) Figure S1. The sites of pyrochlore lattice form a three dimensional network of corner sharing tetrahedra. The up (down) tetrahedra are colored in red (blue). It is easy to see that each up (down) tetrahedron is surrounded by four down (up) tetrahedra. The underlying Bravais lattice is the face centered cubic (fcc) lattice with four sites (sublattices) per unit cell, which are located at the corners of the tetrahedra. Figure S2. To visualize the full tetrahedral group T d , we embed a tetrahedron in a cube and define a coordinate system with the cubic axes.
It is also a common (arguably much more prevalent) practice to write the spin Hamiltonian (S1) in the local coordinates, where the local z axes are defined along the local [111] directions. The bases in the local coordinates of the four sublattices are defined as 31 In the local coordinates, the spin Hamiltonian takes the form 31 where γ i j and ζ i j are unimodular complex numbers, In the form (S6), the spin Hamiltonian has the advantage that when the spin flip interactions are negligible, i.e. in the limit J ±± −→ 0 and J z± −→ 0, it reduces to a local XXZ model, which is studied in Refs. 6, 10 and shown to support quantum spin liquid states.
To obtain the global exchange parameters in (S1) from the local exchange parameters in (S6), we just have to rotate the local bases (x s ,ŷ s ,ẑ s ) such that they align with the global bases (x,ŷ,ẑ). Call these sublattice dependent rotations R s . We then have, for example, The final result is The g tensor in global coordinates, which is sublattice dependent, can be obtained from that in local coordinates by rotations of the bases similar to the consideration in (S8). That is, We list all the g tensors below for completeness 35 .

S3 Mean Field Hamiltonian
On the one hand, the generic nearest neighbor JKΓ model can be expressed solely in terms of the spinon pairing and hopping channels that signify the spin liquid phase, (S12b) (S12c) The Lagrange multipliers µ 1 , µ 2 , µ 3 ∈ R are introduced in (S12a) to enforce the single occupancy constraint (one spinon per site) on average. Note that we have carefully written the various interactions (S12b)-(S12d) in the form from which a mean field decoupling naturally follows, and the mean field energy is bounded from below (i.e. the stability requirement is satisfied).Ô i j are the bond operators as before, while O i j (without the hat) are variational parameters to minimize the mean field energy.
On the other hand, it is shown in Refs. 35,45 , which provides a group theory analysis of the classical model (i.e. the spins are treated as three component vectors with fixed magnitude), that the nearest neighbor spin interactions on a tetrahedral unit in the pyrochlore lattice can be expressed as a summation of bilinears of the magnetic order parameters m X , multiplied by some energy coefficients a X , Each m X is a linear combination of the components of the spins, while each a X is a linear combination of the exchange couplings. The magnetic order parameters represent different q = 0 classical spin configurations on the pyrochlore lattice.
Since each unit cell contains one up and one down tetrahedra, the Hamiltonian is given by summing 2 times (S16) over the unit cells. As discussed in the main text, we keep only the FM (X = T 1,A ) and AFM (X = E) order parameters (their expressions can be found in the main text) as they are the only relevant classical phases to Yb 2 Ti 2 O 7 . Here we quote their respective energy coefficients, We represent the magnetic order parameters in terms of spinon operators, and carry out a mean field decoupling similar to (S15), where R labels the unit cell (not individual site), and m X (without the hat) are now variational parameters. The stability requirement is satisfied as the coefficients a T 1,A and a E are negative in the J 1 − J 2 phase space (with J 3 = −1 fixed) under study. The above formulation allows us to incorporate both the quantum spin liquid and magnetically ordered states into a single mean field Hamiltonian, with the introduction of a weighting factor r ∈ [0, 1], such that the spin liquid component (involving is multiplied by 1 − r, while the magnetic order component (involving m T 1,A and m E ) by r.

S4 Details of the Spin Liquid Ansatzes
We discuss the Z 2 U and U(1)M spin liquid ansatzes in details, especially the interdependence of the spinon hopping and pairing parameters. The allowed forms of these mean field parameters are dictated by the symmetries of the system. Constraint arises when one symmetry element maps a bond to itself, or two different symmetry elements relates two different bonds. In this section, the term mean field Hamitlonian is referred specifically to as the Hamiltonian (S12a), which has only spin liquid channels, after the mean field decoupling.

S4.1 Z 2 Uniform Ansatz
We first introduce the following 2 × 2 matrix whose components are the spinon creation and annihilation operators 30 , (S19) The spin operator can then be expressed as and the mean field Hamiltonian as where u µ i j are 2 × 2 matrices of the mean field ansatzes. For instance, on the bond 01 , with the exchange couplings J, K, Γ < 0, We also have that enforces the single occupancy constraint (S13). In the form (S20), it is apparent that the spinon representation of spin is invariant under an SU(2) gauge transformation We apply the symmetry operations passively, that is, transform the coordinate axes forward (equivalently transform the vectors backward) 46 , such that where X is an element of the space group and R X is the SU(2) spin rotation associated with X. In the representation (S20), the symmetry transformation (S25) is achieved by 30 wheren is a unit vector along the axis of rotation and φ is the angle of rotation associated with X. Therefore, X acts on the mean field Hamiltonian (S21) by where the SU(2) spin rotation R X has been mapped to the SO(3) rotation O X on the Pauli matrices. The Hamiltonian should be left invariant under X by the definition of symmetry. Taking into account the SU(2) gauge redundancy (S24), this implies that the mean field ansatzes should obey the relations u µ=x,y,z where G X (i) is the SU(2) gauge transformation associated with X at site i. To this end, we list the SO(3) matrices O X associated with some representative elements of the O h point group discussed in Section S1, All other space group elements can be constructed from these, e.g. C y 2 = C . The two fold rotation itself can be obtained by twice the four fold improper rotations, e.g. C x 2 = (S x 4 ) 2 . Note that, since spin is a pseudovector, it is invariant under inversion, hence the SO(3) matrix associated with inversion is the identity. The reflections and improper rotations can be viewed as a combination of rotation and inversion, and their corresponding SO(3) matrices only encode the rotation. For example, the reflection σ d across the plane perpendicular to the [011] direction is a rotation by π about the [011] axis followed by inversion about the intersection of the axis and the plane.
On the other hand, time reversal T acts on the mean field Hamiltonian (S21) by Again, with the SU(2) gauge redundancy, that T being a symmetry requires Recall that the collection of the compound operators G X X (the symmetry group {X} now includes both the space group elements and the time reversal) is known as the projective symmetry group (PSG), and we say that the symmetry X is realized projectively if G X is nontrivial.
In the Z 2 U uniform ansatz, for every symmetry X of the system, we set the corresponding SU(2) gauge transformation G X = 1 to be trivial. We now investigate how the various symmetries limit the form of the spinon hopping and pairing parameters χ i j , ∆ i j , E i j , and D i j . First, time reversal symmetry constrains the singlet parameters χ i j and ∆ i j to be real, and the triplet parameters E µ i j and D µ i j to be imaginary, by (S31a) and (S31b). We also have µ 2 = 0. Next, consider the bond 01 , which is mapped to itself under C x 2 . By (S28a) and (S28b), we have u 0 10 = u 0 01 , u x 10 = u x 01 , u y 10 = −u y 01 , and u z 10 = −u z 01 . As the singlet and triplet parameters obey the relations i j , this implies E x 01 = 0 and D x 01 = 0. The bond 01 is also mapped to itself under the reflection σ [011] d , which, similar to the analysis of the effect of C x 2 above, constrains u z 01 = −u y 01 , or E z 01 = −E y 01 and D z 01 = −D y 01 . Finally, we can use C 3 or other symmetries to relate the mean field parameters on other bonds to those on 01 . For the singlet parameters it is easy, χ i j = χ 01 and ∆ i j = ∆ 01 for all bonds i j by (S28a). For the triplet parameters, we give an example below, or E x 02 = −E y 01 , D x 02 = −D y 01 , E y 02 = 0, D y 02 = 0, E z 02 = E y 01 , and D z 02 = D y 01 , by (S28b). The bond parameters on a down tetrahedron are the same as their counterparts on an up tetrahedron, i.e. u µ i j∈down = u µ i j∈up , by inversion symmetry. There is no further constraint from symmetries, and the number of independent mean field parameters χ 01 , ∆ 01 , E y 01 , and D y 01 in the Z 2 U ansatz is four.

S4.2 U(1) Monopole Flux Ansatz
The analysis of the U(1)M ansatz is in some way easier than that of the Z 2 U ansatz because the pairing terms ∆ i j and D i j are zero. There is no need to introduce the matrix (S19) and write down the mean field Hamiltonian in the form (S21). We have instead  Figure S3. (a) The configuration of link fields a i j , which are the arguments of the singlet hopping parameters χ i j (see (S37)), in the monopole flux ansatz. For each link connecting two sites i and j, a i j is equal to π/2 (−π/2) along (against) the direction of the arrow. This gives a flux of π/2 on each elementary triangle. (b) The ansatz changes under a symmetry transformation X, for example X = C x 2 as shown here. To restore the original configuration of link fields, we apply a sublattice dependent gauge transformation G X = ±1, for example G C x 2 (0) = +1, G C x 2 (1) = −1, G C x 2 (2) = −1, and G C x 2 (3) = +1. The compound operators G X X, which leave the mean field ansatz invariant, form the monopole flux PSG (see Table (S1)).
where u µ i j are now numbers that depend on the hopping terms instead of matrices. For instance, on the bond 01 , with the exchange couplings J, K, Γ < 0, In the spinon representation of spins, for an element X of the space group, the symmetry transformation (S25) is achieved by wheren and φ are as defined previously. However, not all of the 48 elements of the full octahedral group O h are respected in the U(1)M ansatz. The 24 elements that correspond to inversion, reflections (including glide symmetries), and improper rotations are broken, while the 24 elements that correspond to proper rotations (including screw symmetries) are realized within the simple PSG constructed in Ref. 34 , where the site dependent gauge transformations G X = ±1. Translational symmetry is also preserved. The proper rotations are 34 , with the coordinate system defined in Fig. S2, e: the identity; 8 C 3 : rotation by ±2π/3 about one of the local [111] axes (the directions along the center to the corners of the tetrahedron); 3 C 2 : rotation by π about one of the cubic axes (x, y and z directions); 6 C 4 : screw symmetry about one of the axes which are (i) parallel to x axis and going through (0, a/4, 0), (ii) parallel to y axis and going through (0, 0, a/4), and (iii) parallel to z axis and going through (a/4, 0, 0) -rotation by ∓π/2 about one of these axes followed by translation by a/4 along that axis; 6 C 2 : screw symmetry about one of the edges (which connects two sublattices) of a tetrahedron -rotation by π about one of the edges followed by translation along that edge.
It is worth noting that the screw symmetries can be obtained by combining the improper rotations or the reflections with inversion, e.g.  Table S1. The projective symmetry group (PSG) of the monopole flux ansatz. The element G X X is denoted by X for simplicity, where X is one of the 24 proper rotations (including screw symmetries) of the O h point group. The action of G X X is shown in (S39). The subscript s of the fermionic operator f s indexes the sublattice.   Finally, we now extend the monopole flux ansatz to include the triplet hopping parameters, which appears in the mean field Hamiltonian of the nearest neighbor JKΓ model on the pyrochlore lattice (S33), using the relation which can be derived in a way similar to (S27). The SO(3) matrices O X of some representative symmetry elements X can be found in (S29) and (S36). Since inversion symmetry is broken, the bond parameters of the up and down tetrahedra no longer obey u µ i j∈up = u µ i j∈down as in the Z 2 U ansatz. We define v µ i j as u µ i j for the bond i j on a down tetrahedron. The form of the singlet hopping parameter χ i j has already been fixed by (S37). For the triplet hopping parameters on the bond 01 , C x 2 constrains u x 10 = −u x 01 , u y 10 = u y 01 , and u z 10 = u z 01 by (S41), which implies E x 01 is imaginary, while E y 01 and E z 01 are real. u µ i j on other bonds are related to u µ 01 by C 3 , while v µ i j are related to u µ i j by C 4 or C 2 . The number of independent mean field parameters is four. There is no further constraint from symmetries. In the absence of pairing channel, for a free fermion hopping Hamiltonian like (S33) at zero temperature, the single occupancy constraint is satisfied (on average) by half filling of the momentum states, so there is no need to introduce extra Lagrange multipliers (though µ 3 is often identified with the Fermi level in the literature).

S5 Comparisons between the Local and Global Minima from the Mean Field Self Consistent Calculations
We tabulate the energy per unit cell E, the phase, and the reduction of magnetic order parameter in magnitude relative to its classical value S/S 0 (see main text), of the local and global minima, which correspond to a pure spin liquid/spin liquid dominant and pure magnetic order/magnetic order dominant phases respectively, at various magnetic field strength B z , for the Z 2 U ansatz, at Gaulin and Coldea parametrizations (see Tables S2 and S3). A representative value r = 0.23 of the weighting factor is chosen. S/S 0 −→ 0 indicates that the magnetic order is very weak and the system is highly quantum, while S/S 0 −→ 1 indicates that the system approaches the classical limit. In other words, the ratio S/S 0 is a good indicator of the quantumness of the system. As B z increases, the energy difference between the local and global minima grows more significant. Once B z exceeds ∼ 0.01|J 3 |, the spin liquid dominant solution becomes so unfavorable that the self consistent calculations always yield the completely magnetic solution. Similar comparisons are made for the U(1)M ansatz in Tables S4 and S5.

S6 Properties of the Spinon Band Structure in the Pure Magnetically Ordered States
We explain three aspects of the spinon band structure in the pure magnetic phase (where the magnetic order parameters are finite and the spin liquid parameters are zero): (i) flatness (dispersionless), (ii) symmetry about zero energy, and (iii) four fold with the coefficients c µ s ∈ R. Since neither spinon hopping nor pairing at two different sites is present, Fourier transform does not introduce any nontrivial phase factor e ik·(R i −R j ) in the second equality of (S42). This explains the flatness of the spinon bands as the energy eigenvalues are independent of the momentum k. Furthermore, written in the basis ( f k,0,↑ , f k,0,↓ , . . . , f k,3,↑ , f k,3,↓ ), the Hamiltonian matrix is an 8 × 8 block matrix whose nonzero blocks are the four 2 × 2 matrices along the diagonal. Diagonalization yields the energy eigenvalues ω ks = ± (c x s ) 2 + (c y s ) 2 + (c z s ) 2 .
The ± sign means that the spinon bands are symmetric about the zero level. Finally, from (S43) we see that the energy eigenvalues depend on the coefficients c µ s only through the second power. We examine the FM and AFM order parameters and find that their respective set of coefficients c µ s satisfies c µ s = ±c µ s for different sublattices s and s . This implies the four fold degeneracy.