# Triplet $p+ip$ paring correlation in doped Kane-Mele-Hubbard model: A quantum Monte Carlo study

## Abstract

By using the constrained-phase quantum Monte Carlo method, we performed a systematic study of the pairing correlations in the ground state of the doped Kane-Mele-Hubbard model on a honeycomb lattice. We find that pairing correlations with symmetry dominate close to half filling, but pairing correlations with symmetry dominate as hole doping moves the system below three-quarters filling. We correlate these behaviors of the pairing correlations with the topology of the Fermi surfaces of the non-interacting problem. We also find that the effective pairing correlation is enhanced greatly as the interaction increases, and these superconducting correlations are robust against varying the spin-orbit coupling strength. Our numerical results suggest a possible way to realize spin triplet superconductivity in doped honeycomb-like materials or ultracold atoms in optical traps.

###### pacs:

71.10.Fd,74.20.Mn,74.20.Rp## I Introduction

Recently, a new research field has emerged in condensed matter physics which is based on the findings that a spin-orbit interaction can lead to topological electronic phase transitions Kane2005 () and the electron-electron interactions can produce these transitions Sun2009 (); Liu2010 (); Raghu2008 (); Wen2010 (); Rachel2010 (); Yu2011 (); Herbut2013 (); HHH2013 (). This new field makes it natural to ask what remarkable new properties and transitions might occur between distinct, competing correlated topological electron phases when the spin-orbit interaction is present. Various researchers have proposed Raghu2008 (); Wen2010 (); Rachel2010 (); Yu2011 (); Herbut2013 (); HHH2013 () the Hubbard model with spin-orbit interactions (the Kane-Mele-Hubbard model) on a honeycomb-like lattice, as a starting point for answering these questions.

Indeed the density of states of the non-interacting version of the model on a honeycomb lattice has special features that point to special properties in its interacting version. In particular, if the hopping is nearest neighbor, a van Hove singularity (VHS) exists in the non-interacting model’s density of states (DOS) at electron filling Ma2010 (). This feature prompted several groups Wen2011 (); Nandkishore2012 (); Thomale2012 () to suggest singlet superconductivity in the interacting model, and several others Schaffer2007 (); Jiang2008 (); Baskaran2010 (); Ma2011 (); Zhong2014 () to suggest pairing in the low doped case. Numerical simulations point to other novel states of matter tied to the strength of the electron-electron interaction and spin-orbit coupling. In particular, if is the strength of the local Coulomb interaction and the hopping amplitude to nearest neighbor sites, the Hubbard model on a honeycomb lattice is known to have a semimetal phase at small and an antiferromagnetic one at large Paiva2005 (), and ground-state quantum Monte Carlo (QMC) simulations performed at half-filling suggest a transition from a semimetal to an insulating antiferromagnetic state insulator when 3.5Herbut2013 (); Paiva2005 (). On the other hand, other quantum Monte Carlo studies, such as the finite temperature determinant quantum Monte Carlo method Ma2010 (), find that strong antiferromagnetic fluctuations dominate around half filling and strong ferromagnetic correlations dominate at less than 3/4-filling. Accordingly, the Kane-Mele-Hubbard model on a honeycomb lattice is an opportunity to study a variety of novel quantum critical and thermal fluctuations of itinerant charge carriers, with doping and interacting strengths being very relevant control parameters that may lead to new topologically-constrained states of matter.

In this paper we study quantum fluctuations in the Kane-Mele-Hubbard model via ground state quantum Monte Carlo simulations. Because of the already observed strong magnetic fluctuations, we ask whether superconducting fluctuations exist in novel forms. A common view, supported by recent experiments Tacon2011 (), is that the magnetic fluctuations provide a glue for pairing and generally lead to unconventional superconductivity Scalapino2012 (). Relative to what we learned from past numerical studies on this model Ma2014 (), at least two important questions remain: What happens to the pairing for dopings below the van Hove singularity? and How does the spin-orbit coupling affect the pairing at any doping? In conjunction with these two questions, we can also ask whether spin triplet superconductivity is possible. A topological superconductor with pairing symmetry is the most elemental way in which a non-Abelian topological state can emerge as the ground state of a many-body system Nayak2010 () and thus provides an ideal platform to construct a possible quantum computer Stern2013 ().

We address the questions stated above by using the constrained-phase quantum Monte Carlo method Ortiz1997 (); Schmidt1999 (); Schmidt2001 (); Schmidt2003 (); Sarsa2003 (); Zhang2003 (). This method allows us to perform simulations not possible with more commonly used quantum Monte Carlo methods. As we show, our results reveal rich physics when both the spin-orbit coupling and the electron-electron interactions are strong. Between half-filling and the filling at the VHS, we find that the pairing dominates. Below the lower band VHS filling, we find strong triplet superconducting correlations. Further, the effective pairing correlation is enhanced greatly as the interaction increases, and these superconducting correlations are robust against varying the spin-orbit coupling. On the basis of these results, we suggest that triplet topological superconductivity might be realizable in doped honeycomb-like materials Novoselov2005 (); Vogt2012 (); Cahangirov2009 () or ultracold atoms in optical traps Gurarie2005 (); Wu2007 (); Tarruell2012 ().

## Ii Model and Methods

The Kane-Mele-Hubbard model Kane2005 () on the honeycomb lattice is

(1) |

() where , , and are the nearest neighbor (NN) hopping energy, the strength of the on-site repulsion, and the next-nearest neighbor (NNN) spin-orbit coupling strength, respectively. Here () annihilates (creates) an electron with spin on site and depending on if the electron makes as “right” or “left” turn when going from site to site on the honeycomb lattice. We assume periodic boundary conditions.

The constrained-phase quantum Monte Carlo method (PCPMC) Ortiz1997 (); Schmidt1999 (); Schmidt2001 (); Schmidt2003 (); Sarsa2003 (); Zhang2003 () is a generalization of the constrained-path method (CPMC) Zhang1995 (); Zhang1997 (); Carlson1999 () and is an analog of the fixed-phase generalization of the fixed-node diffusion Monte Carlo method Ortiz1993 (). The PCPMC method has yielded very accurate results for the ground state energy and other ground state observables for various strongly-correlated lattice models Chang2008 (); Chang2010a (); Chang2010b (); Xu2011 () and for atoms, molecules, and nuclei Schmidt1999 (); Schmidt2001 (); Schmidt2003 (); Sarsa2003 (). The PCPMC method is sometimes called the phaseless method Zhang2003 ().

Briefly, for a Hamiltonian , the constrained-path method, as most ground state quantum Monte Carlo methods, projects to the ground state from some initial state via

(2) |

where is the non-interacting part of the Hamiltonian and is the interacting part. Here, is the hopping and spin-orbit parts of (1) and it is the Hubbard -term. In the constrained-path method, a Hubbard-Stratonovich transformation is applied to each exponential of the interaction. With the consequence of introducing distributions of auxiliary scalar fields into the problem, converting the exponential of the interaction into an exponential of a non-interacting term, depends on imaginary-time dependent external fields. With an initial state being an ensemble of Slater determinant, the constrained-path method propagates one ensemble of Slater determinants into another. The individual propagating determinants are called random walkers. At any step in the projection, the projected state is a weighted sum of the walkers, .

After a sufficient number of steps, the method begins to sample walker weights that represent those of the ground state wave function. For most Fermion simulations these weights are generally not all non-negative and hence do not represent a probability distribution amenable to Monte Carlo sampling. This is called the sign problem.

The constrained-path method approximately handles the sign problem, which is caused by a broken symmetry in the space of Slater determinants, by eliminating any random walker as soon as . The presence of the spin-orbit interaction in the Hamiltonian means the ground state cannot be real. To ensure samples come from a real, non-negative distribution, the constrained-phase approximation generalizes the constrained-path condition: With a phase defined by

two simple forms of constrained-phase method follow Ortiz1997 () from replacing the walker by and eliminating the walker if or by which makes .

lattice with | lattice with in a magnetic field | KMH model with | |||||||
---|---|---|---|---|---|---|---|---|---|

S | S | Double occupancy | Double occupancy | ||||||

ED | -19.581 | 0.73 | 0.506 | -19.783 | -21.715 | 1.932 | 0.0030185 | -14.51 | 0.143 |

CPMC | -19.580(1) | 0.730(1) | 0.507(1) | ||||||

PCPMC | -19.580(4) | 0.731(1) | 0.508(1) | -19.779(1) | -21.712(2) | 1.933(1) | 0.0030188 | -14.35(1) | 0.148(2) |

Here we used the first constraint. Other than this difference in constraint, and the algorithmic details of the constrained-phase method are the same as those in the constrained-path method Zhang1995 (); Zhang1997 (); Carlson1999 (). When no phase problem is present, the constrained-phase method in fact reduces to the constrained-path method. When no sign problem is present, the constrained-path method is exact.

In Table I we show a comparison of the CPMC and PCPMC methods with exact diagonalization (ED) method for several different lattices (square and honeycomb), electron dopings, and spin orbit interaction strengths. Both the CPMC and PCPMC methods agree with each other and with the ED results for energies, double occupancies, and spin-spin correlations. The key point is that the PCPMC allows accurate simulations, as for two cases in Table I, not possible with CPMC or other fixed-node-like methods.

We computed superconducting correlation functions for four different pairing symmetries: (a) the NN, (b) the NNN (c) the NN, (d) and the NNN symmetries, whose form factors are illustrated in Fig. 1. These different pairing symmetries are distinguished by different phase shifts upon 60 or 120 rotations. We now define the vectors and to denote the NN and NNN inter-sublattice connections where or label the different directions. On either sublattice, the NNN-bond and wave pairings have the following form factors,

(3) |

while the singlet NN-bond pairing has the form factor

(4) |

The NN-bond is different for an A and B sublattice: for sublattice A,

(5) |

while for sublattice B, there is a phase shift (Fig. 1c)

(6) |

The NN and the NNN defined above have the same symmetry but different forms in -space and real space. Both are topological nontrivial spin triplet superconducting states. As their form factors indicate, their real (imaginary) parts are symmetric about the x-axis (y-axis) but asymmetric about the y-axis (x-axis). Thus they have a symmetry, which is simply called .

The functional form of the pairing correlation function we computed is

(7) |

where stands for the pairing symmetry and the corresponding order parameter is

(8) |

## Iii Results and discussion

Our numerical simulations are mainly performed on a lattice and a lattice, which is shown in Fig. 2. The honeycomb can be described as two interpenetrating triangular sublattices, which are marked by different colors, blue and yellow, in Fig. 2. The lattices shown preserve most geometric symmetries of the triangular Bravais lattice, and are often adopted in research on graphene. In Fig. 3, we compare the long-range part of pairing correlations for the and pairing symmetries on double-75 lattices at and . The simulations were performed for the closed-shell cases corresponding to fillings of (a) = and (b) =. For finite-size non-interacting problems, these fillings correspond to non-degenerate ground states. Simulations performed at these fillings exhibit much smaller statisical variance in the measurements than those performed at other fillings.

The first is a filling above that of the VHS; the other is below. As previous results in doped graphene without spin orbit coupling Jiang2008 (); Ma2011 (), above the VHS is larger than for most long-range distances between electron pairs. At =, while below, the NNN-bond tends to be larger than the others. Here =. For the present lattice, corresponds to a separation of lattice sites about half the cell size. From another work Ma2010 (), we know that ferromagnetic fluctuations dominate electron fillings below the VHS, and antiferromagnetic correlations dominate around half filling. Thus, the superconducting pairing appears favored by ferromagnetic correlations.

We also show results for larger spin-orbit couplings in Fig.4. For the larger spin-orbit coupling of , the correlations continue to dominate over other pairing forms at , and the correlations continue to dominate over other pairings at . One interesting point is that when compared to those in Fig. 3, the pairing correlations with symmetry are enhanced as the spin-orbit coupling increases, while the pairing correlations are suppressed, but only slightly, by increasing the spin-orbit coupling. Thus the correlations we are finding are quite robust, at least with respect to the spin-orbit interaction.

The observed dependencies of the pairing correlations on filling and spin-orbit coupling can be understood from the Fermi surface topologies. The Fermi surface with different spin-orbit couplings are shown in Fig. 5 at two electron fillings: (a) and (b). In each sub-figure the results for are indicated by a red contour and those for by a blue contour. At the low filling, in Fig. 5(a), the Fermi surface is a small circle around the point, which favors the pairing. At the higher filling in in Fig. 5(b), , the nesting vector of Fermi surface is near to the antiferromagnetic vector, which favors singlet pairing. In both sub-figures, the Fermi surface varies only very slightly as increases. The behaviors of the Fermi surfaces thus correlate with the behaviors of pairing correlations in Figs. 3 and 4, where the pairing correlations an insensitivity to the scale of spin-orbit coupling.

In order to examine the intrinsic pairing interaction in our finite system, we extracted the vertex contribution from . Doing so requires subtracting uncorrelated single-particle contributions, such as , from correlated many-body terms, such as . In Fig. 6, the vertex contribution as a function of distance for different on the double-75 lattice with and is shown. At , the vertex contribution is zero, and at , the tends to be a finite positive value over most of the long range part. Clearly, it is interesting to see that the is enhanced greatly as the interaction increases from to , over most of the long range part vertex (). Such behavior of the vertex contribution suggests effective attractions generated between electrons and the instability toward superconducting pairing in the system. The enhanced pairing strength with the enhancement of the electron correlations indicating the interaction plays a key role in inducing superconductivity. We note that the observed behavior as a function of contrasts that observed for the two-dimensional Hubbard model Zhang1995 (); Zhang1997 ().

To establish more firmly the robust presence of the triplet pairing correlations, we calculated the superconducting correlations for two more pairing symmetries, extended s-wave (Es-wave) and -wave, at low electron fillings for a larger lattice size. The phases of the pairing forms of -wave and -wave are shown in Fig. 7(a) and (b). From Fig. 7,

(9) |

and

(10) |

Figure 8 shows the various pairing correlations as a function of long-range part of the distance on a larger 2108 lattice for the closed shell fillings of (a) and with and . It is clear that the pairing correlations are larger than the other five pairings of different symmetries and this ranking does not change as the strength of the spin-orbit coupling changes. In the Hubbard model, we note that increasing the lattice size suppresses the magnitude of the pairing correlations Zhang1995 (); Zhang1997 (). We also note that in the Hubbard model near half filling it is the antiferromagnetic fluctuations that dominate.

A numerical solution for the full phase diagram for the Kane-Mele-Hubbard model is computationally challenging, especially as a function of doping. Over the range of dopings and lattice sizes studied, we do however find strong triplet superconducting correlations, and this pairing is robust against varying the spin-orbit coupling. Additionally, the effective pairing correlation is enhanced greatly as the interaction increases. Because of these results, we argue that robust spin triplet superconductivity might be present in doped honeycomb-like materials such as doped graphene, silicene, and germane or ultracold atoms in optical traps.

## Iv Acknowledgment

We thank F. Yang, H. Yao and A. Muramatsu for stimulating discussions. T. Ma thanks CAEP for partial financial support. This work is supported in part by NSCFs (Grant. Nos. 11374034, 11334012 and U1530401) and the Fundamental Research Funds for the Central Universities (Grant No. 2014KJJCB26). The work of JEG was supported by the US Department of Energy.

### References

- Kane C. L. and Mele E. J., Phys. Rev. Lett. 95, (2005) 146802.
- Sun K., Yao H., Fradkin E. and Kivelson S. A., Phys. Rev. Lett. 103, (2009) 046811.
- Liu Q., Yao H., and Ma T., Phys. Rev. B 82, (2010) 045102.
- Raghu S., Qi X.-L., Honerkamp C., and Zhang S.-C., Phys. Rev. Lett. 100, (2008) 156401.
- Wen J., Reügg A., Wang C.-C., and Fiete G. A., Phys. Rev. B 82, (2010) 075125.
- Rachel S. and Hur K. L., Phys. Rev. B 82, (2010) 075106.
- Yu S.-L., Xie X. C. , and Li J. X. , Phys. Rev. Lett. 106, (2011) 100403 .
- Assaad F. F., and Herbut I. F., Phys. Rev. X 3, (2013) 031010; Hohenadler M., Toldin F. P., Herbut I. F., and Assaad F. F. Phys. Rev. B 90, (2014) 085146.
- Hung H.-H., Wang L., Gu Z.-C., and Fiete G. A., Phys. Rev. B 87, (2013) 121113(R).
- Ma T., Hu F., Huang Z. and Lin H.-Q., Appl. Phys. Lett 97, (2010) 112501.
- Wen J., Kargarian M., Vaezi A. and Fiete G. A., Phys. Rev. B 84, (2011) 235149.
- Nandkishore R., Levitov L. S., and Chubukov A. V., Nat. Phys. 8, (2012) 158.
- Kiesel M. L., Platt C., Hanke W., Abanin D. A., and Thomale R., Phys. Rev. B 86, (2012) 020507(R).
- Black-Schaffer A. M. and Doniach S., Phys. Rev. B 75, (2007) 134512.
- Jiang Y., Yao D.-X., Carlson E. W., Chen H.-D., and Hu J. P., Phys. Rev. B 77, (2008) 235420.
- Pathak S., Shenoy V. B., and Baskaran G., Phys. Rev. B 81,(2010) 085431.
- Ma T., Huang Z., Hu F., and Lin H.-Q., Phys. Rev. B 84, (2011) 121410(R).
- Zhong Y., Zhang L., Lu H.-T., Luo H.-G., Physica B 462,(2015) 1 .
- Paiva T., Scalettar R. T., Zheng W., Singh R. R. P., and J. Oitmaa, Phys. Rev. B 72, (2005) 085123.
- Tacon M. Le, Ghiringhelli G., Chaloupka J., Sala M. M., Hinkov V., Haverkort M. W., Minola M., Bakr M., Zhou K. J., Blanco-Canosa S., Monney C., Song Y. T., Sun G. L., Lin C. T., De Luca G. M., Salluzzo M., Khaliullin G., Schmitt T., Braicovich L., and Keimer B., Nat. Phys. 7, (2011) 725.
- Scalapino D. J., Rev. Mod. Phys. 84, (2012) 1383.
- Ma T., Yang F., Yao H. and Lin H.-Q., Phys. Rev. B. 90, (2014) 245114.
- Nayak C., Simon S. H., Stern A., Freedman M., Das Sarma S., Rev. Mod. Phys. 80, (2008) 1083.
- Stern A., Lindner N. H., Science 339, (2013) 1179.
- Ortiz G., Gubernatis J. E., and Carlson J. A., Bull. Am. Phys. Soc., March Meeting M11.10 (1997).
- Schmidt K. E. and Fantoni S., Phys. Lett. B 446,(1999) 99.
- Schmidt K. E., Fantoni S., and Sarsa A., Int. J. Mod. Phys. B 15, (2001) 1510.
- Schmidt K. E., Fantoni S., and Sarsa A., Eur. J. Phys. A17, (2003) 469.
- Sarsa A., Fantoni S., Schmidt K. E., and Pederiva F., Phys. Rev. C 46, (2003) 519.
- S Zhang. and Krakauer H., Phys. Rev. Lett. 90, (2003) 136401.
- Novoselov K. S., Geim A. K., Morozov S. V., Jiang D., Katsnelson M. I., Grigorieva I. V., Dubonos S. V., and Firsov A. A., Nature (London) 438, (2005) 197; Zhang Y., Tan Y.-W., Stormer H. L., and Kim P., Nature (London) 438, (2005) 201.
- Vogt P., Padova P., Quaresima C., Avila J., Frantzeskakis E., Asensio M., Resta A., Ealet B., Lay G., Phys. Rev. Lett. 108, (2012) 155501.
- Cahangirov S., Topsakal M., Akturk E., Sahin H., Ciraci S., Phys. Rev. Lett. 102, (2009) 236804.
- Gurarie V., Radzihovsky L., and Andreev A. V., Phys. Rev. Lett. 94, (2005) 230403.
- Wu C. J., Bergman D., Balents L. and Das Sarma S. Phys. Rev. Lett. 99, (2007) 070401; Zhu S.-L., Wang B., and Duan L.-M., Phys. Rev. Lett. 98, (2007) 260402.
- Tarruell L., Greif D.,Uehlinger T., Jotzu G., and Esslinger T., Nature 483, (2012) 302.
- S. Zhang, J. Carlson, and Gubernatis J. E., Phys. Rev. Lett. 74, (1995) 3652.
- Zhang S., Carlson J. and Gubernatis J. E., Phys. Rev. B 55, (1997) 7464.
- Carlson J., Gubernatis J. E., Ortiz G. and Zhang S., Phys. Rev. B 59, (1999) 12788.
- Ortiz G., Ceperley D. M. and Martin R. M., Phys. Rev. Lett. 71, (1993) 2777.
- Chang C.-C. and Zhang S., Phys. Rev. B 78, (2008) 165101.
- Chang C.-C. and Zhang S., Phys. Rev. Lett. 104, (2010) 116402.
- Chang C.-C., Zhang S. and Ceperley D. M. , Phys. Rev. A 82, (2010) 061603(R).
- Xu J., Chang C.-C., Walter E. J., and Zhang S., J. Phys.: Condens. Matter 23, (2011) 505601.
- We do not report results for . For , we found that our estimated error for this subtraction was comparable to our estimate of the average vertex contribution. Computing the vertex contribution requires subtracting two small noisy numbers. Even for , our estimated error is around sixteen to forty percent.