# Itinerant Ferromagnetism and Superconductivity in Doped Bilayer Silicene

###### Abstract

We study the electronic instabilities of doped bilayer silicene using the random phase approximation. In contrast to the singlet superconductivity at the low doping region, we find that the system is an itinerant ferromagnet in the narrow doping regions around the Van Hove singularities, and a triplet superconductor in the vicinity of these regions. Adding a weak Kane-Mele spin-orbit coupling to the system further singles out the time-reversal invariant equal-spin helical pairing as the leading instability. The triplet pairing identified here is driven by the ferromagnetic fluctuations, which become strong and enhance the superconducting critical temperature remarkably near the phase boundaries between ferromagnetism and superconductivity.

###### pacs:

75.10.Lp, 74.20.Rp, 74.25.Dw^{†}

^{†}thanks: yangfan_blg@bit.edu.cn

^{†}

^{†}thanks: ygyao@bit.edu.cn

## I Introduction

Magnetism and unconventional superconductivity (SC) as well as the intimate interplay between them have been the focuses of condensed matter physics for decades due to their rich physics and important applications. Among these subjects, the realizations of itinerant ferromagnetism (FM) and triplet SC are of particular importance in recent years. In general, the triplet SC triplet (), which is connected with topological SC tsc1 (); tsc2 () and becomes hot topic recently, is believed to be driven by ferromagnetic spin fluctuations near the FM order. However, the realization of itinerant FM from the Stoner criterion stoner () usually requires finite and, most of the time, strong electron interactions ston1 (); ston2 (); ston3 (); ston4 (); ston5 (); ston6 (); ston7 (); ston8 (); ston9 (), which is hard to deal with in the weak coupling perturbative approaches. One way to overcome this difficulty is through the introducing of divergent density of states (DOS) at the Van Hove (VH) singularities of the system, which can induce these instabilities without strong electronic interactions. It’s proposed recently that, for a system with its Fermi surface (FS) doped to time-reversal (TR) variant VH saddle points, weak repulsive electron interactions can usually drive itinerant FM and triplet SC VH ().

On another front, as the Si-based counterpart of graphene, silicene has been synthesized recently syn1 (); syn2 (); syn3 (); syn4 (); syn5 (), with experimental evidence showing possible SC in the doped case gap (). Furthermore, bilayer silicene (BLS) has also been available synbi (), with the energetically most favored stacking way between its two layers identified by first-principles calculations bilayer (). Based on the metallic band structure of undoped BLS, the antiferromagnetism (AFM) and the chiral SC tuned by strain have been proposed bilayer (). This intriguing result motivates us to further investigate the electronic instabilities in doped BLS, specifically focusing on the VH doping levels since the divergent DOS there favors the occurrence of electronic instabilities. Paying our attention to VH doping, we notice that, in VH-doped monolayer graphene whose VH saddle points locate at TR invariant momenta, the chiral spin density wave or the chiral pairing has been proposed graph1 (); graph2 (); graph3 (); graph4 (); graphsc (). Similar results have also been found in monolayer silicene silicene (). In contrast, the interesting property of the VH singularities here in BLS lies in that the VH saddle points locate at TR variant momenta. For such VH singularities, the study based on the renormalization group theory has pointed out the possibility of the formation of itinerant FM and triplet SC VH ().

In this paper, we perform the calculations based on the random phase approximation (RPA) to investigate possible electronic instabilities of doped BLS. The main results of our calculations are as follows. In addition to the SC occurring at low doping levels, the itinerant FM and the triplet SC emerge as the leading instabilities of the system in the narrow doping regions around the VH singularities and the vicinity of these regions respectively. In the presence of a weak Kane-Mele spin-orbit coupling (SOC), the helical pairing becomes the leading instability of the system. The emergence of the FM and the triplet SC results from the large DOS and the strong ferromagnetic correlation near the VH singularities. Near the critical doping level separating the FM and triplet SC, the strong ferromagnetic fluctuations will greatly enhance the superconducting critical temperature, which provides possibility to realize this triplet pairing state at experimentally accessible temperatures.

The rest of this paper is organized as follows. In Sec. II, we describe the Hubbard model of BLS, as well as the RPA approach. In Sec. III, we calculate the susceptibilities of the system, and demonstrate the itinerant FM occurring around the VH singularities. In Sec. IV, we study the superconducting pairing symmetries for different doping levels, and propose that the pairing dominates over the one in the vicinity of the ferromagnetic regions. Finally in Sec. V, a conclusion will be reached after discussions on the experimental detection of the novel pairing state proposed here.

## Ii Model and Approach

The lattice structure of BLSbilayer () is shown in Fig. 1(a), which belongs to the point group. While sublattice of the upper silicene layer couples vertically to sublattice of the lower layer with a bond-length Å, the two sublattices and within the same layer () couples to each other with a bond-length Å. Approximately equal bond lengths, together with the bond-angle between the two bonds describes an orbital hybridization more like the type than the planar one. This lattice structure leads to a strong interlayer coupling, and the resulting strong bonding-antibonding splitting between orbitals and pushes them far away from the Fermi level. Thus the low energy subspace formed by orbitals and will take responsibility for the main physics of the system bilayer (). This feature of BLS is obviously different from that of bilayer graphene.

According to Ref. bilayer (), the low energy physics of BLS near the FS can be described by the following 4-band Hubbard model of the system:

(1) |

Here , (), and denote the spin, orbital, and unit cell indices respectively, and is the 4-band tight-binding (TB) Hamiltonian in the basis . The explicit expression of the TB Hamiltonian reads bilayer ()

(2) |

Here with () being the nearest-neighbor vector, eV is the effective on-site energy difference between atoms and , the hopping integrals eV, eV, eV, and eV. Since the basis is mainly composed of the orbital of silicon u1 (), we set eV as a rough estimate of the Hubbard interaction.

The band structure for the above TB Hamiltonian (2) is shown in Fig. 1(b). One feature of the band structure is the meV overlap between the valence and conduction bands near the K-points. For the undoped case, this overlap causes six pairs of small electron and hole pockets around and near the K-points respectively as shown in Fig. 2(a). The undoped system is thus intrinsically metallic and can enter a superconducting state bilayer (). When the system is doped, regardless of electron or hole doping case, the shape of the FS grows gradually from separated electron and hole pockets first to six big merged pockets around the K-points (2(b) and 2(d)), which finally connect to one another at the VH saddle points, causing the Lifshits transition of the FS (2(c) and 2(e)). Defining the doping level by where is the number of electrons per site, we find the doping level for the VH singularities are for electron doping and for hole doping respectively. The Fermi level for these VH dopings are marked in Fig. 1(b), where the flatness of the band structure on the FS leads to the logarithmically divergent DOS near the Fermi level, as shown in Fig. 2(f). We shall focus on these VH dopings in the following study because the divergent DOS around there urges the formation of itinerant FM and the resulting strong ferromagnetic fluctuations in the vicinity of the FM region will induce high-temperature triplet SC.

A special feature of the VH singularities of this material lies in that its VH saddle points locate on the M- axes rather than at the TR invariant M-points as in bilayer graphene as well as monolayer graphene or silicene. Such TR variant VH saddle points are named as “type-II” VH saddle points in Ref. VH (), in contrast to the TR invariant “type-I” VH saddle points. The “type-II” VH singularities is special in that it allows for the formation of triplet SC. If the FS of a system contains TR invariant “type-I” VH saddle points, the triplet pairing will not be energetically favored because its odd parity gap function will have nodes at these TR invariant VH momenta, which is no good for the energy gain. On the contrary, the TR variant “type-II” VH saddle points of BLS locating on the M- axes provide the possibility for the system to enter the triplet pairing state.

To study the electron instabilities of this system represented by the Hubbard model (1), we adopt the standard multi-orbital RPA approach bilayer (); rpa1 (); rpa2 (); rpa3 (); rpa4 (); rpa5 (). We first define and calculate the bare susceptibility tensor . After that, the renormalized spin(s) or charge(c) susceptibilities are obtained in the RPA level. For each doping level, there will be a critical interaction strength . For repulsive , the renormalized spin susceptibility diverges, implying the formation of long-range magnetic order. For , through exchanging the spin or charge fluctuations, we obtain the effective pairing potential . Solving the linearized gap equation for as an eigenvalue problem, we obtain the leading pairing gap function as the eigenvector corresponding to the largest eigenvalue.

## Iii Itinerant Ferromagnetism

The bare () susceptibility tensor of the model is defined as

(3) |

Here denotes the thermal average for , denotes the time-ordered product, and are the sublattice indices. Fourier transformed to the imaginary frequency space, the bare susceptibility can be expressed by the following explicit formulism,

(4) |

Here is the Matsubara frequency, are the band indices, is the Fermi distribution function, and are the eigenvalue and eigenvector of the TB Hamiltonian (2). The Hermitian static susceptibility matrix is defined as . For each , the largest eigenvalue of this matrix represents the static susceptibility of the system in the strongest channel, and the corresponding eigenvector describes the pattern of the dominant intrinsic spin correlation in a unit cell of the system.

In Figs. 3(a)-3(e), we show the -space distributions of the zero-temperature static susceptibility for different doping levels, which reveal the doping evolution of the static susceptibility. In particular, when the doping level changes gradually from zero to the VH doping, regardless of electron or hole doping case, the momenta of the maximum susceptibility evolve from the -point (3(a)) first to the points around it (3(b) and 3(d)), and finally back to the -point again (3(c) and 3(e)). Such a doping evolution of the susceptibility originates from the evolution of the FS mentioned before, and indicates that the intra-sublattice spin correlation of the system changes gradually with doping from ferromagnetic first to antiferromagnetic, and finally back to ferromagnetic again.

From the eigenvector corresponding to the largest eigenvalue of , we finds that the spin correlation within a unit cell is ferromagnetic-like (see Fig. 3(f)) near the VH doping levels, and antiferromagnetic-like (see Fig. 3(g)) near zero doping. Therefore, although the intra-sublattice spin correlations in both the VH doped and undoped systems are ferromagnetic, the inter-sublattice spin correlations in the former and latter cases are ferromagnetic and antiferromagnetic respectively.

When the interaction is turned on, we define the charge () and spin () susceptibilities for the model as

(5) | ||||

(6) | ||||

(7) | ||||

(8) |

where , , are spin indices. For nonmagnetic states, we have . For , we further have .

In the standard RPA approach bilayer (); rpa1 (); rpa2 (); rpa3 (); rpa4 (); rpa5 (), the charge (spin) susceptibility for the model is given by

(9) |

where is a matrix, whose only four nonzero elements are () bilayer (). Clearly, the repulsive Hubbard interaction here suppresses and enhances . When the interaction parameter is weak enough, the RPA works well since all eigenvalues of the denominator matrix in Eq. (9) are positive and hence the matrix itself has an inverse. However, if exceeds a critical value at which the lowest eigenvalue of touches zero, the renormalized spin susceptibility would diverge, which implies the formation of long-range magnetic order.

The doping dependence of the critical interaction strength is shown in Fig. 4. The most obvious feature of Fig. 4 is that drops abruptly to zero near the VH singularities due to the divergent DOS there. For eV adopted in our calculation, the drops below in narrow doping regions around the VH singularities, which will lead to magnetic long-range order. From the ferromagnetic correlation near the VH singularities revealed by shown in Figs. 3(c), 3(e), and 3(f), we conclude that long-range itinerant FM will emerge for in these narrow doping regions.

## Iv Triplet Sc

Away from the above introduced narrow doping region for itinerant FM, the interaction strength is smaller than the critical value . Then through exchanging short-range spin or charge fluctuations between a Cooper pair, exotic superconducting states will emerge in the system as shown in Fig. 5.

More specifically, we consider the scattering of a Cooper pair from the state in the -th () band to the state in the -th () band via exchanging spin or charge fluctuations. This scattering process leads to the following effective interaction vertex rpa5 ():

(10) |

Here, for the singlet channel, we have

(11) |

and for the triplet channel, we have

(12) |

From the effective interaction vertex (IV), we obtained the following linearized gap equation rpa4 () near the superconducting critical temperature :

(13) |

Here the integration is along various FS patches labelled by or , is the Fermi velocity and is the component of along the FS. Solving this gap equation as an eigenvalue problem, one obtains each pairing eigenvalue and the corresponding normalized eigenvector as the relative pairing gap function. The leading pairing symmetry is determined by the corresponding to the largest . The critical temperature is determined by through , where the cutoff energy scales with the low energy bandwidth.

According to its point group, we study the possible pairing symmetries of the system including , , , and -wave ones. The doping dependences of the largest pairing eigenvalues for these pairing symmetries are shown in Fig. 5. At the low doping region, the doubly degenerate and singlet pairings serve as the leading pairing symmetries, consistent with our previous results for the undoped case bilayer (). The gap function of the and symmetries are symmetric and antisymmetric about the -axis and -axis respectively as shown in Figs. 6(a) and 6(b). Below , these two degenerate pairing states will further mix to form the fully-gapped (abbreviated as ) pairings to lower the energy. Physically, the singlet pairing is mediated by the antiferromagnetic spin fluctuations suggested by Figs. 3(a), 3(b), 3(d), and 3(g). More importantly, in the vicinity of the narrow doping region for itinerant FM around the VH singularities, our RPA results identify the doubly degenerate and triplet pairings as the leading pairing symmetries. The gap function of the () symmetry is symmetric about the -axis and antisymmetric about the -axis respectively, with gap nodes on the -axis, as shown in Fig. 6(c)(6(d)). The emergence of the triplet and pairings near the VH singularities is the physical consequence of the strong ferromagnetic spin fluctuations there, as revealed by Figs. 3(c), 3(e), and 3(f).

Since the two -wave pairing symmetries are degenerate, they will probably mix to lower the energy below the critical temperature . To determine this mixture, we set , where and denote the normalized gap functions of corresponding symmetries. Then the mixing coefficients , , and are determined by the minimization of the total mean-field energy

(14) |

Here the chemical potential is determined by the constraint of the average electron number in the superconducting state. Our energy minimization gives and , which leads to the fully-gapped (abbreviated as ) SC. This mixture of the two -wave pairings satisfies the requirement that the gap nodes should avoid the FS to lower the energy. Note that there can be three different components of this triplet pairing with different quantum numbers of the Cooper pair, which are , and . In the absence of SOC, the three spin components are degenerate.

To lift up the degeneracy among the three different spin components of the triplet pairing, we add the following Kane-Mele SOC term to Hamiltonian (1):

(15) |

Here with and being unit vectors along the two bonds that connect next-nearest-neighbors and on the same layer. Such a SOC term lifts up the degeneracy between the component and the components. Our RPA calculations (see the Appendix for the details) yield that the equal-spin helical , pairing wins over the chiral pairing by a small split proportional to , as shown in the insets of Fig. 5 for meV. Such a helical triplet pairing is TR invariant weak topological SC.

It is interesting to note that the pairing eigenvalue of the SC diverges in the doping region near the phase boundaries between FM and SC, due to the divergently strong ferromagnetic spin fluctuations in that region. Although the divergence of is an artifact in the RPA caused by ignorance of the renormalization of the single-particle Green’s function, it is possible that the strong FM fluctuations in the critical regions considerably push up the superconducting critical temperature , which might be experimentally accessible. Taking into account that the doping level is hard to control in practice, we can instead apply tunable strain to the system to change the hopping parameters and the band structure bilayer (). As a result, the VH doping levels and the phase boundaries between FM and SC will shift so that a given doping level can access the phase boundaries to produce the high-temperature triplet SC.

## V Discussion and Conclusion

The unconventional triplet SC proposed here can be detected by various experiments. First of all, as an unconventional superconducting state with the phase of its pairing gap function changing on the FS, the pairing state should show no Hebel-Slichter peak upon the superconducting phase transition in the NMR relaxation rate Hebel (). Secondly, in this triplet pairing state, the Knight shift should not obviously change below the Knight (). To further identify the phase structure of this pairing experimentally, we can fabricate a slice of BLS into a hexagon, and use a dc SQUID to detect the relative phase among different directions in the system squid (). In particular, determined by the -wave symmetry, the phase difference between the opposite (adjacent) edges of the hexagon should be ().

Although the -wave SC is unconventional, the mixing of the and pairings into the complex one leads to a fully-gapped superconducting state, which looks similar to the conventional -wave one in many aspects. For example, near zero temperature, the specific heat, the penetration depth, and the NMR relaxation rate in both fully-gapped pairing states decay exponentially with temperature. What’s more, the STM spectra of both fully-gapped superconducting states should exhibit U-shaped curves. However, all these expected experimental results can be changed by a uniaxial strain applied on the system. More specifically, the mixing proposed here is based on the degeneracy between the and pairing states, and the degeneracy itself originates from the point group of BLS bilayer (). Thus, by applying a uniaxial strain to break the symmetry of the system, we can eliminate the mixing, and leave a single real -wave pairing as the leading instability. Such a -wave pairing can be the or one, which is determined by the axis of the applied strain. Because the resulting or pairing has gap nodes on the FS, the above-mentioned exponential temperature dependence of the experimental quantities in BLS will be replaced by power-law ones. Meanwhile, the U-shaped STM spectrum of BLS will be replaced by a V-shaped one.

In conclusion, we have systematically studied the possible electronic instabilities of doped BLS. The results of our RPA calculations predict that the system is an itinerant ferromagnet in the narrow doping regions around the VH singularities, and a triplet superconductor with a possible high in the vicinity of these regions. With an extra weak Kane-Mele SOC, we further single out the equal-spin helical pairing state as the leading one. This intriguing triplet superconducting state has TR invariant weak topological property, and can harbor the Majorana zero-mode at its boundary tsc2 (); maj1 (); maj2 (); maj3 (), which is useful in the topological quantum computation.

## Acknowledgements

This work is supported by the MOST Project of China (Grants Nos. 2014CB920903, 2011CBA00100), the NSFC of China (Grant Nos. 11274041, 11174337, 11225418, 11334012), and the Specialized Research Fund for the Doctoral Program of Higher Education of China (Grants Nos. 20121101110046, 20121101120046). F.Y. is also supported by the NCET program under Grant No. NCET-12-0038.

## Vi Appendix: RPA with the Kane-Mele SOC

The Kane-Mele SOC term breaks the spin-rotation symmetry, but keeps the spin-rotation symmetry around the -axis, leaving the -component of the total spin to be a good quantum number. Therefore, we define the following susceptibility tensors,

(16) | ||||

(17) | ||||

(18) | ||||

(19) | ||||

(20) | ||||

(21) |

For , we have and

(22) | ||||

(23) | ||||

(24) | ||||

(25) |

With the above expressions of , we consider the scattering of a Cooper pair from the state in the -th () band to the state in the -th () band. This scattering process leads to the following effective interaction vertices:

(30) | ||||

(31) |

Here

(32) | ||||

(33) |

The inversion symmetry, together with the spin-rotation symmetry of our system, enable us to symmetrize the effective interaction vertices into the following channels,

(34) | ||||

(35) | ||||

(36) |

where the index is for the even parity pairing, and for the odd one. From these symmetrized effective interaction vertices, we obtained the following linearized gap equation near the superconducting critical temperature :

(37) |

which replaces Eq. (13) to determine the and the leading pairing symmetry of the system in the presence of the Kane-Mele SOC.

## References

- (1) A. P. Mackenzie and Y. Maeno, Rev. Mod. Phys. 75, 657 (2003).
- (2) N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).
- (3) A. Y. Kitaev, Phys.-Usp. 44, 131 (2001).
- (4) E. Stoner, Philos. Mag. 15, 1018 (1933).
- (5) Y. Nagaoka, Phys. Rev. 147, 392 (1966).
- (6) A. Mielke, J. Phys. A: Math. Gen. 24, L73 (1991).
- (7) A. Mielke, J. Phys. A: Math. Gen. 24, 3311 (1991).
- (8) H. Tasaki, Phys. Rev. Lett. 69, 1608 (1992).
- (9) A. Tanaka and H. Tasaki, Phys. Rev. Lett. 98, 116402 (2007).
- (10) L. Liu, H. Yao, E. Berg, S. R. White, and S. A. Kivelson, Phys. Rev. Lett. 108, 126406 (2012).
- (11) H. Katsura and A. Tanaka, Phys. Rev. A 87, 013617 (2013).
- (12) Y. Li, E. H. Lieb, and C. Wu, Phys. Rev. Lett. 112, 217201 (2014).
- (13) S.-S. Zhang, J. Ye, and W.-M. Liu, arXiv:1403.7031.
- (14) H. Yao and F. Yang, arXiv:1312.0077.
- (15) B. Lalmi, H. Oughaddou, H. Enriquez, A. Karae, S. Vizzini, B. Ealet, and B. Aufray, Appl. Phys. Lett. 97, 223109 (2010).
- (16) P. Vogt, P. DePadova, C. Quaresima, J. Avila, E. Frantzeskakis, M. C. Asensio, A. Resta, B. Ealet, and G. LeLay, Phys. Rev. Lett. 108, 155501 (2012).
- (17) A. Fleurence, R. Friedlein, T. Ozaki, H. Kawai, Y. Wang, and Y. Yamada-Takamura, Phys. Rev. Lett. 108, 245501 (2012).
- (18) L. Chen, C.-C. Liu, B. Feng, X. He, P. Cheng, Z. Ding, S. Meng, Y. Yao, and K. Wu, Phys. Rev. Lett. 109, 056804 (2012).
- (19) L. Meng, Y. Wang, L. Zhang, S. Du, R. Wu, L. Li, Y. Zhang, G. Li, H. Zhou, W. A. Hofer, and H.-J. Gao, Nano Lett. 13, 685 (2013).
- (20) L. Chen, B. Feng, and K. Wu, Appl. Phys. Lett. 102, 081602 (2013).
- (21) B. Feng, Z. Ding, S. Meng, Y. Yao, X. He, P. Cheng, L. Chen, and K. Wu, Nano Lett. 12, 3507 (2012).
- (22) F. Liu, C.-C. Liu, K. Wu, F. Yang, and Y. Yao, Phys. Rev. Lett. 111, 066804 (2013).
- (23) T. Li, Europhys. Lett. 97, 37001 (2012).
- (24) M. L. Kiesel, C. Platt, W. Hanke, D. A. Abanin, and R. Thomale, Phys. Rev. B 86, 020507(R) (2012).
- (25) W.-S. Wang, Y.-Y. Xiang, Q.-H. Wang, F. Wang, F. Yang, and D.-H. Lee, Phys. Rev. B 85, 035414 (2012).
- (26) R. Nandkishore, G.-W. Chern, and A. V. Chubukov, Phys. Rev. Lett. 108, 227204 (2012).
- (27) R. Nandkishore, L. S. Levitov, and A. V. Chubukov, Nat. Phys. 8, 158 (2012).
- (28) L.-D. Zhang, F. Yang, and Y. Yao, Sci. Rep. 5, 8203 (2015).
- (29) C.-C. Liu, H. Jiang, and Y. Yao, Phys. Rev. B 84, 195430 (2011).
- (30) T. Takimoto, T. Hotta, and K. Ueda, Phys. Rev. B 69, 104504 (2004).
- (31) K. Kubo, Phys. Rev. B 75, 224509 (2007).
- (32) K. Kuroki, S. Onari, R. Arita, H. Usui, Y. Tanaka, H. Kontani, and H. Aoki, Phys. Rev. Lett. 101, 087004 (2008).
- (33) S. Graser, T. A. Maier, P. J. Hirschfeld, and D. J. Scalapino, New J. Phys. 11, 025016 (2009).
- (34) T. A. Maier, S. Graser, P. J. Hirschfeld, D. J. Scalapino, Phys. Rev. B 83, 100515(R) (2011).
- (35) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
- (36) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005).
- (37) L. C. Hebel and C. P. Slichter, Phys. Rev. 113, 1504 (1959).
- (38) W. D. Knight, G. M. Androes, and R. H. Hammond, Phys. Rev. 104, 852 (1956).
- (39) D. J. Van Harlingen, Rev. Mod. Phys. 67, 515 (1995).
- (40) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B 78, 195125 (2008).
- (41) G. E. Volovik, The Universe in a Helium Droplet (Oxford Science Publications, New York, 2003).
- (42) X.-L. Qi, T. L. Hughes, S. Raghu, and S.-C. Zhang, Phys. Rev. Lett. 102, 187001 (2009).