Holon Wigner Crystal in a Lightly Doped Kagome Quantum Spin Liquid

TL;DR

This study uses large-scale DMRG to show that light doping in a Kagome quantum spin liquid leads to a holon Wigner crystal, challenging expectations of superconductivity or exotic metals.

Contribution

It provides the first large-scale numerical evidence that light doping of a Kagome spin liquid results in a holon Wigner crystal, not a superconductor or topological metal.

Findings

Doped holes form an insulating charge-density wave with one hole per unit cell.

Spin correlations remain short-ranged, indicating a holon crystal rather than a hole crystal.

Results suggest relevance to doped Herbertsmithite materials.

Abstract

We address the problem of a lightly doped spin-liquid through a large-scale density-matrix renormalization group (DMRG) study of the $t$-$J$ model on a Kagome lattice with a small but non-zero concentration, $\delta$, of doped holes. It is now widely accepted that the undoped ($\delta=0$) spin 1/2 Heisenberg antiferromagnet has a spin-liquid groundstate. Theoretical arguments have been presented that light doping of such a spin-liquid could give rise to a high temperature superconductor or an exotic topological Fermi liquid metal (FL$^\ast$). Instead, we infer that the doped holes form an insulating charge-density wave state with one doped-hole per unit cell - i.e. a Wigner crystal (WC). Spin correlations remain short-ranged, as in the spin-liquid parent state, from which we infer that the state is a crystal of spinless holons (WC$^\ast$), rather than of holes. Our results may be relevant to Kagome lattice Herbertsmithite $\rm ZnCu_3(OH)_6Cl_2$ upon doping.

Figures

Equations

$$H=-t\sum _{⟨ij⟩\sigma }(c_{i\sigma }^{†}{c}_{j\sigma }+h.c.)+J\sum _{⟨ij⟩}\left({\mathbf{S}}_i\cdot {\mathbf{S}}_j-\frac{1}{4}{n}_i{n}_j\right)$$

$${\Phi }_{\alpha \beta }(r)=\frac{1}{L_y}\sum _{y=1}^{L_y}|⟨{\Delta }_{\alpha }^{†}(i_0){\Delta }_{\beta }(i_0+r)⟩|.$$

Full text

Introduction: Broad interest in quantum spin liquid phases (QSLs) was triggered by the notion that they can be viewed as insulating phases with preexisting electron-pairs, such that upon light doping they might automatically yield high temperature superconductivity.Anderson (1987); Kivelson et al. (1987); Rokhsar and Kivelson (1988); Laughlin (1988); Wen and Lee (1996); Balents (2010); Senthil and Lee (2005) It has also been proposed that a doped QSL might form an exotic topologically ordered Fermi liquid state (known as an $\rm FL^{\ast}$ state)Senthil et al. (2003); Punk et al. (2015), or various other topologically ordered versions of other familiar phases.Patel et al. (2016) More broadly, it has been suggested that a host of behaviors of highly correlated electronic systems can be best understood from the perspective of doped spin-liquids.Lee et al. (2006); Fradkin et al. (2015); Tocchio et al. (2013); Gull et al. (2013) However, a “microscopic” theory of QSLs is exceedingly difficult, as they seem to arise only in narrow portions of the generalized phase diagram where more typical broken symmetry states are suppressed by frustration, and in an “intermediate coupling” regime where neither the effective kinetic nor the interaction energy is dominant.

The spin 1/2 antiferromagnet on the Kagome lattice (depicted in Fig.1) with nearest-neighbor (NN) Heisenberg interactions, i.e., $H_{J}$ in Eq.(1), is geometrically frustrated. A number of numerical simulations Jiang et al. (2008); Yan et al. (2011); Jiang et al. (2012); Depenbrock et al. (2012); Gong et al. (2015); Mei et al. (2016) have provided strong evidence that its ground state is a “Z2-QSL” state with exponentially falling spin-spin correlations and a non-zero spin-gap, although some recent studiesRan et al. (2007); Iqbal et al. (2013, 2014); Liao et al. (2016); He et al. (2016) have suggested that the true ground-state in the 2D limit may be a gapless (nodal) QSL. The fact that the observed spin correlation lengths in the earlier references are short compared to the width of the ladders studied leaves little room to doubt that they reflect the properties of a spin-gapped state. None-the-less, it is plausible that there are at least two distinct QSL phases - one gapped and another ungapped - that are very close in energy such that the balance between them can shift as a function of ladder width, geometry, or slight changes in parameters; if this is the case, it could reconcile the two sets of findings while leaving open the issue of which QSL is the ground-state in the 2D limit.

Independent of which QSL has the lowest energy in 2D, in the present study, the fact that the spin correlation lengths we observe are several times shorter than the width of our cyllinders leaves little doubt that we are studying the properties of a doped, fully gapped $Z_{2}$ spin liquid. Experimentally, the celebrated material Herbertsmithite [$\rm ZnCu_{3}(OH)_{6}Cl_{2}$] is a realization of the two-dimensional Kagome antiferromagnetHering and Reuther (2016); Jiang et al. (2017) where the copper ions carry spin S =1/ 2 magnetic moments which condense to form a QSL groundstate. Specifically, experimental evidence of fractional spin excitations has been found in neutron scattering and strong indications of a spin-gap are seen in NMR studies of single crystals. Han et al. (2012); Fu et al. (2015)

The elementary excitations of a Z2-QSL can be constructedRokhsar and Kivelson (1988); Senthil and Fisher (2001) as combinations of a fermionic charge 0 spin 1/2 “spinon,” a charge e spin 0 bosonic holon, and a neutral topological “vison.” The statistics of these particles is a matter of convenience - for instance a fermionic holon can be constructedKivelson (1989) as a boundstate of a bosonic holon and a vison, while a normal spin 1/2, charge e hole can be constructed as a boundstate of a spinon and a holon. (We consider only hole-doping of the QSL, so we can safely ignore negatively charged excitations.)

If we assume that the states achieved by adding a net positive charge density $\delta\ll 1$ per site to a QSL insulator can be described in terms of dilute excitations on top of a persistent background QSL, then a variety of possible ground-state phases are natural to consider: If the lowest energy charged excitations are ordinary holes, then these can either remain itinerant, forming a Fermi liquid, FL∗, (which can easily be seen to be a distinct phase from a usual Fermi liquid in that the Fermi surface only encloses an area corresponding to the density of doped holes, and so violates Luttinger’s theorem) or if they crystalize, they might form an insulating hole Wigner Crystal (hWC) with one doped hole per emergent unit cell. 111Note that the hWC has a spin-1/2 per unit cell, and so would be expected to magnetically order unless the lattice structure of the hWC is such that this higher-order lattice supports a QSL groundstate. If, on the other hand, the lowest energy charged excitations are holons, then they could condense to form a conventional superconducting state with small superfluid density $\sim\delta$, or they could crystalize. We refer to the later state as WC∗, which is distinguishable from hWC in that there are no low energy spin degrees of freedom. One could also imagine that the lowest energy charged excitations are holon pairs (or equivalently, spin-singlet hole pairs), which if they crystalize would form an insulating WC of Cooper pairs. Still more complicated phases, including supersolid phases and conducting charge density waves, could occur if a fraction of the charged excitations crystalize while others remain itinerant, or by condensing fermionic holons or bosonic spin 1/2 charge e holes.

Principal Results: We find that the spin-spin correlation function (Fig. 2) is remarkably insensitive to doping. Indeed, the spin correlation length, $\xi_{s}\lesssim 2$ lattice constants, is small compared to the circumference of the cylinders studied and to the mean separation between doped holes. The fact that they look little different than those for $\delta=0$ is consistent with viewing the system as a lightly doped QSL. The expectation value of the charge density (Figs. 4 and 5) is inhomogeneous, and the amplitude of the charge density variations is relatively insensitive to system size, implying that the ground-state in the thermodynamic limit spontaneously breaks translation symmetry. Moreover, for the most part, the charge density appears to favor a triangular lattice with one doped hole per unit cell. Given the fact that there is no significant doping-induced magnetic order, we may identify this state as a WC∗. However, the precise crystal structure of the WC∗ in the thermodynamic limit is not something we can infer with confidence, as in some cases, depending on the value of $\delta$ and the circumference of the cylinder, we find a stripe crystal rather than a triangular lattice. Given the strong evidence of crystalization, it is perhaps not surprising that we find that all superconducting correlations are extremely short-ranged (Fig. 3), with a correlation length $\xi_{SC}\lesssim 1.3$. 222It is worth noting that the phases of a putative doped QSL was considered in another recent DMRG study of the $t$-$J$ model on the square lattice with a ratio of second and first neighbor exchange interactions, $J_{2}/J_{1}\approx 1/2$ and with $\delta=1/8$. While the systems considered there were much narrower than in the present study, clear evidence was found of a tendency to form a WC of singlet hole (or holon) pairs.Dodaro et al. (2016)

Model Hamiltonian: We employ the density-matrix renormalization group (DMRG) White (1992) to investigate the ground state properties of the hole-doped Kagome antiferromagnet. The simplest model which captures the strong correlation physics on this lattice is the $t$-$J$ model (depicted in Fig.1) which is defined by the Hamiltonian

[TABLE]

where $c_{i\sigma}^{+}$ ($c_{i\sigma}$) is the electron creation (annihilation) operator with spin-$\sigma$ on site $i$. $\vec{S}_{i}$ is the spin operator and $n_{i}=\sum_{\sigma}c_{i\sigma}^{+}c_{i\sigma}$ is the electron number operator. $\langle ij\rangle$ denotes NN sites and the Hilbert space is constrained by the no-double occupancy condition, $n_{i}\leq 1$. At half-filling, i.e., $n_{i}=1$, the $t$-$J$ model reduces to the spin-$\frac{1}{2}$ antiferromagnetic Heisenberg model on the Kagome lattice.

The lattice geometry used in the DMRG simulations is depicted in Fig.1, where $\mathbb{e}_{1}$ and $\mathbb{e}_{2}$ denotes the two basis vectors of the Kagome lattice. We consider Kagome cylinders with periodic (open) boundary condition in the $\mathbb{e}_{2}$ ($\mathbb{e}_{1}$)-direction. A cylinder geometry introduced Ref.Kolley et al. (2015) (which we will refer to as YC) is used such that one of the three bond orientations is along the $\mathbb{e}_{2}$-axis. Here, we focus on cylinders with width $L_{y}$ and length $L_{x}$, where $L_{y}$ and $L_{x}$ are the numbers of unit cells (and $2L_{x}$ and $2L_{y}$ are the number of sites) along the $\mathbb{e}_{2}$ and $\mathbb{e}_{1}$ directions, respectively. Notice that the unit cells at the right boundary of the cylinder contain only two sites (A and C) in order to reduce the boundary effects due to sharp edges. Following Ref.Yan et al. (2011); Kolley et al. (2015), we refer to the cylinders as $\rm YC$-$2L_{y}$, whose total number of sites is $N=L_{y}(3L_{x}+2)=N_{u}+2L_{y}$, where $N_{u}$ denotes the number of sites inside intact unit cells. In this paper, we will focus primarily on $\rm YC$-6 and $\rm YC$-8 cylinders, i.e., $L_{y}$=3 and $4$, with $L_{x}$=12 to 24. We have also considered $\rm YC$-10 cylinders, i.e., $L_{y}$=5, and found similar results (see 4(e)). As usual, the doping level of the system away from half-filling is defined as $\delta=N_{h}/N_{u}$, where $N_{h}$ is the number of holes. Although $N_{u}\neq N$ so that the average value of $\delta$ differs slightly from $\tilde{\delta}=N_{h}/N$, deep in the bulk, i.e., relatively far from the open boundaries, $\tilde{\delta}=\delta$.

For the present study, we focus on the lightly doped case with doping level $0\leq\delta\leq 11\%$. We set $J$=1 as an energy unit and consider $t$=3. The results also hold for other $t$. We perform up to 50 sweeps and keep up to $m$=10000 states in the DMRG block with a typical truncation error $\epsilon\sim 10^{-6}$ for $\rm YC$-6 cylinders, $\epsilon\sim 10^{-5}$ for $\rm YC$-8 cylinders and $\epsilon\sim 5\times 10^{-5}$ for $\rm YC$-10 cylinders. This leads to excellent convergence for the results that we report here which have been extrapolated to $m=\infty$ limit (see Supplementary Information).

Spin-Spin correlations: To describe the magnetic properties of the ground state of the systems, we calculate the spin-spin correlation functions defined as ${\rm F(r)}=\frac{1}{L_{y}}\sum_{y=1}^{L_{y}}|\langle\mathbb{S}_{0}\cdot\mathbb{S}_{r}\rangle|$. Here $\mathbb{S}_{0}$ denotes the spin operator on the reference $A$-site in the middle of the cluster, while $\mathbb{S}_{r}$ runs over both $A$ and $B$ sites along the $\mathbb{e}_{1}$ -direction with the distance $r$ between them. At half-filling, i.e., $\delta=0$, the ground state of the system has been shown to be a QSL state Jiang et al. (2008); Yan et al. (2011); Jiang et al. (2012); Depenbrock et al. (2012); Kolley et al. (2015) with short-range spin-spin correlations. This is also confirmed by our present study where the spin-spin correlations $\rm F(r)$ for $\rm YC$-6 in Fig. 2(a) and $\rm YC$-8 cylinders in Fig. 2(b) both decay rapidly with distance $r$, and can be well fitted by an exponential function ${\rm F(r)}\sim e^{-r/\xi_{s}}$ with short correlation lengths $\xi_{s}$=1.1$\sim$1.3 lattice spacings.

Upon doping, i.e., $\delta>0$, we find that the spin-spin correlations still remain short-ranged. The spin-spin correlations $\rm F(r)$ as a function of doping level $\delta$ at different system sizes are shown in Fig.(2). For all cases, we find that the spin-spin correlations $\rm F(r)$ decay exponentially with short correlation lengths $\xi_{s}$, although $\xi_{s}$ slightly depends on the doping level $\delta$ and lattice geometry. For both $\rm YC$-6 and $\rm YC$-8 cylinders, we find that $\xi_{s}$=1$\sim$2 lattice spacings as shown in the inset of Fig.(2) (a). Our results show that the ground states of the Kagome antiferromagnet upon doping remain spin singlet with short-range spin-spin correlations, similar to the QSL state at half-filling.

Superconducting correlation: We have also investigated the possiblity of superconductivity. Since the ground state remains a spin-singlet state upon doping, we focus on spin-singlet superconductivity. The primary diagnostic of superconducting order is the pair-field correlator defined as

[TABLE]

Here, $\Delta^{\dagger}_{\alpha}(i)$ is the spin-singlet pair-field creation operator given by $\Delta_{\alpha}^{\dagger}(i)=\frac{1}{\sqrt{2}}\left(c^{\dagger}_{i\uparrow}c^{\dagger}_{i+\alpha\downarrow}-c^{\dagger}_{i\downarrow}c^{\dagger}_{i+\alpha\uparrow}\right)$, where $\alpha$ denotes the bond type (see Fig.1), i.e., a, b or c, with bond vectors defined as $\mathbb{a}=\mathbb{e}_{1}/2$, $\mathbb{c}=\mathbb{e}_{2}$ and $\mathbb{b}=(\mathbb{e}_{2}-\mathbb{e}_{1})/2$. $i_{0}$ is the index of the reference bond in the middle of the cluster, and $r$ is the distance between two bonds along the $\mathbb{e}_{1}$-direction.

In the present study, we find that $\Phi_{aa}$, $\Phi_{bb}$, $\Phi_{cc}$, $\Phi_{ab}$, $\Phi_{bc}$ and $\Phi_{ca}$ all decay exponentially for both $\rm YC$-6 and $\rm YC$-8 cylinders (see Fig.3(a) and (b)). In all cases we have studied, for large separations along the cylinder, $1\ll|r|\ll L_{x}$ where $\mathbb{r}=r{\hat{e}}_{1}$, $\Phi$ can be well expressed as $\Phi_{\alpha,\beta}(r)\sim e^{-|\mathbb{r}|/\xi^{\alpha\beta}_{sc}}$ from which we derive the values for superconducting correlation length $\xi^{\alpha\beta}_{sc}$. As shown in Fig.3(c), $\xi_{sc}=0.5\sim 1.3$ lattice spacings over all doping levels $0<\delta\leq 11\%$ we have explored. Therefore, our results suggest that there is no (quasi-) long-range superconductivity in doped QSL on the Kagome lattice.

Charge density wave order: Finally, we consider the charge density profile $n_{h}(x,y)=1-n(x,y)$, where $n(x,y)$ is the electron density on site $i=(x,y)$ with coordinations $x$ $\&$ $y$. Fig.4 shows some examples of $n_{h}(x,y)$ at different hole doping levels $\delta$ for both $\rm YC$-6 and $\rm YC$-8 cylinders. Clear CDW ordering is observed, although its pattern depends on both the lattice geometry (either $\rm YC$-6 or $\rm YC$-8) and doping level. There is unidirectional CDW order at low doping levels for $\rm YC$-6 cylinders, Fig.4(a), but this appears to be special for $\rm YC$-6 geometry. For higher doping level for $\rm YC$-6 cylinders and all doping levels for $\rm YC$-8 cylinders, the CDW order resembles a two-dimensional Winger crystal, Fig.4(b)-(d).

Approximately, the doped system can be divided into new larger emergent unit cells, i.e. each containing one of the red stripes in Fig.4(a) or one of the red-spots in Fig.4(b)-(e); the number of emergent unit cells is equal to the number of doped holes at all doping levels. This is not a crystal of Cooper pairs.

To determine whether the CDW order survives in the thermodynamic limit, we further calculate the averaged rung charge density defined by $n_{h}(x)=\frac{1}{L_{y}}\sum_{y=1}^{L_{y}}n_{h}(x,y)$. Examples of $n_{h}(x)$ at different $\delta$ for both $\rm YC$-6 and $\rm YC$-8 cylinders are plotted in Fig.5. Here, the existence of long-range CDW order in the ground state can be determiend by fitting the amplitude $A_{cdw}$ of the oscillation of the charge density $n_{h}(x)$ and extrapolating the value to $L_{x}=\infty$. To minimize the boundary effect, we have removed a few points from both ends in the fitting process. Examples of the extrapolation are given in Fig. 5 for system sizes $L_{x}=12\sim 40$. The observed finite amplitude $A_{cdw}$ in the thermodynamic limit establishes the presence of long-range CDW order.

Discussion: In light of our observations, it is worth asking if there is an intuitive reason that the holons crystalize, rather than forming one of the possible quantum fluid states that have been proposed. It is already clear from previous numerical studies that the QSL in the Kagome system is nearly degenerate with a number of possible valence-bond-crystalline phases. It is thus natural to imagine that the holon is a highly structured particle, surrounded by a “polaronic” cloud of valence-bond-crystal like correlations.

In Fig.S3 in the Supplemental Section, we show that two doped holes in cylinders of moderate length ($L_{x}$=12 & 16) induce a strong and extended pattern of valence bond crystalline order in their neighborhoods. The most obvious consequence of this is that the holon effective mass is strongly renormalized (increased). Moreover, the induced valence bond order implies the existence of moderate range effective interactions between pairs of holons, which if they are repulsive can naturally lead to crystallization. Note that in Fig.S2, we show that the spin-gap in the presence of two doped holes (extrapolate to the $L_{x}\to\infty$ limit) is of the same order, although probably smaller than in the undoped ladder; this further corroborates our identification of this as a two holon state.

Since the Kagome antiferromagnet, i.e., $H_{J}$ in Eq.(1), has been shown to be a realistic model to describe the celebrated Herbertsmithite Hering and Reuther (2016); Jiang et al. (2017), our results may be directly relevant to the real material upon doping. Consistent with the present results, a recent experimental study has reported the absence of superconductivity in doped Kagome systems either with doped electrons Kelly et al. (2016) or with doped holes Julien et al. (2013). While we have only studied finite cylinders, it is plausible that the results are representative of the thermodynamic limit, given the fact that the size of the cylinders, including both width $L_{y}$ and length $L_{x}$, are much larger than both the spin-spin and superconducting correlation lengths.

Acknowledgement: We would like to thank Senthil Todadri, Young Lee, John Tranquada, Masaki Oshikawa, Arun Paramekanti, Roderich Moessner, Steve White, Zheng-Yu Weng, Hong Yao, Shenxiu Liu and especially Leon Balents for insightful discussions, and John Dodaro for making Fig.S3 (b) and (d). This work was supported by the Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, under Contract DE-AC02-76SF00515. Parts of the computing for this project was performed on the Sherlock cluster.

Appendix A A: Extrapolated Correlation Funcation

In Fig. S1(a), we show the spin-spin correlation function $\rm F(r)$ for $\rm YC$-8 cylinders with $L_{x}$=18 at $\delta=\frac{1}{27}$ hole doping concentration, keeping $m=2187\sim 8000$ states. The navy blue line with star symbol shows the second-order polynomial extrapolation of the correlation function to $\epsilon\rightarrow 0$, i.e., $m=\infty$. The inset shows the correlation length $\xi_{s}$, fitted using the exponential function $\rm F(r)\sim e^{-r/\xi_{s}}$, as a function of truncation error $\epsilon$. The red line is a second-order fit of the data which gives the accurate $\xi_{s}$ in the $\epsilon\rightarrow 0$ limit.

In Fig. S1(b), we show the superconducting correlation function for the same system. For simplicity, we define the averaged superconducting correlation function as $\rm\Phi(r)=[\Phi_{aa}(r)+\Phi_{bb}(r)+\Phi_{cc}(r)]/3$. The navy blue line with star symbols show the second-order polynomial extrapolation of $\rm\Phi(r)$ to $\epsilon\rightarrow 0$, and the inset shows the correlation length $\xi_{sc}$ as a function of $\epsilon$. The red line denotes the second-order fit.

Appendix B B: Spin Gap

In Fig. S2, we shown the spin gap, defined as $\Delta_{s}(N_{h})=E_{0}(N_{h},S_{z}=1)-E_{0}(N_{h},S_{z}=0)$ where $E_{0}(N_{h},S_{z})$ denotes the ground state energy for system with total spin $S_{z}$ doped with $N_{h}$ holes, for $\rm YC$-6 and $\rm YC$-8 cylinders at half-filling ($N_{h}$=0) and doped with 2 holes ($N_{h}$=2). At half-filling, $\Delta_{s}$ is finite for both cylinders, which is consistent with previous studies Jiang et al. (2008); Yan et al. (2011); Jiang et al. (2012); Depenbrock et al. (2012). When doped with 2 holes, $\Delta_{s}$ also appears to be finite for $\rm YC$-6 cylinder whose value is close to half-filling. However, for $\rm YC$-8 cylinder, we find that $\Delta_{s}$ appears to be much smaller than half-filling.

Appendix C C: Valence Bond Order

In order to study the correlations that are associated with the dressed holon, we have studied the properties of two doped holes in a system of a moderate size so that the charge density associated with the two holons are well localized near opposite ends of the cylinder, far enough from each other so that they hardly overlap. To look for local valence-bond crystalline correlations, we measure the difference in the valence bond strength defined as $B_{d}(i,j)=\langle\vec{S}_{i}\cdot\vec{S}_{j}\rangle_{N_{h}=2}$ - $\langle\vec{S}_{i}\cdot\vec{S}_{j}\rangle_{N_{h}=0}$, where $i$ and $j$ are nearest neighbor sites and the subscripts indicate that the expectation values are taken in the ground-states with $N_{h}$=2 and $N_{h}=0$ doped holes, respectively. We show both the expectation value of the charge density (Figs. S3 (a) & (c)) and $B_{d}(x,y)$ (Figs. S3 (b) & (d)) for a $\rm YC$-6 system with $L_{x}$=16 and a $\rm YC$-8 system with $L_{x}$=12. Clearly, the holon is an extended object with a diameter of roughly 5 lattice constants in which strong local valence bond crystalline correlations are manifest.