Triplet p-wave pairing correlation in low-doped zigzag graphene nanoribbons

We reveal an edge spin triplet p-wave superconducting pairing correlation in slightly doped zigzag graphene nanoribbons. By employing a method that combines random-phase approximation, the finite-temperature determinant quantum Monte Carlo approach, and the ground-state constrained-path quantum Monte Carlo method, it is shown that such a spin-triplet pairing is mediated by the ferromagnetic fluctuations caused by the flat band at the edge. The spin susceptibility and effective pairing interactions at the edge strongly increase as the on-site Coulomb interaction increases, indicating the importance of electron-electron correlations. It is also found that the doping-dependent ground-state p-wave pairing correlation bears some similarity to the famous superconducting dome in the phase diagram of a high-temperature superconductor, while the spin correlation at the edge is weakened as the system is doped away from half filling.

Introduction: Triplet superconductivity (SC) has been the focus of modern condensed matter physics due to its possible connection with topological quantum information and computation [1][2][3][4][5][6][7][8][9].It is proposed that gapless majorana bound state would localize at the end of one dimensional spinless p−wave superconductor [1], which could be used to practical realization of topological quantum computation [10,11].Such Majorana bound state has been proposed to be realized in real material by superconducting proximity effect [12][13][14] and evidence of its existence was reported by experiments recently [15].Here we explore the possibility of intrinsic triplet SC, which is generated by electronic correlation.
In this paper, we reveal a possible edge spin triplet p-wave superconducting pairing correlation in slightly doped zigzag graphene nanoribbons with appropriate interactions.Graphene, a single layer of carbon, has generated immense interest ever since its experimental discovery [16,17].Recently, experimental advances in doping methods have allowed to change the type of carriers, electrons or holes [18,19], opening the doors for exotic phases like SC and magnetism induced by repulsive interactions.For instance, it was shown by the two-stage renormalization group calculation that unconventional SC is induced by weak repulsive interactions in honeycomb Hubbard models away from half-filling [20], and a topological d+id SC in a heavily doped system [21][22][23].At graphene edges the density of states may be peaked due to the presence of edge-localized states close to the Fermi level [24].Especially at extended zigzag edges this leads to a phenomenon called edge magnetism, where various theories [25,26] predict ferromagnetic (FM) intraedge and antiferromagnetic (AFM) interedge correlations.From these discoveries, a question which naturally arises is, whether there is triplet SC mediated by the FM spin correlations on each edge in the doped zigzag graphene nanoribbons?
In the present work, we establish the p-wave supercon- ducting pairing correlation at the edge of zigzag graphene nanoribbons by using combined random phase approximation (RPA) [27][28][29][30][31], the finite temperature determinant quantum Monte Carlo (DQMC) [32,33] and the ground state constrained path quantum Monte Carlo (CPQMC) [34,35] methods.Our unbiased results show that, both the ferromagnetic spin correlation and the effective p-wave superconducting pairing correlation are The Fermi level of the half-filled system is marked by the red dash lines in both figures.
enhanced greatly as the interaction increases.Model: The ribbon geometry considered here is depicted in Fig. 1, in which the blue and white circles represent sublattice A and B respectively, and the transverse integer index 1, 2, ..., L y defines the width of the ribbon while 1, 2, ..., L x at the zigzag edge defines the length.Assuming the ribbon to be infinite in the x direction but finite in the y direction, we produce a graphene nanoribbon with zigzag edges, and its electronic and magnetic properties can be well described by the following Hubbard model [17] where c † iσ is the electron creation operator at site i and with spin polarization σ =↑, ↓, U labels the on-site repulsive interaction, and µ is the chemical potential.Here the t and t terms describe the nearest neighbor (NN) and next nearest neighbor (NNN) hoppings, respectively.In the following study, we adopted t = −0.1t,which is consistent with experiments [36].In our calculation, we employ periodic boundary conditions at x direction and open boundary conditions at the zigzag edge.In Fig. 2, the carrier distribution (a) as a function of site index at U = 2.0t and (b) from edge → bulk → edge with different interactions is shown.It is clear to see that most charge carrier distributes along the edge, and the increasing interaction pushes much more charge carrier to edges.
The band structure of a six-chain nanoribbon system is shown in Fig. 3(a).Here, as the system is periodic in the x-direction, the momentum k x is a good quantum number.From Fig. 3(a), one finds a flat band bottom with energies locating near the Fermi level (≈ −0.2t) of the half-filled system.Physically, such a flat band bottom is caused by the edge states, which leads to the DOS peak at around −0.2t shown in Fig. 3(b).
RPA study: Guided by the idea that triplet SC may be mediated by the strong FM spin fluctuations in the  system, we performed a RPA-based study on the possible pairing symmetries of the system.Multi-orbital RPA approach [27][28][29][30][31], which is a standard and effective approach for the case of weak coupling limit, is engaged in our study for small U (< 0.01t).Various bare susceptibilities of this system are defined as where l i (i = 1, 2L y ) denote orbital (sublattice) indices.The largest eigenvalue χ(q x ) of the susceptibility matrix χ (0) l,m (q x ) ≡ χ (0)l,l m,m (q x , iν = 0) as function of q x is shown in Fig. 4(a) for three different dopings near halffilling.As a result, the susceptibility for the doping δ = 3% with chemical potential µ = −0.2tpeaks at zero momentum, which suggests strong FM intra-sublattice spin fluctuations in the system.Further more, from eigenvector of the susceptibility matrix, one can obtain the pattern of the dominating spin fluctuation in the system, which is shown in Fig. 4(b).Obviously, the dominating spin fluctuation, which mainly locates on the two edges, is FM on each edge and AFM between the two edges.When µ deviates from −0.2t, the susceptibility peaks deviate from zero, as shown in Fig. 4(a), suggesting weaker FM spin fluctuations in the system.
With weak Hubbard-U , the spin (χ s ) or charge (χ c ) susceptibilities in the RPA level are given by Clearly, the repulsive U suppresses χ c but enhances χ s .Thus, the spin fluctuations take the main role of mediating the cooper pairing in the system.In the RPA level, via exchanging the spin fluctuations represented by the  spin susceptibilities, the cooper pairs near the FS acquire an effective interaction V eff [27,31], from which one solves the linearized gap equation to obtain the leading pairing symmetry and its critical temperature T c .Focusing on the low-doping regime for which the chemical potential µ locates within the flat band bottom, we obtained the largest pairing eigenvalues λ for the singlet and triplet pairings as function of µ for a 6-chain ribbon with weak interaction U = 0.001t, as shown in Fig. 5(a).Interestingly, while both pairings possess their largest eigenvalues at µ = −0.2t(3% doping) due to the DOS peak there (as shown in Fig. 3(b)), the triplet pairing wins over the singlet one in the low doping regime at the flat band bottom.Physically, the triplet pairing in this regime is mediated by FM spin fluctuations shown in Fig. 4. In Fig. 5(b), the results for U = 0.005t are shown.Comparing (a) and (b), it's obvious that stronger interaction leads to qualitatively the same but quantitatively stronger pairing correlations than weak interaction.In Fig. 5(c) and (d), the results for a 4-chain ribbon and 8-chain ribbon are shown with U = 0.001t.The results for all these cases are qualitatively similar.
Note that we have chosen very weak U in our RPA calculations, as for U > U c ≈ 0.007t (for µ = −0.2t), the divergence of the spin susceptibility invalidates our calculations.Physically, such a spin susceptibility divergence will not lead to a magnetically-ordered state since the Mermin and Wagner's theorem prevents a one dimensional system to form long range order.Instead, short-ranged FM spin correlations here might mediate triplet superconducting pairing correlations.We leave the study for the case of U > U c to the following DQMC and CPQMC approach, which is suitable for strong coupling problems.
QMC Result: As FM fluctuations play an essential role, we first study the magnetic correlations.Specifi-cally, we compute the spin correlation S i,j = S i • S j in the z direction, and define the uniform spin susceptibility at zero frequency, To investigate the SC property, we compute the effective pairing interaction and study the distance dependent pairing correlation.The effective pairing interaction is extracted from the pairing susceptibility, with Here α stands for the pairing symmetry, f α (δ l ) is the form factor of pairing function, and the vectors δ l (l = 1, 2) denote the NNN site along the edge.In order to extract the effective pairing interaction, the bubble contribution p α (i, j) should be evaluated by replacing c † i↓ c j↓ c † i+δ l ↑ c j+δ l ↑ in Eq. ( 5) with c † i↓ c j↓ c † i+δ l ↑ c j+δ l ↑ .Thus we have the effective pairing interaction P α = p α − p α .The corresponding pairing correlation is defined as Considering the special structure of graphene zigzag nanoribbons shown in Fig. 1, the interesting pairing correlation in such system is the pairing between sites on the same sublattice, and two form factors shall be studied p-wave: f p (δ l ) = e i(l−1)π , l = 1, 2 In Fig. 6 (a), the edge spin susceptibility χ is shown as a function of temperature with different U at δ=0.02.The edge χ defined here is to make a summation over the sites on the edge, as that marked as larger circles in Fig. 1.It is interesting to see that χ increases as the temperature decreases, which indicates a dominant FM fluctuations on the zigzag edge.The χ increases as U increases, indicating that the electronic correlation is important on the magnetic excitation in such system.The uniform spin susceptibility χ B for the whole system is also shown, which decreases slightly as the temperature decreases.To further reveal the FM correlation on the zigzag edge, the spin-spin correlation along the edge is shown in Fig. 6 (b).One can see that, the spin correlation S 1i (i = 2, 3, • • • ) along the edge is always positive, suggesting FM correlation.One may also see that the spin correlation is weakened as the system is doped away from half fillings, agreeing with the results indicated by RPA shown in Fig. 4(a).In Fig. 7, we plot the effective pairing interaction P α as a function of temperature for different U and electron fillings on a lattice with 2 × 4 × 12 sites.Clearly in Fig. 7, the intrinsic pairing interaction P α is positive and increases with the lowering of temperature.Such a temperature dependence of P α suggests effective attractions generated between electrons and the instability toward SC in the system at low temperatures.Moreover, Fig. 7 shows that the intrinsic pairing interaction for pwave symmetry enhances with larger U , indicating the enhanced pairing strength with the enhancement of the electron correlations.As for the another extensive-s pairing symmetry, our DQMC results yield negative intrinsic pairing interactions (not shown here).
Numerical approaches like DQMC, however, had its own difficulties, typically being limited to small sizes, high temperature, and experience the infamous fermion sign problem, which cases exponential growth in the variance of the computed results and hence an exponential growth in computer time as the lattice size is increased and the temperature is lowered [32].In general, to determine the superconducting pairing symmetry by numerical calculation for finite size models, we have to look at the distance dependent pair-correlation function at zero temperature.To shed light on this critical issue, it is important to discuss the results obtained from CPQMC method [34,35]  many other ground state observables for large system [34].In Fig. 8, we compare the pairing correlations on lattices with 2 × 6 × 24 sites for different electron fillings at U = 2.0t.Here, the simulations are performed for the closed-shell cases.The distance dependent pairing correlations for δ = 3/144 0.021 (dark triangle), δ = 5/144 0.035 (red circles ), and δ = 7/144 0.049 (blue square) are shown.One can readily see that C p (r) decreases as the distance increases, and the decay of the distance dependent pairing correlations is different for different dopings.In the inset of Fig. 8, the pairing correlation C p (r = 12) at the largest distance is shown as a function of doping.In the filling range we investigates, C p (r = 12) is not monotonic function of the doping and there exists an optimal doping (around 0.035 electron/site) where the magnitudes of C p (r = 12) are maximized.This result is consistent with that of RPA, where the doping dependent pairing correlation bears some similarity to the famous superconducting "dome" in the phase diagram of high-temperature superconductors [37], while the spin correlation at the edge is weakened as the system is doped away from half filling.
Discussions: We performed combined RPA and quantum Monte Carlo study of the edge magnetic and pairing correlation in low doped zigzag graphene nanorib-bons.Our studies show that the triplet edge p-wave SC occurs as the ground state of our model system.The optimal doping is around 0.03, which can be easily understood as the DOS peaks there, and it is in the current experimental capability for graphene-based material.Our accurate numerical results establish the properties of the p-wave superconducting correlation in zigzag graphene nanoribbons, and will be important for any experimental scheme aimed at detecting the p-type superconducting state, as such scheme will likely be based on the distinctive properties of the edge.

FIG. 1 .FIG. 2 .
FIG. 1. (Color online) A piece of a honeycomb lattice displaying zigzag edges with Ly = 4 which defines the width of the ribbon in transverse direction and Lx=12 indicating the length in longitudinal direction.The lattice sites at zigzag edge are much larger than the sites in the bulk indicates the charge carrier are moving along the edge.

FIG. 3 .
FIG. 3. (Color online) Band structure (a) and DOS (b) of a six-chain nanoribbon system.Note that the flat band bottom locating at about −0.2t in (a) leads to the DOS peak in (b).The Fermi level of the half-filled system is marked by the red dash lines in both figures.

FIG. 6 .
FIG. 6. (Color online) (a)The edge χ as a function of temperature at δ=0.02 for different U , and the uniform χB for U = 2.0t is also present.(b) The spin correlation as a function of site index along the edge for U = 2.0t at δ=0.02 and δ=0.04, as well as U = 1.0t at δ=0.02.
on larger lattice.In a variety of benchmarking calculations, the CPQMC method has yielded very accurate results of the ground state energy and FIG. 8. (Color online) The p-wave superconducting pairing correlation as a function of the distance r on a lattice with 2× 6 × 24 sites.Inset: the doping dependent pairing correlation at r = 12.