Designing Quantum Spin-Orbital Liquids in Artificial Mott Insulators

Quantum spin-orbital liquids are elusive strongly correlated states of matter that emerge from quantum frustration between spin and orbital degrees of freedom. A promising route towards the observation of those states is the creation of artificial Mott insulators where antiferromagnetic correlations between spins and orbitals can be designed. We show that Coulomb impurity lattices on the surface of gapped honeycomb substrates, such as graphene on SiC, can be used to simulate SU(4) symmetric spin-orbital lattice models. We exploit the property that massive Dirac fermions form mid-gap bound states with spin and valley degeneracies in the vicinity of a Coulomb impurity. Due to electronic repulsion, the antiferromagnetic correlations of the impurity lattice are driven by a super-exchange interaction with SU(4) symmetry, which emerges from the bound states degeneracy at quarter filling. We propose that quantum spin-orbital liquids can be engineered in artificially designed solid-state systems at vastly higher temperatures than achievable in optical lattices with cold atoms. We discuss the experimental setup and possible scenarios for candidate quantum spin-liquids in Coulomb impurity lattices of various geometries.


Wavefunction of the strong coupling subcritical regime
The wave function Ψ(r) of a two dimensional massive Dirac fermion moving around a Coulomb impurity satisfies, where V (r) = −g/r, r > a −g/a, r ≤ a is the regularized Coulomb potential. Here we set = v = 1. The two-component wave function is he Dirac equation with the presence of an external potential becomes where κ = j − 1 2 .

Solution for r > a
The wave functions assumes the form (j indexes will be dropped) where ρ = 2λr, β = √ m 2 − 2 , andγ = j 2 − g 2 . After some algebra, we get These equations can be decoupled both of which are Kummer's differential equation with the solution y = cF(a; c; x). Here we call the confluent Hypergeometric function of the first kind 1 F 1 (a; c; x) as F(a; c; x), and notice that 1 F 1 (a; c; x = 0) = 1. The solutions are Weak coupling regime When g < j, γ is real. Integrability of the wavefunction at ρ → ∞ requires that d 1 = d 2 = 0. From Eq.(5), one can determine the ratio To simplify the notation we callm = mg/β,˜ = g/β, and That leads to the solution This solution is regular at ρ → 0, and the short distance cut-off can be set to zero.

Strong coupling regime
When the coupling g > 1 2 , γ becomes imaginary. With the requirement of imposing a small distance cut-off, the condition that the wave function behaves well at ρ = 0 is not necessary, so we should include both ±γ branches into the solutions. The ratio between the two And in this case we hope the wave functions die off at ρ → +∞, which can also serve to settle down the ratio between γ−branch and (−γ)-branch. The formula can be used here is the asymptotic form of the hypergeometric function, The second part is required if the gamma function Γ(a) is infinite (when a is a negative integer) or Re(z) is nonpositive. In our case, we could exclude these two conditions, and only keep the second term. Therefore for large |z|, the dominating part (which is growing) of F(z) is and we ask for some condition to cancel this term. For aF (ρ; γ) + bF (ρ; −γ) we need the following two terms to be finite at ρ → +∞ From 9, We can assign The solution for F (±) is

Solution for r ≤ a
In the r < a region, we define the Dirac equation becomes where E ± = + g a ± m. These equations can be decoupled into The solutions are and

Energy
The energy can be determined by matching the inside solution and the outside one, formally through Out(j, , r, g) r=a = Ins(j, , r, g) r=a For given j, g, and at r = a, we can determine the energy [2] V.M. Pereira, V.N. Kotov, and A.H. Castro Neto, Phys. Rev. B 78, 085101 (2008)