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

We study the effects of doping the Kitaev model on the honeycomb lattice where the spins interact via the bond-directional interaction $J_K$, which is known to have a quantum spin liquid as its exact ground state. The effect of hole doping is studied within the $t$-$J_K$ model on a three-leg cylinder using density-matrix renormalization group. Upon light doping, we find that the ground state of the system has quasi-long-range charge-density-wave correlations but short-range single-particle correlations. The dominant pairing channel is the even-parity superconducting pair-pair correlations with $d$-wave-like symmetry, which oscillate in sign as a function of separation with a period equal to that of the spin-density wave and two times the charge-density wave. Although these correlations fall rapidly (possibly exponentially) at long distances, this is never-the-less the first example where a pair-density wave is the strongest SC order on a bipartite lattice. Our results may be relevant to ${\rm Na_2IrO_3}$ and $\alpha$-${\rm RuCl_3}$ upon doping.

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. [1,2] The first example of PDW is the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state [3,4] 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 ∼ µ 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 Bi 2 Sr 2 Ca 2 O 8+x [5][6][7] as well as the dynamical inter-layer decoupling observed in 1/8 hole doped La 2 BaCuO 4 . [8,9] 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 state [1,[8][9][10][11], 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 chain [12] and the extended two-leg Hubbard-Heisenberg model [13]. The evidence of PDW is also observed in the t-J model with four-spin ring exchange on a 4-leg triangular lattice [14] and an extended Hubbard model with a staggered spin-dependent magnetic flux per plaquette on a 3-leg triangular lattice. [15] However, there is no evidence of PDW ordering found in more standard models even with second neighbor interactions. [16][17][18][19] 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. [20][21][22] The QSLs can be viewed as insulating phases with preexisting electron pairs such that upon light doping they might automatically yield high temperature superconductivity. [23][24][25][26][27][28][29][30] 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-wave [31] and d ± id-wave superconductivity [32], 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 arXiv:2008.03858v1 [cond-mat.str-el] 10 Aug 2020 ground state. [33] Moreover, it has potential experimental realizations in magnets with strong spin-orbit coupling such as Na 2 IrO 3 and α-RuCl 3 . [34][35][36][37][38][39][40][41][42][43][44][45] 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. [46][47][48][49] These include the p-wave superconductivity [46,47], topological superconductivity [48], and Fermi liquid state. [49] 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 DMRG [50]. 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 e 1 direction, is λ c = a 0 /3δ, where a 0 is the length of unit cell. This corresponds to an ordering wavevector K = 3πδ/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. [46][47][48] Moreover, the dominant pairing channel oscillates in sign as a function of distance which is consistent with the striped PDW. [2] Its wavelength λ sc = 2a 0 /3δ is the same as that of the spin density wave (SDW) correlations λ s = 2a 0 /3δ and two times of that of the CDW λ c = a 0 /3δ. Correspondingly, the SC ordering wavevector Q = 3πδ/2 is half of the ordering wavevector K = 3πδ/a 0 of the CDW. Although quasilong-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 ξ sc ≥ 3a 0 . This is comparable with the cylindrical width and is notably larger than that of doping the Kagome QSL. [51] Similar with doping the QSL on the triangular lattice [31], 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 DMRG [50] to investigate the ground state properties of the hole-doped Kitaev model on the honeycomb lattice defined by the Hamiltonian Here c † iσ (c iσ ) is the electron creation (annihilation) op- erator with spin-σ on site i = (x i , y i ), S γ i is the γcomponent of the S = 1/2 spin operator on site i, where γ = x, y, z labels the three different links of the hexagonal lattice as illustrated in Fig.1. ij denotes nearest-neighbor (NN) sites and the Hilbert space is constrained by the no-double occupancy condition, n i ≤ 1, where n i = σ c † iσ c iσ 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. [33,52] The lattice geometry used in our simulations is depicted in Fig.1, where e 1 = ( √ 3, 0) and e 2 = (1/2, √ 3/2) denote the two basis vectors. We consider honeycomb cylinders with periodic and open boundary conditions in the e 2 and 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 andL x = 2L x are the number of sites) along the e 2 and e 1 directions, respectively. The total number of sites is N = 2 × L y × L x . The hole doping concentration is defined as δ = 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 ∼ 48a 0 (i.e.,L x = 24 ∼ 96), at doping levels δ = 1/12 and 1/9. We find similar results on four-leg cylinders which are provided in the Supplemen-tal Materials (SM). We keep up to m = 8000 number of states in each DMRG block with a typical truncation error ∼ 10 −6 . This leads to excellent convergence for our results when extrapolated to m = ∞ limit.
Pair-density wave: To test the possibility of superconductivity, we have calculated the equal-time SC pairpair correlations. A diagnostic of the SC order is the SC pair-pair correlation function, defined as is the even-parity SC pair-field creation operator, where the bond orientations are labelled as α = x, y, z ( Fig.1). (x 0 , y) is the reference bond with x 0 ∼L x /4 to minimize the boundary effect and r is the distance between two bonds in the e 1 direction. We have also calculated the odd-parity SC correlations to test the possibility of pwave 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 Φ yy (r) for doping δ = 1/12 and 1/9. The SC correlation shows clear spatial oscillation for both doping levels, which can be well fitted by Φ yy (r) ∼ f (r) * φ yy (r) for a large region of r as we will discuss below. Here f (r) is the envelope function and φ 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) ∼ e −r/ξsc , as shown in Fig.2a-b. The extracted correlation length is ξ sc ≥ 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) ∼ r −Ksc , with an exponent K sc > 4 for both doping levels δ = 1/12 and 1/9. We have also measured other types of SC correlations Φ αβ (r) which are provided in the SM. While they are slightly weaker than Φ yy (r) due to the broken symmetry induced by the cylindrical geometry, they have very similar decaying behavior with Φ 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 Φ yy (r) is characterized by the normalized correlation φ yy (r) = Φ yy (r)/f (r), as shown in Fig.2c. It is clear that φ yy (r) varies periodically in real space and can be well fitted by φ yy (r) ∼ sin(Qr + θ) for both cases. This is consistent with the PDW state with vanishing spatial average of φ yy (r). The ordering wavevector of the PDW is Q = 3πδ/2a 0 as indicated by the peak position of the Fourier transformation φ yy (q) of φ yy (r) (Fig.2d) and θ is a fitting phase factor. The corresponding wavelength is λ sc = 2a 0 /3δ, which is λ sc = 8a 0 for δ = 1/12 and λ sc = 6a 0 for δ = 1/9. As we will see below that the relation λ sc = λ s = 2λ c can be clearly seen in our results as expected from that of striped PDW. Here, λ c and λ 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 ρ K ∆ * Q ∆ −Q in the free energy, where ρ K and ∆ 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 ρ Q ∆ * 0 ∆ −Q with coexisting dominant uniform ∆ 0 and secondary stripe pairing ∆ Q .
To identify the pairing symmetry, we have calculated the SC correlations using both real-valued and complexvalued DMRG simulations. We first rule out the d ± idwave 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., Φ xy (r)/Φ yy (r) in the insets of Fig.2a-b, where a clear out of phase can be observed as Φ xy (r)/Φ 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) = Ly y=1 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 e 1 direction with the wavelength λ c = a 0 /3δ, i.e., λ c = 4a 0 and λ c = 3a 0 for δ = 1/12 and δ = 1/9, respectively. Therefore, there are two thirds of a doped hole in each CDW unit cell. Moreover, it is clear that the relation λ sc = 2λ c holds for both δ = 1/12 and δ = 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 powerlaw decay can be extracted by fitting the charge density oscillation (Friedel oscillation) induced by the open boundaries of the cylinder [53] n(x) = n 0 + δn * cos(Kx + θ)x −Kc/2 .
Here x is the distance in the e 1 direction from the open boundary, n 0 the average density and K is the ordering wavevector. δn and θ 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 δ = 1/12 and δ = 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 I. 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. [18] For a given cylinder of length L x , the CDW amplitude A cdw (L x ) can be obtained by fitting the centralhalf region of the charge density profile n(x). [17][18][19] For quasi-long-range charge order, the amplitude should follow A cdw (L x ) ∝ L −Kc/2 x 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 I.
Spin-spin and single-particle correlations: To describe the magnetic properties of the system, we have further calculated the spin-spin correlation function F γ (r) = S γ i0 S γ i0+r . Here r is the distance between two sites in the e 1 direction and γ = x, y, z. i 0 = (x 0 , y) is the reference site with x 0 ∼L x /4. For the pure Kitaev model without doping, it is known that F γ (r) is nonzero only for NN sites. [33,54] 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 δ = 1/12 and δ = 1/9 on three-leg cylinders, we find that both F x (r) and F z (r) decay exponentially as F (r) ∼ e −r/ξs , where the spin-spin correlation length is given in Table.I. 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 λ s that is the same as that of the SC correlation, i.e., λ s = λ 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 σ (r) = c † i c i+r to test the possibility of Fermi liquid state. [49] It is clear that the single-particle Green function decays exponentially at long distances as G σ (r) ∼ e −r/ξ G (see SM for details). Here ξ G is the correlation length which is summarized in Table I. 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) = −Tr [ρ x ln ρ x ], where ρ x is the reduced density matrix of subsystem with length x. For 1D critical system described by conformal field theory, it has been established [55,56] that for an open finite system of length L x , Here c is central charge, k F denotes the Fermi momentum, a andc are model dependent parameters.  shows an example of S(x) for three-leg cylinder at doping level δ = 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 δ = 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πδ/2; the connected spin correlations fall exponentially with a correlation length ξ s ∼ a 0 . 3) There are power-law CDW correlations with an exponent K c ∼ 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 ξ sc ∼ 3a 0 or possibly with a high power law K sc ≥ 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 ≈ 4π/δ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. [18,19] As other terms such as the Heisenberg interaction and spin-orbit couplings are also present in real materials such as Na 2 IrO 3 and α-RuCl 3 , [34][35][36][37][38][39][40][41][42][43][44] it will be interesting to study these systems as well.

Results for four-leg cylinders
We provide more results for four-leg cylinder, where we have observed similar features of the striped PDW discussed in the main text. Similar with three-leg cylinders, the density profiles on four-leg cylinders also exhibit similar spatial oscillation (Fig.A1a), with the period λ c = a 0 /4δ corresponding to an ordering vector K = 4πδ/a 0 , where a 0 is the length of unit cell.
The SC pair-pair correlation function Φ yy (r) on fourleg cylinder of length L x = 24 and L x = 36 at δ = 1/12 doping are shown in Fig.A1b. Similar with threeleg cylinders, the SC pair-pair correlation Φ yy (r) decays exponentially whose envelope can be well fitted by the function f (r) ∼ e −r/ξsc with a correlation length ξ sc = 1.37 (8) and ξ sc = 1.51 (3) for L x = 24 and L x = 36, respectively. The spatial oscillation of Φ yy (r) is characterized by the normalized pair-pair correlation φ yy = Φ yy (r)/f (x) shown in Fig.A1c. The ordering wavevector of Φ yy (r) is Q = 2πδ/a 0 as indicated by the peak position of the Fourier transformed φ yy (q) in Fig.A1d. The relationship Q = K/2 between the SC ordering vector Q and CDW ordering vector K also holds true for four-leg cylinders, as expected for the striped PDW state.
On four-leg cylinders, all the spin-spin correlations F γ (r), where γ = x, y, z, decay exponentially at long distances, e.g., F z (r) as shown in Fig.A1e. The envelope of the strongest component F y (r) (see Fig.4b) can be well fitted by the exponential function f (r) ∼ e −r/ξs with a correlation length ξ s = 2.76(8) from length L x = 36 cylinder. The single-particle Green's function also decays exponentially at long distances as G σ (r) ∼ e −r/ξ G (Fig.A1f) with the correlation length ξ G = 2.48 (5).

Power-law fitting of SC pair-pair correlation
Here we double-check the alternative power-law fitting for the SC pair-pair correlation. Fig.A3 shows the same SC pair-pair correlations Φ yy (r), i.e., Fig.2 in the main text, in double-logarithm scale. At long distances, it is quite clear that Φ yy (r) deviates from the power-law fitting labelled by the linear dashed linesf (r) ∼ r −Ksc with the exponent K sc = 4.0(4) and 4.3(2) for δ = 1/9 and δ = 1/12 dopings, respectively. This deviation is also reflected in the vanishing normalized correlatioñ φ yy (r) = Φ yy (r)/f (r) at long distance and the rounded peak ofφ yy (q) located at the momentum q = 3πδ/2a0 shown in Fig.A3c and d.

All types of even parity pair-pair correlations
In Fig.A4, we list all types of even-parity SC pair-pair correlations Φ αβ (r) on three-leg cylinders of length L x = 48 at doping δ = 1/9 and 1/12. It is clear that they all decay exponentially, similar with Φ yy (r) discussed in main text, although with slightly different amplitude at long distances. This could be attributed to the fact that the cylinder geometry breaks the C 3 rotational symmetry of the honeycomb lattice.

Single-particle Green function
To check the possibility of other metallic state such as the Fermi liquid state, We have also measured the singleparticle Green function G σ (r), as shown in Fig.A5. It is clear that G σ (r) ∼ e −r/ξ G decays exponentially, which is inconsistent with the Fermi liquid state. For instance, the correlation lengths at doping δ = 1/9 and 1/12 are ξ G = 1.7(1) and 1.71(3) on three-leg cylinders of length L x = 48. The correlation length ξ G in the long-cylinder limit after finite-size scaling using L x = 18 ∼ 48 is slightly longer, which is given in Table I of the main text.