Possible Electric-Field-Induced Superconducting States in Doped Silicene

Silicene has been synthesized recently, with experimental evidence showing possible superconductivity in the doped case. The noncoplanar low-buckled structure of this material inspires us to study the pairing symmetry of the doped system under a perpendicular external electric field. Our study reveals that the electric field induces an interesting quantum phase transition from the singlet chiral d + id′-wave superconducting phase to the triplet f-wave one. The emergence of the f-wave pairing results from the sublattice-symmetry-breaking caused by the electric field and the ferromagnetic-like intra-sublattice spin correlations at low dopings. Due to the enhanced density of states, the superconducting critical temperature of the system is enhanced by the electric field remarkably. Furthermore, we design a particular dc SQUID experiment to detect the quantum phase transition predicted here. Our results, if confirmed, will inject a new vitality to the familiar Si-based industry through adopting doped silicene as a tunable platform to study different types of exotic unconventional superconductivities.

T he past few decades have witnessed the glorious development of unconventional superconductivity (SC) 1 . While most of the so far confirmed unconventional superconductors are singlet pairing including, say, the cuprates with d-wave pairing and the iron-pnictide with s 6 pairing 2,3 , no triplet superconductor has been confirmed except for the Sr 2 RuO 4 system which shows evidence for possible p 1 ip9 pairing 4 . Recently, however, triplet pairing is catching more and more attentions partly due to its possible connection with topological SC [5][6][7][8] , which has inspired great enthusiasm in searching for triplet superconductors. Here we predict a tunable quantum phase transition (QPT) from the singlet d 1 id9 superconducting state to the triplet f-wave one in doped silicene via applying a perpendicular external electric field on the system. Although both pairing symmetries break the time reversal symmetry, they are quite different. While the former shows chiral complex gap structures which is a hot topic recently [9][10][11][12][13][14][15] , the latter has a real pairing gap which changes sign with every 60u rotation. The f-wave pairing has been proposed in the study of the quasi-1D organic system 16 , the triangular lattice system for Na x CoO 2 17 , and the spinless fermion system on the honeycomb optical lattice 18 . However, unambiguously confirmed experimental evidence for it is still lack now. If the f-wave SC predicted here is confirmed experimentally, it will be the first time to realize such an intriguing high-angular-momentum pairing state.
Silicene, a single atomic layer of Si forming a 2D honeycomb lattice, can be regarded as the Si-based counterpart of graphene. It has attracted a lot of research interest since being successfully synthesized recently [19][20][21][22][23] . Similar lattice structures between graphene and silicene bring them similar band structures with linear dispersion near the Fermi surface (FS), which further lead to similar physical properties between them. The most important structural difference between the two layered systems lies in that, while the graphene layer forms a regular flat plane, the silicene layer instead takes the form of noncoplanar low-buckled (LB) structure, with the sublattices A and B forming two separate planes. The LB structure of silicene causes many fascinating consequences, such as the enhanced quantum spin Hall effect 24,25 , the quantum anomalous Hall effect and valley polarized quantum Hall effect in external electric field [26][27][28][29] , and possible d 1 id9 chiral SC in bilayer silicene 30 . On the other hand, this LB structure provides the possibility to use gated silicene as a tunable source of the perfectly spin-polarized electrons 31 . Interestingly, a recent tunnelling experiment reported electronic gap in doped silicene 32 , probably caused by SC, which has attracted a lot of research interests 33,34 .
In this paper, we report our study of the pairing symmetry of doped silicene under a perpendicular external electric field, which differentiates the on-site energies of the two sublattices of this LB system. Based on the random phase approximation (RPA) to the Hubbard model of the system, we reveal a tunable QPT from the singlet d 1 id9 chiral SC under weak field to the triplet f-wave one under strong field, without the necessity of long-range Coulomb interaction [35][36][37] . The physics behind this interesting QPT lies in that, under the strong sublattice-symmetry-breaking electric field, only one sublattice is left as the low energy subspace, and the ferromagnetic-like intra-sublattice spin correlations favor triplet f-wave pairing in this subspace. Due to the enhanced density of states (DOS) near the Fermi level, the superconducting critical temperature of the system is enhanced by the applied electric field remarkably. To detect the QPT predicted here, we further design a particular experiment with a phase-sensitive dc SQUID, which can distinguish between the d 1 id9 and f-wave pairings. Our study will open up a new era to utilize the familiar Si-based material as a tunable platform to study the competition and QPT among different types of exotic unconventional SCs.

Results
Model. The LB honeycomb lattice of silicene is shown in Figs. 1(a) and 1(b). Since the two sublattices A and B lie in two parallel planes separated from each other, an applied perpendicular electric field induces an on-site energy difference D between them. We adopt the following Hubbard model as an appropriate start point to study the low energy physics of the system: Here the t-term (with t 5 1.12 eV 25 ) describes the nearest-neighbor (NN) hoppings, and {1 ð Þ pi~z 1 {1 ð Þ for sublattice A (B). The Hubbard interaction strength U 5 2 eV is adopted in the following, which is near the value obtained from the first principle calculations 38 . The next-nearest-neighbor (NNN) hoppings are not included here because they are qualitatively unimportant. The resulting particle-hole symmetry enables us to focus the following discussions only on the electron-doped case.
The band structure of the system with the dispersion e + k~+ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi D 2 2 zt 2 1ze ikx ze {iky 2 s is shown in Fig. 1(c). Clearly, the applied field induces a band gap D near the K-point and the nearby low energy band structure is flattened, although it does not modify the shape of the FS for a given doping concentration. Note that the electric field breaks the sublattice symmetry of the system. In particular, the on-site energy of sublattice A, i.e., V A 5 D/2, is closer to the Fermi level shown for the electron-doped case, and consequently this sublattice will dominate the low energy physics of the system. This effect turns out to be very important for our following discussions.
Susceptibilities. According to the standard RPA approach 30,39-41 , we define the free susceptibility (U 5 0) of the model as where p, q, s, t 5 1, 2 is the sublattice index. The Hermitian static susceptibility matrix is defined as The largest eigenvalue of this matrix represents the static susceptibility of the system in the strongest channel, and the corresponding eigenvector determines the pattern of the applied detecting field in that channel. In addition, the eigenvector also describes the pattern of the dominant intrinsic spin fluctuations in the system.
In Fig. 2(a)-2(c), we show the k-space distributions of the zero temperature static susceptibility of the system in the strongest channel mentioned above. From Fig. 2(a) to 2(c), one can clearly observe the doping evolution of the static susceptibility. In particular, when the doping increases gradually from zero to the Van Hove (VH) doping x 5 1/4, the momenta of the maximum susceptibility evolves from the C-point (2(a)) first to a small circle around it (2(b)), and finally to the M-points (2(c)). Such an doping evolution of the susceptibility originates from the evolution of the FS which grows gradually from the K-points (the Dirac-points) first to small pockets around them, and finally to the connected large FS with perfect nesting at the VH singularity. From 2(a) to 2(c), one verifies that the dominant spin fluctuations on each sublattice of the system changes gradually with doping from ferromagnetic-like to antiferromagnetic-like. In Fig. 2(d), a typical pattern of the dominant spin fluctuations of the system is shown at 10% doping for D 5 1 eV. Most prominently in Fig. 2(d), the magnetic moments are asymmetrically distributed between the two sublattices, with the magnitude of the moment on sublattice A obviously larger than that on sublattice B. This is consistent with the above band structure analysis which suggests that sublattice A will dominate the low energy physics of the system. Our calculation reveals that such asymmetry is enhanced by both electric-field and doping.
When the Hubbard interaction is turned on, the charge (c) or spin (s) susceptibility is given by RPA as where (U) is a 4 3 4 matrix, whose only two nonzero elements are U ð Þ mm mm~U m~1,2 ð Þ 30 . Clearly, the repulsive Hubbard interaction suppresses x (c) and enhances x (s) . When the interaction strength U   is greater than a critical value U c , the spin susceptibility x (s) of the model diverges, which implies the instability towards the long-range spin-density-wave (SDW) order. Below the critical interaction strength U c , the spin fluctuations take the main role in mediating the Cooper pairing.
Pairing symmetries. Through exchanging the renormalized susceptibilities (3), one obtains the effective interaction vertex V ab (k, k9) 41 between two Cooper pairs on the FS, which is then plugged into the following linearized gap equation near the superconducting critical temperature T c 41 : Here the integration is along various FS patches labelled by a or b, v b F k' ð Þ is the Fermi velocity and k' jj is the tangential component of k9 along the FS. Solving Eq. (4) as an eigenvalue problem, one obtains the largest eigenvalue l, which determines T c via T c < te 21/l , and the corresponding eigenvector D a (k) as the leading gap form factor.
The phase diagram on the x-D plane obtained by our RPA calculations is shown in Fig. 3(a). Clearly, except for the SDW phase near the VH doping x 5 1/4, the singlet chiral d x 2 {y 2 zid xy and triplet nodeless f-wave pairings beat other instabilities and serve as the leading instability in different regimes of the phase diagram. The gap function of the d xy symmetry (for doping x 5 0.15 and D 5 0) shown on the FS in Fig. 3(b) is antisymmetric about the x and y axes, and that of its degenerate partner d x 2 {y 2 (not shown) is symmetric about these axes. Their mixing in the form of d 1 id9 minimizes the ground state energy. The main part of the d 1 id9 pairing in real space is distributed on NN-bonds as shown in Fig. 3(d). The gap function of the fwave pairing (for doping x 5 0.15 and D 5 1 eV) is shown in Fig. 3(c). The main part of the f-wave pairing in real space is distributed on NNN-bonds as shown in Fig. 3(e). Clearly, the gap function of this time-reversal-breaking f-wave pairing changes sign with every 60u rotation either in k-space or in real space. As shown in Fig. 3(a), the leading instability for D 5 0 is d 1 id9 pairing at low dopings and SDW order near the VH doping. This result is consistent with previous calculations on the graphene system 9-14 , which shares the same Hubbard model as here. In Fig. 3(f), we show the doping dependence of the SDW critical interaction U c . Clearly, when the doping x * > 0:2, U c is less than U 5 2 eV, which accounts for the emerging of the SDW state in Fig. 3(a). The most interesting and important discovery here is that the triplet f-wave pairing, which is mediated by the NNN-bond ferromagnetic spin fluctuations, wins over the d 1 id9 pairing and rises as the leading pairing symmetry when sufficiently strong electric-field is applied on the system at low dopings. Physically, such an electric-field-induced QPT originates from the sublattice-symmetry-breaking effect mentioned before, which has selected the sublattice A as the low energy subspace of the system. All the relevant low energy physics, including the pattern of the spin fluctuations shown in Fig. 2(d), the effective pairing interaction mediated by these spin fluctuations, and the pairing itself take place mainly in this subspace in strong electric field. Consequently, at low dopings, the inter-sublattice pairing d 1 id9 symmetry shown in Fig. 3(d) has to give way to the intra-sublattice pairing f-wave symmetry shown in Fig. 3(e).
The superconducting critical temperature T c < te 21/l of the f-wave pairing state is controllable, and can be remarkably enhanced by the applied electric field. As shown in Fig. 3(g) for doping x 5 0.15, the eigenvalue l of the f-wave pairing can be enhanced to about 0.18 when D is tuned to 1 eV. Although the RPA approach tends to overestimate T c , we still have space to enhance the T c of the system to the experimentally accessible range by tuning the doping level closer to the VH doping or increasing the electric field. Physically, such enhancement is attributed to the increase of the DOS near the Fermi level (see the inset of Fig. 3(g)) caused by the flattening of the band structure under the electric field (see Fig. 1(c)). The increase of the DOS not only enhances the number of Cooper pairs near the FS, but also enhances the pairing interaction. The T c of the d 1 id9 pairing can also be enhanced by sufficiently strong electric field due to this DOS enhancement, but it is lower than that of the f-wave pairing as shown in Fig. 3(g). For a vanishingly small U, the f-wave SC has been proposed even in the electric-field-free system 42 , but with extremely low T c .
Detecting the QPT. The QPT from d 1 id9 to f-wave pairings predicted here can be detected by the dc SQUID, a phase-sensitive device which has been adopted in determining the pairing symmetries of such superconducting systems as cuprates 43 and Sr 2 RuO 4 44 . Our basic scheme is shown in Figs. 4(a) and 4(d), where a slice of silicene is fabricated into a hexagonal shape, allowing the relative phase among different directions in the system to be detected. Two superconductor-normal metalsuperconductor (SNS) Josephson tunneling junctions are formed on the opposite (4(a)) or adjacent (4(d)) edges of the hexagon, which are connected by a loop of a conventional s-wave superconductor, forming a bimetallic ring with a magnetic flux W threading through the loop.
As a result of the interference between the two branches of Josephson supercurrent, the maximum total supercurrent (the critical current) I c in the circuit modulates with W according to Here I 0 is the critical current of one Josephson junction, W 0 5 h/2e is the basic flux quantum, and d ab is the intrinsic phase shift inside the silicene system between pairs tunneling into the system in two direc-   symmetry is characterized by the minima of I c (which can be nonzero when the self-inductance of the loop is not negligible) at zero flux for both configurations 4(a) and 4(d). For the d 1 id9 symmetry at zero flux, while the critical current I c is a maximum in configuration 4(a), it is not an extreme point in 4(d). Thus, by observing the modulation of the SQUID response vs applied magnetic flux, one can distinguish between the two pairing symmetries and hence detect the QPT.

Discussion
The non-monotonic doping dependence of the critical electric field in Fig. 3(a) originates from the competition between two opposite effects. On the one hand, the sublattice-symmetry-breaking effect, which favors the f-wave pairing, is enhanced by doping since the Fermi level sits in the middle of V A and V B at zero doping (see Fig. 1(c)). On the other hand, as shown in Fig. 2(a)-2(c), the intrasublattice spin correlations in the system evolves from ferromagnetic-like at low dopings, which favors triplet f-wave pairing, to antiferromagnetic-like near the VH doping, which favors singlet d 1 id9 pairing.
Besides the nodeless f-wave symmetry predicted here, there is also another nodal f9-wave symmetry which is mainly composed of intrasublattice pairings with a bond length of three times of the lattice constant. This f9-wave symmetry is not favored here, for its gap nodes along the C-K lines do not avoid the FS. When a sufficiently strong Kane-Mele spin-orbit coupling (SOC) 45,46 is added into model (1), the Fermi pockets of the system would shift from near the Kpoints to near the M-points, and the f9-wave symmetry would be favored, for its gap nodes avoid the FS. In a system with a medium SOC, the two f-wave pairings could be nearly degenerate, and their mixing in the form of f 1 if9 might be realized to minimize the energy. As a result of its nontrivial topological property, this intriguing triplet chiral superconducting state can harbor the Majorana zero-mode at its boundary 6,[47][48][49] , which is useful in the topological quantum computation. Although the SOC in the present silicene system is too small 24,25 to favor the f 1 if9 pairing, we can expect the system with a LB honeycomb lattice similar to silicene and a medium SOC, such as BiH 50,51 , could serve as the platform for this novel high-angular-momentum chiral pairing state. We leave this subject for future studies.
As a result of the rapid development of modern experimental techniques, the electric field strong up to 0.3 V/Å has already been applied to the research works of materials 52 , which has approached what we need here. Besides, the monolayer of silicene has been synthesized on top of the substrates of Ag, ZrB 2 , Ir, etc [19][20][21][22][23] . Due to the low-buckled structure of silicene, the distance between the sublattice A and the substrate is different from that between the sublattice B and the substrate. As a result, the interaction between the silicene monolayer and the substrate provides an effective electric field, which causes the sublattice-symmetry-breaking and the difference between the on-site energies of the two sublattices 53 . Since such an effective electric field is an internal electric field of chemical origin, it can generally be as strong as what we need here. Therefore, our proposal of the f-wave pairing is feasible in practice.
In conclusion, we have systematically studied the pairing phases of doped silicene in a perpendicular external electric field. The results of our RPA study predict that with the enhancement of the electric field, the system will experience a QPT from singlet d 1 id9 superconducting state to triplet f-wave one, and the superconducting critical temperature of the system will be enhanced due to the increase of the DOS. Our model needs neither long-range Coulomb interaction nor the situation of vanishingly small Hubbard-U, and thus is more realizable than other proposed ones.

Methods
Susceptibilities. According to the standard RPA approach 30,39-41 adopted in our study, we first define the free susceptibility (U 5 0) of the model (1) as in Eq. (2). Direct calculation yields the following explicit expression of the free susceptibility, where e a k and j a k are the a-th eigenvalue and eigenvector of the single particle Hamiltonian of the system respectively, and n F is the Fermi-Dirac distribution function.
When the interaction is turned on, we define the charge (c) and spin (s) susceptibilities as For U 5 0, we have x (c) 5 x (s) 5 x (0) . In the RPA level, x (c(s)) is given by Eq. (3).
Effective interaction and gap equation. Consider the scattering of a Cooper pair from the state (k9, 2k9) in the b-th (b 5 1, 2) band to the state (k, 2k) in the a-th (a 5 1, 2) band. This scattering process can be described by the following effective interaction, Here the projective interaction vertex V ab (k, k9) is given by the effective vertex C pq st k,k' ð Þthrough The mean-field decoupling of the effective interaction (9) in the Cooper channel gives rise to the self-consistent gap equation. Near the critical temperature T c , this equation is linearized as Eq. (4), solving which one obtains the leading pairing symmetries and the corresponding T c . If the leading pairing symmetries include two degenerate doublets such as d x 2 {y 2 and d xy , we can express the gap function of the system as Here d a x 2 {y 2 k ð Þ and d a xy k ð Þ represent the normalized gap functions of corresponding symmetries. Then the mixing coefficients C 1 , C 2 , and C 3 are determined by the minimization of the total mean-field energy with the constraint of the average electron number in the superconducting state.