Helical Majorana fermions in d_{x^2-y^2} + i d_{xy}-wave topological superconductivity of doped correlated quantum spin Hall insulators

Large Hubbard U limit of the Kane-Mele model on a zigzag ribbon of honeycomb lattice near half-filling is studied via a renormalized mean-field theory. The ground state exhibits time-reversal symmetry (TRS) breaking d + i d'-wave superconductivity. At large spin-orbit coupling, the Z2 phase with non-trivial spin Chern number in the pure Kane-Mele model is persistent into the TRS broken state (called spin-Chern phase), and has two pairs of counter-propagating helical Majorana modes at the edges. As the spin-orbit coupling is reduced, the system undergoes a topological quantum phase transition from the spin-Chern to chiral superconducting states. Possible relevance of our results to adatom-doped graphene and irridate compounds is discussed.

Introduction. Searching for topological states of quantum matters constitutes one of the central and fundamental issues in condensed matter systems. The growing interest in topological insulators (TIs), which support gapless edge (or surface) states protected by timereversal symmetry (TRS) while the bulk remains insulating [1,2], is one prime example. Of particular interest are topological superconductors which support gapless selfconjugate, charge-neutral fermionic quasi-particle excitations [3]. These excitations which reflect non-trivial topological bulk properties are localized at the edges, known as Majorana fermions (MFs).
Much effort has been put in searching for signatures of Majorana fermions in solid state materials. Onedimensional semiconductor nano-wires with strong SO coupling under a magnetic field proximity to a s-wave superconductor have been proposed theoretically to host MF at both ends of the wire [4][5][6], and also studied experimentally [7][8][9][10][11]. Similar ideas have been proposed in 2D systems where chiral MFs exist at the edges of spintriplet, p−wave (odd-parity) superconductors [12][13][14].
While realization of the above systems relies on TRS breaking by the Zeeman field, time-reversal invariant topological superconductors (TRITOPs) [15][16][17][18] has recently been proposed to host two time-reversal pairs of helical MFs at edges in repulsively interacting SO coupled nano-wire proximity to either a s-wave [6,16] or a d−wave [19] superconductor at each end of the wire. Proposals to realize TRITOPs in 2D systems include the spin-triplet p x ± i p y superconductors [15], the bi-layer Rashba system [20], and in exciton condensates [21].
In this letter, we suggest a novel mechanism for realizing helical Majorana fermions in 2D topological superconductors with spin-singlet and TRS-breaking pairing gap-by directly doping correlated 2D quantum spin Hall insulators (QSHIs or 2D TIs) [22,23]. A paradigmatic model for QSHIs is the Kane-Mele (KM) model [23,24], which shows a non-trivial Z 2 topological (or spin Chern) number and supports helical topological edge states protected by TRS. The half-filled KM model with strong electron correlations is in the Mott-insulating (MI) phase [25], while superconductivity appears upon doping. Attractive candidates to realize correlated KM model include: graphene with enhanced SO coupling (∼ 20meV ) by doping with heavy adatoms, such as indium or thallium [26], Iridium-based honeycomb compounds X 2 IrO 3 (X = N a or Li) with strong SO coupling and electron correlations [25,27].
Via renormalized mean-field theory (RMFT) approach [28], we find the spin-singlet TRS breaking d x 2 −y 2 +i d xywave superconductivity to appear at the ground state where the chiral edge states have been shown to occur [29]. Surprisingly, for sufficiently large SO coupling, instead of chiral edge states, we find zero-energy helical MFs to appear at each ribbon edge. We further show that the system exhibits a non-trivial Z 2 topological invariant, which supports helical MFs at edges despite the presence of TRS breaking superconductivity. This seemingly un-expected feature comes as a result of persistence of spin-Chern phase in the pure Kane-Mele model in the superconducting state at finite dopings. The competition between SO coupling and chiral d−wave superconductivity leads to a novel spin-Chern to chiral topological quantum phase transition. Nearest-neighbor and next-nearestneighbor lattice vectors areδa=1,2,3, ai=1,2 with an unit length ofā, a, respectively. We set a = 1 here. The gray shaded region represents for the super-unit-cell of the zigzag ribbon, which repeats itself along x−axis. Mele t-J (KM-tJ) model is given by [23]: where α =↑, ↓ stands for the spin index, i, j and i, j refer to the nearest-neighbor (NN) and next-nearestneighbor (NNN) sites, respectively (see Fig. 1 (a)). Here, ν ij = 1 for i, j ∈ A and ν ij = −1 for i, j ∈ B in the SO coupling; S i refers to the electron spin operator on site i, defined as: S i = 1/2 α,β=↑,↓ c †α i σ αβ c β i , n i = α c †α i c α i is the electron density operator, and the anti-ferromagnetic spin-exchange couping J ∼ 4t 2 U can be derived via the second-order perturbation from the Kane-Mele Hubbard model in the limit of a strong onsite Coulomb repulsion U t [30]. The H J term has been known to favor the spinsinglet pairing. To address superconductivity of the model, we apply RMFT based on Gutzwiller projected single-occupancy constraint due to the proximity of the Mott insulating ground states, known to well describe the ground state of d−wave cuprate superconductors [28]. As a result, the hopping t term acquires a reduction factor g t at a finite doping δ: t → tg t with g t = 2δ/(1 + δ), while the spin-exchange J term gets enhances by a factor g s : J → g s J with g s = 4/(1 + δ) 2 . The spin-exchange J term within RMFT reads: , and αβ = i σ y αβ . Based on the C 6v symmetry of the underlying lattice, the pairing symmetry of ∆ δa=1,2,3 may take the following forms [31]: Fig. 1) [32]. The Fourier transformed pairing gap ∆ k for a periodic 2D lattice reads: ∆ k = a=1,2,3 ∆ δa e i k· δa .

The mean-field Hamiltonian
The Hamiltonian Eq. (.2) possesses both the Particlehole (PH) symmetry: with τ x being the Pauli matrix on particle-hole space where K stands for complex conjugation as well as sublattice symmetry : The matrixĥ k , describing the KM model, shows TRS: However, M k breaks the TRS for d + i d superconducting order parameter: ∆ d+id , and The mean-field free energy reads: Ns with E k > 0 being positive eigenvalues and N s the number of sites. We diagonalized the mean-field Hamiltonian H k on a finite-sized zigzag ribbon with N s = N/2 zigzag chains and N = 56 is set as the total number of sites along y−axis throughout the paper.
Results. The mean-field variables are solved selfconsistently by minimizing the free energy both for a periodic lattice and a finite-sized ribbon. On a 2D periodic lattice, the results as a function of doping are shown in Fig. 1(b) and (c). Compared to the TRS extended s−wave, we find d x 2 −y 2 +i d xy −wave pairing is the ground state [32]. Same pairing symmetry has been reported in superconducting phase of the doped graphene in the absence of the spin-orbit coupling [29,31,33,34], which was argued to support two co-propagating chiral edge states at low energies with a non-trivial topological winding number N T KN N = ±2 [35]. The superconducting transition temperature T c is estimated as T c ∼ g t ∆ 0 .
On a finite-sized zigzag ribbon and at a generic doping, the Bogoliubov quasi-particle dispersion shows four doubly-degenerate bulk bands (due to the S z symmetry of our Hamiltonian) grouped in two pairs (see Fig.  2(a)); it satisfies the particle-hole symmetry with 2π periodicity: 3 ) (see Fig. 2 [32]. Surprisingly, in the regime of a strong SO coupling and weak pairing (∆ d+id √ 3t SO ), we find the low energy excitations of our model support helical MFs at edges instead of chiral edge states as expected for a chiral d−wave superconductor. On a finite-sized zigzag ribbon, we find two Dirac-dispersed lines intersecting at momenta k M F x ∼ 2π/3, 4π/3 where the Bogoliubov quasi-particle excitation energy vanishes, E(k M F x ) = 0 (see Fig. 2(a) and Fig. 3(a)) [36]. Near each of these gapless points, two pairs of two-fold degenerate states are generated via intersecting the two Dirac lines by a constant energy at two momentum points k M F x1(x2) , denoted as: are co-propagating (see Fig.  3 , exhibit an exponential decay from both edges into the bulk, where R/L refers to the right/left moving state, and u(i),ũ(i),ṽ(i), v(i) are the corresponding matrix elements.
These features are clearly different from the copropagating chiral edge states realized either in the chiral d-wave superconductivity in doped graphene or by proximity of a s−wave superconductor to a quantum anamolous Hall insulators [38]. Instead, the edge states we find fit well to the helical MFs described by the linearly-dispersed Hamiltonian defined by the Bogoliubov quasi-particle operators γ R(L)τ k as: with α = L, R refers to the Bogoliubov quasi-particle destruction operator defined by the coherence factors u τ k,ī ,ũ τ k,ī ,ṽ τ k,ī , v τ k,ī , corresponding to the right-moving quasi-particle with "pseudo-spin" τ =↑ (↓ ) , and the summation over repeated site indexī = 1, · · · N is implied; similarly for γ Lτ k . The pair of the degenerate wavefunctions Ψ can be expressed as Ψ R(L),↑(↓) (ī), formed by the coherence factors: ; similarly for the other doublet Ψ R(L),↓(↑) (ī). It is clear from Fig. 3 (b) that the edge states (Ψ R(L),↑(↓) , Ψ L(R),↓(↑) ) (as well as the Bogoliubov operators (γ R↑(L↓) k , γ L↓(R↑) k )) form pairs (see pink (blue) curve in Fig. 3(b) for |Ψ R,↑ (i)| 2 (|Ψ L,↓ (i)| 2 )). Furthermore, these Bogoliubov operators with linear dispersion satisfy γ −k (−E) = γ † k (E) withk ≡ k − π via PH symmetry (see top and bottom panels of Fig. 3(b)). Hence, they can be regarded as examples of helical MFs at edges [15]; the MF zero-modes occur at k = k M F x (ork = 0) where γ ατ k=0 = γ †ατ k=0 . These helical MFs are protected by an additional P-H symmetry: γ −k (E) = γ † k (E) (see Fig. 3(b)). Our seemingly unexpected results have roots in the competition between TRS SO couping and TRS breaking chiral d−wave superconductivity. We find the TRS protected Z 2 QSH insulating phase of the pure undoped Kane-Mele model persists up to a finite doping and a finite pairing gap in the chiral d−wave superconducting phase.
To gain more understanding, we compute the Z 2 topological invariant (spin Chern number) N W = 1 2 (C n − Cn) with C n = 1 2πi k∈BZ F 12 (k) where the integral is done in the first Brillouin zone (BZ), the field strength F 12 (k) and the associated Berry's connection A(k) are defined as: F 12 (k) ≡ ∂ 1 A 2 (k)−∂ 2 A 1 (k) and A µ = n(k)|∂ µ |n(k) with |n(k) being the normalized eigenvector of the nth band [39]. Cn is defined similarly withn being the corresponding degenerate band associated with the nth band via the transformation S z =↑→↓. In the strong SO coupling regime, t SO ∆ 0 and for a sizable range of doping around 1/2-filling, the spin-Chern phase prevails and we find C n = −Cn or N W = ±1 [37], same result as that for the Z 2 TRITOPs supporting helical MFs at edges [1,15,40,41]. In this limit, our system is well approximated by the effective spin-singlet p x ± i p y superconductivity near the two Dirac points K ± . This can be seen when re-expressing the superconducting pairings in terms of the electron operators ψ ±,k which diagonaliz the tight-binding KM Hamiltonian [32,37]. In this basis, the intra-band pairing ∆ −− k ψ †↑ −,k ψ †↓ −,−k dominates at ground state (see Fig. 2 (c) and Ref. 37). Near K ± points with q ± = K ± + (±q x + q y ), we find ∆ −− q± ∼ ±∆ 0 (q y ± i q x ), resembling the case of a TRI-TOPs. In the opposite limit for t SO ∆ 0 or sufficiently large doping where the chiral d-wave pairing dominates, however, we recover the chiral superconductivity: N W = 0 and C n = Cn = 1, equivalent to the case of doped graphene [29,31]. A novel spin-Chern-to-chiral topological quantum phase transition is found as ∆ 0 /t SO or µ/t SO is varied [42]. The generic phase diagram by tuning µ (in a non-self-consistent way) at a fixed ∆ 0 /t SO is shown in Fig.4 (a). For ∆ 0 t SO , we find the critical values of µ being at µ = ±µ c ∼ ± √ 3t SO . The bulk band gap closes only at the phase transition, while it remains open in either phase [42][43][44] (see Fig. 4(b)). The similar persistence of spin-Chern phase in a TRS breaking magnetic field has been observed experimentally in Ref. 45. The spin-Chern phase of our system belongs to class D topological superconductors, distinct from the TRITOPs [19,46,47].

Conclusions.
We demonstrate here for the first time a 2D spin-singlet topological superconductor with non-trivial spin Chern number in doped correlated Kane-Mele model, which supports helical counter-propagating Majorana zero modes despite the d+i d superconducting pairing gap breaks TRS. This seemingly unexpected feature comes as a result of persistence of spin-Chern phase of the pure Kane-Mele model in the superconducting state upon doping. As T → 0, distinct differential conductance spectrum for each pair of Majorana zero mode through differential Andreev conductance in the normalmetal/superconducting (NS) junction is expected [29,31].
We thank Y. Oreg  in the bulk. This can be understood as the electronic density of states (DOS) of the un-doped pure KM ribbon is enhanced at the edges due to the presence of the helical edge states. As a result, the superconducting pairing strength is enhanced at edges. However, the inhomogeneity of superconducting gap at the edges concerns only with the bulk states at a sizable negative energy. Since we focus in this work on the Majorana fermions near zero energy, far away from the superconducting states at edges, we therefore ignore this spatial inhomogeneity and consider only the superconducting gap with an uniform magnitude within our RMFT approach.

B. Doping evolution of the Bogoliubov excitation spectrum
Here, we provide the doping evolution of the Bogoliubov excitation spectrum as shown in Fig. 2. The spin-Chern phase with helical Majorana modes survive at low dopings. For J/t = 0.1, t SO /t = 0.5, it survives for 0 < δ < 0.2. Above δ = 0.2, the superconducting gap vanishes and the system is in the normal state. On the other hand, the chiral phase will prevail in the opposite regime t SO ∆ 0 for low dopings. However, the chiral superconductivity has been addressed extensively (see Refs. 29,35 in the main text). We will leave this issue out of the focus of our present work.

III. Z2 TOPOLOGICAL INVARIANT
In this section, we provide details on obtaining the Z 2 topological number (or the spin Chern number) in the doped Kane-Mele t-J model following the approach in Ref. 2. First, the TKNN number, a topological (Chern) number, for the n−th band is defined as 1 : where A(k) ≡ n(k)|∇ k |n(k) , |n(k) is the Bloch state in the n-th band up to a normalization coefficient, which satisfies the Schroedinger equation H(k)|n(k) = E n (k)|n(k) . The existence of A(k) implies that there is a non-trivial magnetic field across the first Brillouin Zone (FBZ). The FBZ of the honeycomb lattice we choose for the integration Eq. (III.1) is illustrated in Fig. 3, the parallelogram (area enclosed by red lines) spanned by vectors q 1 and q 2 in momentum space.
In order to numerically calculate the TKNN number covered by the entire FBZ (area enclosed by the solid red lines) in Fig. 3, we first discretize it into N − 1 × N − 1 small flakes of parallelograms such as the smaller parallelogram abcd (enclosed by solid black lines) in Fig.  3. The position of a discretized lattice point (that is the vertices of the small parallelograms ) can be specified by a vector k i,j = ( k i1 , k j2 ) , for (i, j) ∈ [1 , N ] with its components where q ix and q iy are the components of the vectors q 1 and q 2 with unit lattice spacing |ā| = 1 , they are given by The phase difference ∆θ 1 (k ij ) of Ψ n (k) ≡ k|n(k) travelling from a lattice point k i,j to k i+1,j can be computed to be The subscript "1" of ∆θ 1 and U 1 (k ij ) represents that Ψ n travels along the direction of q 1 . Sum over ∆θ of the four sides of the small parallelogram abcd in Fig. 3 can one obtain the total phase shift of the Bloch wavefunction on the n-th band after circling around that small parallelogram. Compare with Eq. (III.1), we define the lattice U (1) gauge field strength F n ij by  For the same band we can further define a topological number, so-called the spin Chern number, as: where Cn is the TKNN number of the other degenerate eigenvector of the same band. Numerically, we obtain Cn = −C n , thus N W is either 1 or −1. We checked numerically that for the parameters used in Fig. 3 of the main text |N W | approaches to 1.
In the other extreme limit, ∆ t SO , we found that N T KN N ≡ C n + Cn = ±2, and N W = 0, as expected for d + i d superconductivity in doped graphene 3 . Further studies suggest the existence of a quantum phase transition at a critical ratio of ∆/t SO separating the spin-Chern phase with N W = ±1 from the chiral (TKNN) phase with N T KN N = ±2 (see Ref. 3).

IV. EFFECTIVE TRS SUPERCONDUCTIVITY NEAR THE DIRAC POINTS
In this section, we derive the low-energy effective TRS superconducting pairing of our model near Dirac points.
We re-express the d + i d superconducting order parameter in terms of the electron operators ψ ±,k which diagonalizes the tight-binding KM Hamiltonian with the corresponding eigenvalue E ± = ±( 2 k + γ 2 k − µ) (see Ref. 5):