A Z spin-orbital liquid state in the square lattice Kugel-Khomskii model
We argue for the existence of a liquid ground state in a class of square lattice models of orbitally degenerate insulators. Starting with the SU(4) symmetric Kugel-Khomskii model, we utilize a Majorana Fermion representation of spin-orbital operators to access novel phases. Variational wave functions of candidate liquid phases are thus obtained, whose properties are evaluated using variational Monte Carlo. These states are disordered, and are found to have excellent energetics and ground state overlap () when compared with exact diagonalization on 16 site clusters. We conclude that these are spin-orbital liquid ground states with emergent nodal fermions and Z gauge fields. Connections to spin 3/2 cold atom systems and properties in the absence of SU(4) symmetry are briefly discussed.
In correlated insulators, the degrees of freedom that remain at low energies are spin and orbital degeneracy. At low temperatures one usually obtains an ordered state described essentially by a classical variable, the Landau order parameter. Ground states that are not described by the Landau framework are expected to possess strikingly new properties. While they are known to occur in one dimensional systems, an important question is whether they arise in bulk two and three dimensional materials. Theoretical studies have largely focused on quantum spin systems. While model Hamiltonians for spin liquids exist, one needs special conditions like strong frustration to ensure that the spins do not order. Otherwise, even for spin 1/2 quantum fluctuations are not typically strong enough to destroy order. On the other hand, orbital degeneracy in insulators can enhance quantum fluctuationsZaanen (); Khallulin (); BosscheZhangMila (), destroying order and possibly lead to a spin-orbital liquid state. An experimental illustration is provided by the insulating spinels MnScS and FeScS. The former, a pure S=5/2 system, magnetically orders below Kelvin. The latter, a spin S=2 system, is identical in most respects except that it involves a twofold orbital degeneracy. In contrast to the spin system, it is found to remain spin disordered down to the lowest temperatures of mKelvin Loidl1 ().
Here, we theoretically study a simple spin 1/2 square lattice model with the minimal twofold orbital degeneracy. Such a spin-orbital Hilbert space is realized for e.g. with the orbitals of Ag, Cu or (low spin) Ni in an octahedral environment. We focus on a model that captures the effect of enhanced quantum fluctuations from orbital degeneracy, and argue for a liquid ground state in this case. We give theoretical arguments for the stability of this phase as well as well as a numerical Monte Carlo study of a variational wave function. The latter is found to have extremely good energetics for our Hamiltonian, and allows us to characterize this phase beyond the simple fact that it is disordered. The low energy collective excitations of the liquid state are captured by an emergent Z gauge field, coupled to Dirac like fermionic excitations with fractional spin and orbital quantum numbers. The ground state wavefunction is strongly entangled, and can be thought of as a product of three Slater determinant wave-functions.
Realistic spin-orbital Hamiltonians tend to be rather complicated with several exchange couplings that are strongly direction dependent. Moreover, a linear coupling between the orbital degrees of freedom and the Jahn-Teller phonons can quench coherent orbital dynamics. However, for sufficiently strong exchange interactions, the coupling to phonons can be ignored, and the orbitals can be taken to be quantum degrees of freedom. Since we are interested here in the general effects of enhanced quantum fluctuations from orbital degeneracy, we follow BosscheZhangMila () and others in considering a model that treats all four states on a site symmetrically, i.e. the SU(4) symmetric spin orbital model. This will allow for comparisons and is a useful starting point. We show later that our essential conclusions are unchanged on perturbing away from this high symmetry point.
The high symmetry SU(4) model may, in fact, have a direct physical realization. We point out that a model of e orbitals on certain high symmetry lattices, like the diamond lattice, can have spin orbital Hamiltonian that are nearly SU(4) symmetric. A different setting for this physics has been opened up by the recent experimental developments on the trapping and cooling of alkaline earth atoms. Confining these to the sites of an optical lattice leads to SU(N) symmetric magnetic models. The nuclear spin here provides the flavors, and, given the weak dependence of scattering lengths on nuclear spin, leads to SU(N) symmetric exchange interactions, which, for fermionic atoms, will be antiferromagneticHermele (). Fermionic alkaline-earth like atoms of Yb have been cooled to quantum degeneracy Fukuhara1 (), while the Mott state of the bosonic Yb has been recently realized Fukuhara2 (). A different realization may be provided by spin 3/2 cold atoms, such as Cs confined to the sites of an optical lattice. At unit filling, one has four states per site, and although the physical symmetry is only that of spin rotations, the small difference in scattering lengths imply only a weak breaking of SU(4) symmetry. It was pointed out in Wu () that even including these differences only breaks the symmetry down to SO(5)Z in the low energy limit.
Ii The Model:
where are the spin-1/2 Pauli matrices and are the Pauli matrices acting on the two degenerate orbital states. We consider the antiferromagnetic(AFM) case () and set hereafter. The high symmetry of this model implies that the three spin operators , three orbital operators and nine spin orbital operators all appear with equal weight. These are fifteen generators of SU(4). It should be noted that symmetry does not uniquely define a model, one also needs to specify the representation of the symmetry group appearing at each site. Here the fundamental representation appears and an SU(4) singlet can only be formed between four sites.
The model (1) was numerically studied in BosscheZhangMila () using exact diagonalization (the model suffers from a sign problem in the spin-orbital basis) on system sizes of up to 44. The ground state () is an SU(4) singlets at zero wavevector. Simple trial states, such as with spin-orbital order , or a box singlet state, have much higher average energy ( and respectively) pointing to the importance of quantum fluctuations. This motivates the study of spin orbital liquids as candidate ground states.
The idea of resonating valence bonds PWA (); RKS (), provides an intuitive picture of the quantum spin liquid. A more formal approach that is easier to generalize is the slave particle formalism, where the spin is decomposed into ‘partons’, e.g. , where the are the Pauli matrices, and creates a boson (Schwinger boson) or fermion with spin at site , and the constraint is imposed at every site. One then makes a mean field decomposition to obtain a quadratic Hamiltonian, and the constraint is then imposed by projecting the wave function to obtain a variational ground state. The state so obtained is a candidate spin-liquid wave function. This procedure can be generalized in a straightforward way to the spin-orbital Hilbert space at hand, by introducing a four component and the constraint above at every site ShenZhang (); Mishra (). The physical spin and orbital operators are again bilinears of . However, there are some drawbacks to this straightforward generalization. The bosonic parton representation cannot treat the SU(4) symmetric point, while the fermionic parton theories necessarily lead to Fermi surface states which can be hard to stabilize as ground states.
Below, we will sidestep these difficulties by showing that this problem admits a third, physically distinct, parton representation in terms of Majorana Fermions. This representation, which has not previously been applied to two dimensional systems, offers us many advantages. Besides being more economical (in terms of expanding the Hilbert space in the minimal fashion), it leads to liquid states with a gauge group, whose low energy physics is well understood and know to exist as stable phases. We emphasize here that the projected wavefunction obtained from this Majorana parton representation of the SU(4) model is distinct from the ’Schwinger’ fermion representation, involving four fermionic operators.
Iii Majorana Parton Formulation:
We first point out the group isomorphism SU(4)SO(6). The 15 generators of the latter are dimension six antisymmetric real matrices , where and . An operator representation of this algebra is obtained by introducing six Majorana fermions which satisfy the anticommutation relations . The operators , where summation over repeated indices is assumed, reproduce the commutation relations for SO(6) generators. The Majorana Fermions transform as SO(6) vectors.
We now use the group isomorphism to obtain a representation of the spin-orbital operators, in terms of Majorana fermions. It is helpful to write the set of Majorana fermions as a pair of three component vectors where, e.g. and we have introduced site indices . The spin and orbital operators can then be written in the compact form:
and , which automatically obey the expected algebra. Note, the sign of the Majorana fermion operators can be changed without affecting the physical operators. This Z redundancy is connected to the fact that the Hilbert space is now enlarged - since each Majorana fermion corresponds to degrees of freedom, we have states whereas there are only physical states per site. The excess states can be removed by implementing a Z constraint at each site: first define the operator commutes with the physical operators (which are fermion bilinears) and is idempotent .
Hence implementing the constraint in (3), restricts us to the physical Hilbert space. This operator generates the gauge transformation on the Majorana fermions .
The model in (1) can be written in these variables as:
When supplemented by the constraint (3), this is an exact rewriting of the model. Here, the symmetry of the model is explicit. The quartic nature of the Hamiltonian requires an approximation. We begin with a mean field treatment and use it to generate variational states in which the constraint is treated exactly.
In the context of models with only spin 1/2 (and no orbital degrees of freedom) we note that a representation utilizing three Majorana fermions per site, where the spin operator is given by an expression identical to that in Eqn. 2, has been studied Tsvelik (). However, we point out that this is distinct from our current formalism since a single site constraint cannot be applied to generate the physical Hilbert space. Nevertheless, this provides an alternate parton approach - for example, with an even number of sites, one can view half of the spins as ‘orbital pseudo-spins’, then the spin problem will be artificially converted to a spin-orbital problem, (although without SU(4) symmetry). Then the gauge fixing, and construction of complex fermions and variational wave functions can be proceeded in the same fashion as in this paper. But this artificial discrimination of spin and ‘orbital pseudo-spin’ usually will superficially break lattice symmetry. Another way to construct a Hilbert space for the spin 1/2 only model, is to introduce a fourth Majorana fermion on every site, and the set of four fermions satisfies a product constraint as in Eqn. 3 Kitaev (). This has the benefit of being formulated with a unique single site constraint. However, unfortunately it turns out that these four fermions are just the real and imaginary parts of the two Schwinger fermion operators, and do not lead to a new representation. The constraint is the familiar one of requiring single occupancy of the Schwinger fermions.
Mean Field Theory and Gutzwiller Projection: With real mean field parameters (), we have:
In self consistent mean field theory . For convenience we combine the 6 Majorana fermions into 3 complex fermions: which are more intuitive although the SO(6) symmetry is no longer explicit. The constraint then is: (while the odd values of the site occupation are forbidden). Writing , we have i.e. the mean field theory simply involves fermions hopping with pure imaginary amplitudes. Such a band structure is automatically particle-hole symmetric, which leads to half-filled bands for each of the fermions. Note, despite the imaginary hoppings the mean field ansatz is time reversal symmetric if the hopping is bipartite. The mean field wave function is simply a product of three identical Slater determinants. While the specific Slater determinant depends on the mean field ansatz, we make a few general observations below. If we consider a system with sites, required to obtain an SU(4) singlet state, each Slater determinant is a function of particle coordinates, corresponding to half filling. Gutzwiller projecting the mean field state into the constrained Hilbert space, yields a physical spin-orbital wave-function. In the fermion representation, a site can either have no fermions (denoted by ), or two fermions, in which case there are three states, . These are related to the spin orbital basis states via:
Given a configuration specified by the locations of the , and states (at sites }, and respectively, where are distinct positions), the spin-orbital wave function assigns an amplitude to it. Note, the locations of the states are automatically specified. For an SU(4) singlet we need equal numbers, , of the four types of sites, so etc. After the Gutzwiller projection we obtain:
Thus the projected spin-orbital wave function is a product of three Slater determinants with a lot of entanglement. We now apply this formalism to specific models.
Iv One Dimensional Chain
This SU(4) symmetric nearest neighbor model in 1D is very well understood and serves as a benchmark for our technique. The only symmetric mean field ansatz is a uniform , leading to a dispersion . We construct the resulting projected wave function for site chain with and antiperiodic boundary conditions, and evaluate its properties using variational Monte Carlo. The energy per site from the projected wave functions extrapolated to the thermodynamic limit is , not far from the exact resultSutherland (), . The leading term in the asymptotic spin correlation function is , consistent with theoretical and numerical predictions Azaria (); Affleck (); Pati (). Note, these desireable properties of the wavefunction only arise after projection. (see TABLE 1).
Interestingly, the wavevector of the dominant correlations do not correspond to a natural wavevector of the mean field dispersion, and arises entirely from projection. In contrast, this wavevector is easier to understand on projecting a quarter filled band, which arises in the standard fermionic representation of spin-orbital operators . Remarkably, one can show that the projected wave functions arising from this representation and the Majorana fermion representation discussed above, are identical in one dimension. We stress that this is a special feature of one dimension, and in higher dimensions, the two will lead to physically distinct states. Details of the proof can be found in Appendix A.
V Square Lattice:
The mean field states on the square lattice can be distinguished by the gauge invariant flux through the elementary plaquettes e.g. for the plaquette . Translation and Time reversal symmetry dictates that this flux must be uniform and can be either or . This leads to two distinct mean field states the uniform and flux state ansatz. The uniform states ansatz is , where is the position of lattice sites. The mean field dispersion is for all the three flavors, and has a square Fermi surface. However, the uniform ansatz state has higher energy than the -flux ansatz both in mean field theory and after projection, so we focus on the -flux ansatz.
The -flux ansatz is , as shown in FIG. 1. The unit cell in mean field theory is doubled, with and sublattices as shown.
After a Fourier transform, the mean field Hamiltonian (5) is:
where the sum over is over the (for lattice) -points in the reduced() Brillouin zone(BZ), and . The above result can be further diagonalized by a Bogoliubov transformation and produce the two branches of the mean field dispersion . This dispersion has two Dirac nodes at and , with isotropic dispersion in their vicinity. Including flavor indices, we thus have 6 two-component Dirac fermions.
We use anti-periodic boundary conditions in both directions for lattices( even) lattices. The -points are then , which avoids zero energy modes. Filling th negative energy states gives us a Slater determinant mean field wave function for each of the three fermion species. The Gutzwiller projected wave function is then easily written down as (6), in terms of this Slater determinant. Evaluating its properties however requires a numerical variational (determinantal) Monte Carlo approach Gros (); Ceperley (). We generate a random initial basis state having significant overlap with the mean field wave function. Random pairs of sites are selected and updated with the Metropolis rejection rule. 10000 ‘thermalization’ sweeps( pairwise updates) are performed before measurements of physical quantities. Measurements are done in 100000 sweeps. The entire process is repeated 10 times to ensure stability of results.
Energetics: The energy from the projected wave function is listed in TABLE 2 for up to 20. We notice that for the lattice our total energy is close to the ground state energy obtained in the exact diagonalization studyBosscheZhangMila (). More importantly, our energy lies below the first excited state energy obtained in that study, which already implies a significant overlap () between our wave function and the exact ground state wave function. Note, our ‘variational’ wave function has no variational parameter - which makes this agreement more remarkable, especially given that there are 24024 SU(4) singlet states already at this system sizeBosscheZhangMila ().
|exact ground stateBosscheZhangMila ()||-17.35|
|exact first excited stateBosscheZhangMila ()||-16|
|projected SO(6) Majorana fermion mean field||-16.57|
|projected SU(4) Schwinger fermion mean field||-6.38|
|orbital ferromagnetic, spin AFM stateLi (); Gross ()||-14.46|
|box(plaquette) ordered stateLi ()||-12|
|four-sublattice SU(4) Neel stateLi ()||0|
We have checked several other simple states on this lattice, they all show much higher energy, compared to the first excited state in the exact spectrum. The comparison is shown in TABLE 3.
Wavefunction properties: We now provide evidence that the resulting wave function is a spin-orbital liquid. First, we would like to establish that it has no conventional order, to clearly show it is not a conventional state. Next, the specific type of liquid state being proposed - with emergent Dirac fermions and Z gauge fluxes - needs to be established.
We first check for spin-orbital order. Given the SU(4) symmetry, it is sufficient to compute the correlations which are found to be rapidly decaying in space. The structure factor for various lattice size was computed, e.g. FIG. 2 shows the result for the lattice. A broad maximum at is seen, but no Bragg-peak develops. We thus exclude magnetic-orbital-ordering.
A more likely order is an SU(4) singlet state that breaks lattice symmetry. This is the analog of the Valence Bond Solid order for SU(2) magnets. However, SU(4) singlets require at least four sites so box crystalline orders may arise. Two such natural orders were proposed by Li, et al.Li (). In both, the SU(4) singlets are formed on 1/4 of the elementary plaquettes, but these are either arranged in a square lattice, or in a body centered rectangular lattice with aspect ratio 2. Both box orders have bond energies modulated at the wavevectors or . We first check if our wave function has these correlations by defining . Then for a lattice should scale as , if long range order is present. For example, for the perfect square box state which is a product state of SU(4) singlets, . In the absence of order however, this quantity will scale as . For we did not observe scaling but rather good scaling, as shown in Column 3 of TABLE 2 indicating no sign or box order upto site systems. An independent check is provided by modulating the mean field parameters to realize the box orders. The average energy of the projected state is then compared against the unmodulated wave function. We find that the energy always increases, for both kinds of orders, pointing to the stability of the unmodulated state. Within mean field theory alone, however, the state is locally unstable to box modulation AffleckMarston (). The more reliable projected energy study however point to the opposite conclusion.
Motivated by a recent proposal of chiral SU(N) states in large-N limitHermele (), we also consider a chiral state on the square lattice. We add to the -flux state pure imaginary hopping of fermions on the diagonal bonds in such a way that each triangle has flux. This is a particular mass term of the Dirac fermions and it opens a gap in the mean field dispersion, similar to the box order mentioned above. Indeed the mean field energy decreases with the diagonal hoppings. However after projection the energy always increase after adding this term, indicating stability against this chiral order. This distinction between mean field and projected mean field energetics has been observed in other projected wave function studies as wellYingRan ().
We expect the spin-orbital liquid to be a nodal Z state, i.e. it contains emergent Z gauge fields and nodal Dirac fermions that behave like free particles at low energies. Establishing this directly is more challenging - it is well known that observing the Z topological order of projected wave functions in the presence of gapless gauge charged fermions is tricky IvanovSenthil () and left to future work. Free nodal fermions would lead to spin and orbital correlations that decay as , which we check for by computing in a size system. The fast decay limits us to . The results are shown in column 4 of TABLE 2. There are strong commensuration effects which reduce the correlation when is an even number. However, for the other three values of the correlation seems to show the required scaling. Another indirect evidence for such fermionization is the nature of these spin correlations in the presence of a Zeeman field , which leads to a shifted chemical potential for one fermion specie. We find that the projected wave function now has a ring of incommensurate correlations around (FIG. 3 and FIG. 4).
Breaking SU(4) symmetry It is natural to ask if the physical conclusions derived above are stable when enlarged SU(4) symmetry of our model is lost. Since the gauge fluxes are gapped, a weak perturbation cannot lead to confinement. Also, the gapless nodal fermions are actually protected by discrete symmetries - one needs to break lattice symmetry to gap the nodes, as can be seen from an analysis of the fermion bilinear terms(see Appendix B for details). The only physical difference that arises in the lower symmetry case is that the chemical potential of the fermions may not be at the nodal points. Hence the SU(4) symmetry is not essential to our conclusions. In futureFWAV (), we will apply this analysis to realistic Hamiltonians with reduced symmetry and search for liquid phases in those regimes.
Vi Physical Realizations
Since most natural spin orbit Hamiltonians are not near the SU(4) symmetric point, it is natural to ask where we might expect to find a model where the SU(4) symmetry is even approximately realized. As discussed in the introduction, cold atom systems in optical lattices provide some promising direction for realization, if sufficient cooling of those magnetic Hamiltonians can be achieved. In this section, we point out that even in solid state systems, approximate SU(4) symmetry may be achieved, on certain high symmetry lattices. In particular, we point out that on the diamond lattice, if only nearest neighbor exchange is considered, the interactions are close to the SU(4) point. Exchange interaction arises from a combination of hopping and on site interaction. The main observation is that due the high symmetry of the diamond lattice, hopping matrix elements must be SU(4) symmetric. The onsite interactions deviate from SU(4) symmetry due to, for e.g., the Hunds interaction. However, these are typically a fraction of the overall repulsion leading to nearly SU(4) symmetric exchange.
For a system of orbitals on the diamond lattice with full lattice symmetry and without spin-orbital coupling, we first prove that the electron hopping matrix elements on nearest neighbor bonds have SU(4) symmetry. Denote the creation operators of the two orbitals as and respectively, where is spin index, is site index. The hopping amplitude on bond is generically a matrix , and the process is described by the term . Consider a bond from origin along the direction, the reflection does not change this bond, however the orbitals transform non-trivially . Since this reflection is a physical symmetry, the electronic Hamiltonian should be invariant under its action, thus we get . Similarly consider a threefold rotation , it does not change the bond as well, but the orbitals transform as . Then we get . These two conditions on ensure that is proportional to identity matrix. Thus we have proved that the hopping on direction preserves both orbital and spin, namely is given by the term . By lattice symmetry we conclude that all other nearest neighbor bonds have this property. Therefore the nearest neighbor hoppings have SU(4) symmetry. One should note that this proof cannot be extended to next nearest neighbor and other generic hoppings.
However, the Coulomb interaction of these orbitals generically does not have SU(4) symmetry. The onsite Coulomb interaction is given by the Kanamori parametersKanamoriParameters (), , and approximately . The SU(4) symmetry is present only if and . We usually expect that , then SU(4) is an approximate symmetry of the Hubbard model and thus an approximate symmetry of the derived spin-orbital exchange model.
This type of systems may be realized in certain A-site spinels, where the magnetic ions form a diamond lattice, and when only the e orbitals are active. One caveat is that the spinel structure allows for the next nearest neighbor exchange strength to be fairly large and even comparable to the nearest neighbor one, which may significantly break the SU(4) symmetry. Interestingly, the experimentally discussed ‘spin-orbital’ liquid candidate, FeScS, is an e system on the diamond lattice Loidl1 (). However, it differs in two important respects from the ideal model considered here. First, there is a magnetic moment on each site that is Hunds coupled to the e fermion, and second, the further neighbor exchange interactions are believed to be substantial in this material.
We acknowledge support from NSF-DMR0645691, and discussions with M. Hermele.
Appendix A Proof of the equivalence between projected SO(6) Majorana mean field state and projected SU(4) Schwinger fermion mean field state for 1D chain.
Consider a -site chain with periodic boundary condition. The mean field wave function for the Majorana fermion representation is
where is fermion vacuum, is the Fourier transform of real space fermion operator . For a physically allowed real space configuration mentioned above , the overlap with the mean field wave function is
where is the Slater determinant for one fermion specie with the following matrix elements ()
Therefore the determinant is
where . And the overlap is
The mean field Hamiltonian for the standard ‘Schwinger’ fermion representation is
where is the particle-hole conjugate of the Schwinger fermion . The particle-hole transformation on the fourth specie is required for the following projected wave function to represent a bosonic spin-orbital wave function. The quarter-filling mean field wave function for the standard fermionic representation is (we assume is even for simplicity)
where is the Fourier transform of the real space SU(4) ‘Schwinger fermions’ , and is the fermion ‘vacuum’ that can be annihilated by and the particle-hole conjugate of the fourth specie .
A physically allowed real space configuration in this representation is still labeled by three sets of distinct numbers ,
The overlap of this configuration with the mean field wave function is the product of four Slater determinants,
The Slater determinants has the following matrix elements(), . The matrix elements of the Slater determinant is ()
Therefore the determinant of is
And the determinant of is
Finally the overlap between this basis state and the mean field wave function is (note that )
Therefore these two projected wave functions for 1D chain are identical. The crucial property utilized was that the Slater determinants appearing are Vandermonde determinants in one dimension. This property does not hold in higher dimensions.
Appendix B Projective symmetry group analysis of fermion bilinears in the -flux state on square lattice
For the -flux ansatz in FIG. 1 on an infinite lattice the lattice group symmetries are realized as follows(flavor index omitted):
Define a 4-compoent field
linearize the dispersion around the Dirac point, the low-energy Hamiltonian becomes
where are Pauli matrices acting on the 2D space. The original fermion in real space can be represented as where . Then we have the transformation property of the field
where are Pauli matrices acting on the 2D space of the two Dirac nodes and . The low energy Hamiltonian is invariant under these transformations. In the following we will prove that these symmetries prohibit mass term and velocity anisotropy in the low energy theory.
Consider a general mass term where is a constant non-trivial(not proportional to identity) Hermitian matrix. and translation symmetries require that with be a constant Hermitian matrix. reflection symmetry requires , but this violate the rotation symmetry and is forbidden.
Consider a general velocity anisotropy term with constant matrices and . Again and translation symmetries require that . reflection symmetry requires and , therefore . Also consider 180 rotation symmetry , it requires , thus . Now we have with constants . Use the 90 rotation symmetry we get