Evidence of pair-density wave in doping Kitaev spin liquid on the honeycomb lattice

The pair-density wave (PDW) is a superconducting (SC) state in which the order parameter varies periodically in space in such a way that its spatial average vanishes.Berg et al. (2009); Agterberg et al. (2020) The first example of PDW is the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) stateFulde and Ferrell (1964); Larkin and Ovchinnikov (1965) when a Zeeman magnetic field, $H$, is applied to a s-wave superconductor so that the Fermi surface is spin-split. The SC order has a wave vector $Q\sim\mu_{B}H/E_{F}$ which is typically very small. The LO version of this state is accompanied by an induced magnetization density wave and a charge density wave (CDW) with ordering wavevectors $K=2Q$. Intense interest in a somewhat different sort of PDW state has emerged due to recent discoveries in underdoped cuprate superconductors, where a direct observation of PDW has been made experimentally via local Cooper pair tunneling and scanning tunneling microscopy in ${\rm Bi_{2}Sr_{2}Ca_{2}O_{8+x}}$Hamidian et al. (2016); Ruan et al. (2018); Edkins et al. (2019) as well as the dynamical inter-layer decoupling observed in $1/8$ hole doped ${\rm{La}_{2}BaCuO_{4}}$.Berg et al. (2007); Agterberg and Tsunetsugu (2008) While similar in having oscillatory SC order and associated $K=2Q$ CDW order, this PDW is conjectured to be stable in zero magnetic field (zero net magnetization), have an ordering vector that is independent of $H$ (at least for small or vanishing $H$), and moreover can either have no associated magnetic order, or possibly have spin density wave order (SDW) with the same ordering vector $Q$.

Although much is known about the properties of the PDW stateBerg et al. (2007); Agterberg and Tsunetsugu (2008); Berg et al. (2009); Lee (2014); Jian et al. (2020), there are very few microscopic models which are shown to have PDW ground states. These include the one-dimensional (1D) Kondo-Heisenberg model with 1D electron gas coupled to a spin chainBerg et al. (2010) and the extended two-leg Hubbard-Heisenberg modelJaefari and Fradkin (2012). The evidence of PDW is also observed in the t-J model with four-spin ring exchange on a 4-leg triangular latticeXu et al. (2019) and an extended Hubbard model with a staggered spin-dependent magnetic flux per plaquette on a 3-leg triangular lattice.Venderley and Kim (2019) However, there is no evidence of PDW ordering found in more standard models even with second neighbor interactions.Dodaro et al. (2017); Jiang et al. (2018); Jiang and Devereaux (2019); Jiang et al. (2020)

Theoretically, it has been proposed that superconductivity can also emerge in doping quantum spin liquids (QSLs), which are exotic phases of matter that exhibit a variety of novel features associated with their topological character.Balents (2010); Savary and Balents (2016); Broholm et al. (2020) The QSLs 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 et al. (1989); Lee et al. (2006, 2007); Fradkin et al. (2015) Indeed, recent numerical studies using density-matrix renormalization group (DMRG) have provided strong evidences that lightly doping the QSL and chiral spin liquid on the triangular lattice will naturally give rise to nematic $d$-waveJiang (2019) and $d\pm id$-wave superconductivityJiang and Jiang (2020), respectively. Although dominant SC correlations have been observed in both systems, there is still no evidence for PDW ordering.

In this paper, we define a $t$-$J$-like extension of the Kitaev model on the honeycomb lattice (Fig.1) so as to address the question of whether SC emerges upon light doping. The Kitaev model, i.e., $J_{K}$ term in Eq.(1), is exactly solvable and has a gapless spin liquid as its exact ground state.Kitaev (2006) Moreover, it has potential experimental realizations in magnets with strong spin-orbit coupling such as ${\rm Na_{2}IrO_{3}}$ and ${\rm\alpha}$-${\rm RuCl_{3}}$.Jackeli and Khaliullin (2009); Jiang et al. (2011); Singh et al. (2012); Hwan Chun et al. (2015); Kim and Kee (2016); Banerjee et al. (2017); Trebst (2017); Takagi et al. (2019); Hickey and Trebst (2019); Jiang et al. (2019); Yokoi et al. (2020); Chern et al. (2020) This provides us a unique theoretical opportunity to test the physics of doping QSLs and may also give us some hints for understanding the mechanism of high temperature superconductiivity in the cuprates. Theoretically, doping Kitaev spin liquid (KSL) has been studied and distinct metallic states were proposed.You et al. (2012); Hyart et al. (2012); Schmidt et al. (2018); Mei (2012) These include the $p$-wave superconductivityYou et al. (2012); Hyart et al. (2012), topological superconductivitySchmidt et al. (2018), and Fermi liquid state.Mei (2012) However, controlled results of the sort that can be obtained using density-matrix renormalization group (DMRG) are still lacking concerning the phase(s) that arise upon doping the KSL.

Principal results: In the present paper, we study the lightly doped Kitaev model in Eq.(1) on the honeycomb lattice using DMRGWhite (1992). Based on DMRG calculations on three-leg cylinders, we find that upon light-doping the KSL state, the system exhibits power-law CDW correlations at long distances corresponding to a local pattern of partially-filled charge stripes. For three-leg cylinders, the wavelength of CDW, i.e., the spacings between two adjacent charge strips in the $\mathbf{e}_{1}$ direction, is $\lambda_{c}=a_{0}/3\delta$, where $a_{0}$ is the length of unit cell. This corresponds to an ordering wavevector $K=3\pi\delta/a_{0}$ with two thirds of a doped hole per CDW unit cell.

We find that the even-parity SC correlations are the most pronounced SC correlations, far dominant compared with the odd-parity SC or topological superconductivity.You et al. (2012); Hyart et al. (2012); Schmidt et al. (2018) Moreover, the dominant pairing channel oscillates in sign as a function of distance which is consistent with the striped PDW.Agterberg et al. (2020) Its wavelength $\lambda_{sc}=2a_{0}/3\delta$ is the same as that of the spin density wave (SDW) correlations $\lambda_{s}=2a_{0}/3\delta$ and two times of that of the CDW $\lambda_{c}=a_{0}/3\delta$. Correspondingly, the SC ordering wavevector $Q=3\pi\delta/2$ is half of the ordering wavevector $K=3\pi\delta/a_{0}$ of the CDW. Although quasi-long-range SC correlations are not seen in the range of parameters we have studied, short-range correlations are fairly strong with the corresponding correlation length $\xi_{sc}\geq 3a_{0}$. This is comparable with the cylindrical width and is notably larger than that of doping the Kagome QSL.Jiang et al. (2017) Similar with doping the QSL on the triangular latticeJiang (2019), we find that the pairing symmetry of SC correlations is also consistent with $d$-wave. To the best of our knowledge, this is the first observation of PDW in doping the KSL.

Model and Method: We employ DMRGWhite (1992) to investigate the ground state properties of the hole-doped Kitaev model on the honeycomb lattice defined by the Hamiltonian

Here $c^{\dagger}_{i\sigma}$($c_{i\sigma}$) is the electron creation (annihilation) operator with spin-$\sigma$ on site $i=(x_{i},y_{i})$, $S^{\gamma}_{i}$ is the $\gamma$-component of the $S=1/2$ spin operator on site $i$, where $\gamma=x,y,z$ labels the three different links of the hexagonal lattice as illustrated in Fig.1. $\langle ij\rangle$ denotes nearest-neighbor (NN) sites and the Hilbert space is constrained by the no-double occupancy condition, $n_{i}\leq 1$, where $n_{i}=\sum_{\sigma}c^{\dagger}_{i\sigma}c_{i\sigma}$ is the electron number operator. At half-filling, i.e., $n_{i}=1$, Eq.(1) reduces to the Kitaev model, which is known to have a gapless spin liquid ground state that can be gapped out into a non-Abelian topological phase by certain time-reversal symmetry perturbations.Kitaev (2006); Zhu et al. (2018)

The lattice geometry used in our simulations is depicted in Fig.1, where $\mathbf{e}_{1}=(\sqrt{3},0)$ and $\mathbf{e}_{2}=(1/2,\sqrt{3}/2)$ denote the two basis vectors. We consider honeycomb cylinders with periodic and open boundary conditions in the $\mathbf{e}_{2}$ and $\mathbf{e}_{1}$ directions, respectively. We focus on cylinders of width $L_{y}$ and length $L_{x}$, where $L_{y}$ and $L_{x}$ are the number of unit cells ($L_{y}$ and $\tilde{L}_{x}=2L_{x}$ are the number of sites) along the $\mathbf{e}_{2}$ and $\mathbf{e}_{1}$ directions, respectively. The total number of sites is $N=2\times L_{y}\times L_{x}$. The hole doping concentration is defined as $\delta=N_{h}/N$, where $N_{h}$ is the number of doped holes. We set $J_{K}=1$ as an energy unit and consider $t=3$. In this paper, we focus primarily on three-leg cylinders, i.e., $L_{y}=3a_{0}$, of length $L_{x}=12a_{0}\sim 48a_{0}$ (i.e., $\tilde{L}_{x}=24\sim 96$), at doping levels $\delta=1/12$ and $1/9$. We find similar results on four-leg cylinders which are provided in the Supplemental Materials (SM). We keep up to $m=8000$ number of states in each DMRG block with a typical truncation error $\epsilon\sim 10^{-6}$. This leads to excellent convergence for our results when extrapolated to $m=\infty$ limit.

Pair-density wave: To test the possibility of superconductivity, we have calculated the equal-time SC pair-pair correlations. A diagnostic of the SC order is the SC pair-pair correlation function, defined as

where $\Delta^{\dagger}_{\alpha}(x,y)=\frac{1}{\sqrt{2}}[\hat{c}^{\dagger}_{(x,y),\uparrow}\hat{c}^{\dagger}_{(x,y)+\alpha,\downarrow}-\hat{c}^{\dagger}_{(x,y),\downarrow}\hat{c}^{\dagger}_{(x,y)+\alpha,\uparrow}]$ is the even-parity SC pair-field creation operator, where the bond orientations are labelled as $\alpha=x,y,z$ (Fig.1). ($x_{0},y$) is the reference bond with $x_{0}\sim\tilde{L}_{x}/4$ to minimize the boundary effect and $r$ is the distance between two bonds in the $\mathbf{e}_{1}$ direction. We have also calculated the odd-parity SC correlations to test the possibility of $p$-wave or topological superconductivity. However, we find that they are much weaker than the even-parity SC correlations (see SM for details), suggesting that $p$-wave or topological superconductivity is unlikely. Therefore, we will focus on the even-parity SC correlation in this paper.

Fig.2 shows the SC pair-pair function $\Phi_{yy}(r)$ for doping $\delta=1/12$ and $1/9$. The SC correlation shows clear spatial oscillation for both doping levels, which can be well fitted by $\Phi_{yy}(r)\sim f(r)*\phi_{yy}(r)$ for a large region of $r$ as we will discuss below. Here $f(r)$ is the envelope function and $\phi_{yy}(r)$ is a spatial oscillatory function. At long distances, the envelope function $f(r)$ is consistent with an exponential decay, i.e., $f(r)\sim e^{-r/\xi_{sc}}$, as shown in Fig.2a-b. The extracted correlation length is $\xi_{sc}\geq 3a_{0}$, which is comparable with the cylindrical width $L_{y}=3a_{0}$. Alternatively, the SC correlation at long distances can also be fitted by a power law (see SM for details), i.e., $f(r)\sim r^{-K_{sc}}$, with an exponent $K_{sc}>4$ for both doping levels $\delta=1/12$ and $1/9$. We have also measured other types of SC correlations $\Phi_{\alpha\beta}(r)$ which are provided in the SM. While they are slightly weaker than $\Phi_{yy}(r)$ due to the broken symmetry induced by the cylindrical geometry, they have very similar decaying behavior with $\Phi_{yy}(r)$. Although the SC correlations can be fitted by either functions, it is clear that the SC susceptibility will not diverge in the thermodynamic limit.

The spatial oscillation of the SC correlation $\Phi_{yy}(r)$ is characterized by the normalized correlation $\phi_{yy}(r)=\Phi_{yy}(r)/f(r)$, as shown in Fig.2c. It is clear that $\phi_{yy}(r)$ varies periodically in real space and can be well fitted by $\phi_{yy}(r)\sim{\rm sin}(Qr+\theta)$ for both cases. This is consistent with the PDW state with vanishing spatial average of $\phi_{yy}(r)$. The ordering wavevector of the PDW is $Q=3\pi\delta/2a_{0}$ as indicated by the peak position of the Fourier transformation $\phi_{yy}(q)$ of $\phi_{yy}(r)$ (Fig.2d) and $\theta$ is a fitting phase factor. The corresponding wavelength is $\lambda_{sc}=2a_{0}/3\delta$, which is $\lambda_{sc}=8a_{0}$ for $\delta=1/12$ and $\lambda_{sc}=6a_{0}$ for $\delta=1/9$. As we will see below that the relation $\lambda_{sc}=\lambda_{s}=2\lambda_{c}$ can be clearly seen in our results as expected from that of striped PDW. Here, $\lambda_{c}$ and $\lambda_{s}$ are the wavelengths of the CDW and SDW, respectively. According to Ginzburg-Landau theory, the development of $K=2Q$ charge oscillation corresponds to the cubic term $\rho_{K}\Delta^{*}_{Q}\Delta_{-Q}$ in the free energy, where $\rho_{K}$ and $\Delta_{Q}$ are the charge density and PDW order parameters with corresponding momenta $K$ and $Q$, respectively. This is distinct from the CDW modulated superconductivity with term $\rho_{Q}\Delta^{*}_{0}\Delta_{-Q}$ with coexisting dominant uniform $\Delta_{0}$ and secondary stripe pairing $\Delta_{Q}$.

To identify the pairing symmetry, we have calculated the SC correlations using both real-valued and complex-valued DMRG simulations. We first rule out the $d\pm id$-wave symmetry as we find that both the wavefunction and SC correlations are real while their imaginary parts are zero. We have further analyzed the relative phase of different SC correlations, e.g., $\Phi_{xy}(r)/\Phi_{yy}(r)$ in the insets of Fig.2a-b, where a clear out of phase can be observed as $\Phi_{xy}(r)/\Phi_{yy}(r)<0$. Therefore, we conclude that the pairing symmetry of the PDW is consistent with $d$-wave.

Charge density wave: In addition to SC correlations, we have also measured the charge density profiles to describe the charge density properties of the system. The rung charge density $n(x)=\sum_{y=1}^{L_{y}}n(x,y)/L_{y}$ is shown in Fig.3, where $x$ is the rung index of the cylinder and $n(x,y)$ is the local charge density on site $i=(x,y)$. It is clear that the charge density varies periodically in real space along the $\mathbf{e}_{1}$ direction with the wavelength $\lambda_{c}=a_{0}/3\delta$, i.e., $\lambda_{c}=4a_{0}$ and $\lambda_{c}=3a_{0}$ for $\delta=1/12$ and $\delta=1/9$, respectively. Therefore, there are two thirds of a doped hole in each CDW unit cell. Moreover, it is clear that the relation $\lambda_{sc}=2\lambda_{c}$ holds for both $\delta=1/12$ and $\delta=1/9$, which is consistent the striped PDW state. Different with SC correlations, we find power-law decay of the charge density correlation at long distances. The Luttinger exponent $K_{c}$ of the power-law decay can be extracted by fitting the charge density oscillation (Friedel oscillation) induced by the open boundaries of the cylinderWhite et al. (2002)

Here $x$ is the distance in the $\mathbf{e}_{1}$ direction from the open boundary, $n_{0}$ the average density and $K$ is the ordering wavevector. $\delta n$ and $\theta$ are the model dependent constants. Examples of the fitting using Eq.(3) on the a sub-lattice are given in Fig.3a for both doping $\delta=1/12$ and $\delta=1/9$, where four data points near the boundary are removed to minimize boundary effect for a more reliable fit. The extracted exponent $K_{c}$ is given in Table 1. It is clear that $K_{c}<1$ which demonstrates the dominance of the charge density correlations.

Alternatively, we can estimate $K_{c}$ from the amplitude $A_{cdw}(L_{x})$ of the charge density modulation.Jiang and Devereaux (2019) For a given cylinder of length $L_{x}$, the CDW amplitude $A_{cdw}(L_{x})$ can be obtained by fitting the central-half region of the charge density profile $n(x)$.Jiang et al. (2018); Jiang and Devereaux (2019); Jiang et al. (2020) For quasi-long-range charge order, the amplitude should follow $A_{cdw}(L_{x})\propto L_{x}^{-K_{c}/2}$ with similar $K_{c}$ with that obtained from Friedel oscillation in Eq.(3). This is indeed the case as shown in Fig.3c-d, where the exponent $K_{c}$ is given in Fig.3b and Table 1.

Spin-spin and single-particle correlations: To describe the magnetic properties of the system, we have further calculated the spin-spin correlation function $F_{\gamma}(r)=\langle S^{\gamma}_{i_{0}}S^{\gamma}_{i_{0}+r}\rangle$. Here $r$ is the distance between two sites in the $\mathbf{e}_{1}$ direction and $\gamma=x,y,z$. $i_{0}=(x_{0},y)$ is the reference site with $x_{0}\sim\tilde{L}_{x}/4$. For the pure Kitaev model without doping, it is known that $F_{\gamma}(r)$ is nonzero only for NN sites.Kitaev (2006); Baskaran et al. (2007) Upon doping, the $Z_{2}$ flux on each hexagonal plaquette is no longer a conserved quantity, and the spin-spin correlation functions become nonzero even at long distance. For both doping levels $\delta=1/12$ and $\delta=1/9$ on three-leg cylinders, we find that both $F_{x}(r)$ and $F_{z}(r)$ decay exponentially as $F(r)\sim e^{-r/\xi_{s}}$, where the spin-spin correlation length is given in Table.1. On the contrary, $F_{y}(r)$ appears to be long-range ordered as shown in Fig.4. This is unexpected but allowed theoretically since there is no continuous spin symmetry in the system. Moreover, we find that all spin-spin correlations show clear spatial oscillation with a wavelength $\lambda_{s}$ that is the same as that of the SC correlation, i.e., $\lambda_{s}=\lambda_{sc}$, which is consistent with the striped PDW. However, our results suggest that the long-range correlation $F_{y}(r)$ is special to three-leg cylinder. On the contrary, $F_{y}(r)$ decays exponentially on the wider four-leg cylinders (Fig.4b), which is similar with both $F_{x}(r)$ and $F_{z}(r)$. However, the evidences of the striped PDW are robust for both three and four-leg cylinders with similar sign-changing SC correlations (see SM for details).

In addition to superconductivity, we have also measured the single-particle Green function $G_{\sigma}(r)=\langle c^{\dagger}_{i}c_{i+r}\rangle$ to test the possibility of Fermi liquid state.Mei (2012) It is clear that the single-particle Green function decays exponentially at long distances as $G_{\sigma}(r)\sim e^{-r/\xi_{G}}$ (see SM for details). Here $\xi_{G}$ is the correlation length which is summarized in Table 1. As a result, our study suggests that the Fermi liquid state is unlikely in the lightly hole-doped KSL.

Central charge: Our results suggest that the ground state of the lightly hole-doped KSL has quasi-long-range CDW correlation with gapless charge mode. To show this, we have calculated the von Neumann entanglement entropy $S(x)=-{\rm Tr}\left[\rho_{x}\ln\rho_{x}\right]$, where $\rho_{x}$ is the reduced density matrix of subsystem with length $x$. For 1D critical system described by conformal field theory, it has been establishedCalabrese and Cardy (2004); Fagotti and Calabrese (2011) that for an open finite system of length $L_{x}$,

Here $c$ is central charge, $k_{F}$ denotes the Fermi momentum, $a$ and $\tilde{c}$ are model dependent parameters. Fig.5a shows an example of $S(x)$ for three-leg cylinder at doping level $\delta=1/9$. Here we have calculated $S(x)$ by dividing the system into two parts with smooth boundary through both $x$ and $y$ bonds. The extracted central charge $c$ is given in Fig.5 as a function of $L_{x}$ at doping levels $\delta=1/9$ and $1/12$. It is clear that the central charge quickly converges to $c=1$ with the increase of $L_{x}$, although it deviates notably from $c=1$ on short cylinders due to finite-size effect. Therefore, our results show that there is one gapless charge mode which is consistent with the quasi-long-range CDW correlations. We have obtained similar results on four-leg cylinders which are provided in the SM.

Summary and discussion: Admittedly, the DMRG calcualtions are carried out on finite length cylinders. However, based on the results we have obtained we conjecture that the exact ground state for an infinite long three-leg cylinder has the following properties: 1) There is a single gapless charge mode and a gap (which produces exponentially falling correlations) for all spin carrying excitations. 2) Long-range SDW order with the ordered moment oriented in the $y$ direction - that is along the circumference of the cylinder - and an ordering vector $Q=3\pi\delta/2$; the connected spin correlations fall exponentially with a correlation length $\xi_{s}\sim a_{0}$. 3) There are power-law CDW correlations with an exponent $K_{c}\sim 2/3$ and an ordering vector $K_{c}=2Q$. 4) There are strong even-parity d-wave-like PDW correlations with an ordering vector $Q$ which fall either exponentially with a correlation length $\xi_{sc}\sim 3a_{0}$ or possibly with a high power law $K_{sc}\geq 4$. 5) All other forms of SC correlations - those corresponding to odd-parity pairing or spatially uniform even-parity pairing - are much weaker in comparison with the PDW correlations.

There are many aspects of these observations that are surprising, and will need to be understood theoretically. The PDW correlations are sufficiently short-ranged that one would infer (based on any reasonable conjecture concerning their time dependence) that the corresponding PDW susceptibility would be finite. Thus, at present, we can only conclude that the present results are suggestive of a possible PDW ordered state for the 2D (infinite leg) version of this model. However, it is notable that the favored forms of order are remarkably reminiscent of those conjectured to be present in the cuprate high temperature superconductor, LBCO, where direct evidence exists of stripe SDW order with wave vector $Q\approx 4\pi/\delta a$ and CDW order with ordering vector $K=Q/2$, and indirect evidence has been adduced for PDW order with ordering vector $Q$.

In this paper, we primarily focus on the lightly doped Kitaev model, it will be interesting to study the higher doping case as well as the extend Kitaev model with further neighbor hopping which is shown to be essential to enhance the superconductivity on the square lattice.Jiang and Devereaux (2019); Jiang et al. (2020) As other terms such as the Heisenberg interaction and spin-orbit couplings are also present in real materials such as ${\rm Na_{2}IrO_{3}}$ and $\alpha$-${\rm RuCl_{3}}$,Jackeli and Khaliullin (2009); Jiang et al. (2011); Singh et al. (2012); Hwan Chun et al. (2015); Kim and Kee (2016); Banerjee et al. (2017); Trebst (2017); Takagi et al. (2019); Hickey and Trebst (2019); Jiang et al. (2019); Yokoi et al. (2020) it will be interesting to study these systems as well.

Acknowledgments: We would like to thank Senthil Todadri for insightful discussions, and especially Steve Kivelson for insightful discussion, invaluable suggestions and generous help to understand the results and improve the manuscript. 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.