Instabilities on graphene’s honeycomb lattice with electron-phonon interactions
We study the impact of electron-phonon interactions on the many-body instabilities of electrons on the honeycomb lattice and their interplay with repulsive local and non-local Coulomb interactions at charge neutrality. To that end, we consider in-plane optical phonon modes with wavevectors close to the point as well as to the points and calculate the effective phonon-mediated electron-electron interaction by integrating out the phonon modes. Ordering tendencies are studied by means of a momentum-resolved functional renormalization group approach allowing for an unbiased investigation of the appearing instabilities. In the case of an exclusive and supercritical phonon-mediated interaction, we find a Kekulé and a nematic bond ordering tendency being favored over the -wave superconducting state. The competition between the different phonon-induced orderings clearly shows a repulsive interaction between phonons at small and large wavevector transfers. We further discuss the influence of phonon-mediated interactions on electronically-driven instabilities induced by onsite, nearest neighbor and next-to-nearest neighbor density-density interactions. We find an extension of the parameter regime of the spin density wave order going along with an increase of the critical scales where ordering occurs, and a suppression of competing orders.
Electrons in graphene feature many unusual properties that can be captured within relatively simple theoretical frameworks based on a single-particle description of the electrons close to the Fermi levelnovoselov (); castroneto2009 (). The remarkable success of the interplay between experiment and single-particle theory for phenomena such as the half-integer quantum Hall effectnovoselov2005a () or the Klein paradoxkatsnelson2006a () leads to the conclusion that electron-electron interactions in pristine graphene only play a quantitative, however, not a qualitative role. On the charge neutral honeycomb lattice, due to the vanishing density of states for energies close to the Fermi level, qualitative changes from interactions such as strongly correlated electronic phases can only appear beyond a critical interaction strength sorella1992 (); khvesh2001 (); herbut2006 (). In this case, however, depending on the type of the interaction, possible occurrences of exotic states of matter such as Quantum Spin Hall phasesraghu2008 (); daghofer2013 () and even spin liquids are under considerationmeng (); sorella2012 (). Doped graphene features a non-vanishing density of states at the Fermi level which enhances the role of electronic interactions as compared to the charge neutral situation and can give rise to, possibly unconventional superconductivitydoniach2007 (); honerkamp2008 (). At least, a supercurrent in graphene has been inducedheersche2007 () by means of a contact of a graphene sample to superconducting electrodes.
This raises the question under which circumstances graphene can give rise to intrinsic superconductivity. In this context, the role of electron-phonon interactions for different types of superconducting states has been investigateduchoa2007 (); basko2008 (); kopnin2008 (); einenkel2011 () with a focus on the effects of in-plane phonons that were identified in the Raman spectra of graphenepiscanec2004 (). Further ordering patterns as, e.g., a Kekulé order due to the electron-phonon coupling have been considerednomura2009 (); kharitonov2012 (). Generally, it is a difficult task to identify the leading ordering tendency given a large variety of possible ordering patterns, in particular when various interaction effects compete.
In this work, we investigate the ordering tendencies of electrons on the honeycomb lattice when electron-phonon mediated electronic interactions from in-plane optical phonons as well as short-ranged Coulomb interactions are present. Therefore, we employ a functional renormalization group approach in a momentum-resolved patching scheme for the vertex function. This method provides an unbiased investigation of the appearing instabilities and has been proven to be a reliable tool for the study of a large range of two-dimensional solid state systems with strongly-correlated phases, e.g. high- superconductors, such as cuprates, pnictides and has been estabilished before for investigations on the honeycomb lattice, see Refs. metzner2012, ; thomale2013, for recent reviews.
This paper is organized as follows. In Sec. II, we introduce our model in terms of a tight-binding Hamiltonian with nearest-neighbor hopping and density-density interactions. Phonon-modes are included upon expansion of the hopping amplitude in the displacements and integrated out to give a contribution to the electron-electron interaction. In Sec. III, we describe the functional renormalization group (fRG) method and discuss the -patch scheme as well as the employed approximations. We present results on the ordering tendencies in Sec. IV, first by discussing exclusive in-plane optical phonons to analyze their isolated effect, see Sec. IV.1. Depending on the ratio of the strengths of the electron-phonon couplings for phonon modes close to or , we find that the leading instability is of Kekulé or nematic type. The interplay with short-ranged Coulomb interactions is studied and the impact of the electron-phonon coupling is discussed in Sec. IV.2. We draw conclusions in Sec. V.
Ii Model Hamiltonian
We consider a tight-binding model of electrons on the bipartite two-dimensional honeycomb lattice with nearest-neighbor hopping
where annihilates (creates) an electron in unit cell on sublattice with spin and analogous for sublattice . The first sum includes all neighboring sites denoted by . They are connected by the nearest-neighbor hopping amplitude which in graphene has been estimated to be eV. After Fourier transformation with and , the tight-binding Hamiltonian reads
with , where labels the primitive lattice vectors together with zero, i.e. . Explicitly, the are given by , and , where is the lattice constant. Diagonalization of gives two bands with two inequivalent, linear band crossing points at the Brillouin zone corners, the Dirac cones at and . With the spin resolved density operator we account for repulsive onsite, nearest and next-nearest neighbor interactions
An estimate for the interaction parameters can be obtained from constrained random phase approximationswehling2011 (). Diagonalizing the single-particle Hamiltonian provides an orbital makeup for the interaction terms, i.e. a momentum-dependent vertex in the band representation determined by four band indices , and three independent momenta .
ii.1 Inclusion of phonon modes
To determine the coupling of electrons and lattice displacements, we expand the hopping amplitude in the displacement fields , based on the assumption that it depends on the distance between neighboring sites, i.e.
The expansion depends on the bond direction pointing along one of the three nearest neighbor vectors. The expansion parameter is determined by ab initio calculations to be with piscanec2004 (); lazzeri2008 (); cappelluti2012 (). We introduce the phonons by the usual quantization of the Fourier transformed displacement fields using the explicit expressions and for site in sublattice and its nearest neighbor in sublattice . Further,
The carbon mass is denoted by and is the annihilation (creation) operator of a phonon in mode with momentum . Corresponding dispersions and polarizations are given by and . This inclusion of lattice distortions in terms of the phonon operators leads to the following electron-phonon coupling in orbital momentum space
Integrating out the phonons in the functional integral representation gives an effective electronic interaction
mediated by the phonons. The multiindex collects the fermionic or bosonic Matsubara frequency and the wave vector.
In the following, we will need expressions for the phonon dispersion and polarization. In principle, a calculation of the phonon spectrum would give eigenvectors with - and -components of the displacements for a given lattice site that vary with the phonon-wavector . In this two dimensional, bipartite system, the eigenvectors correspond to four possible polarizations being orthogonal to each other. In DFT calculations it was shown that the optical modes with wavevector close to and the highest-energy modes close to give the main contributions to the electron-phonon coupling strengthpark2008 (). This is why, we will concentrate on them and use only their energy and polarization in the phonon-mediated interaction, see Sec. IV. This approximation simplifies the study a lot, but should not affect the qualitative results, as the smaller variation of the phonon energy due to the dispersion of the opctial modes in the denominator of Eq. (7) does not have a strong impact on the effective interaction.
Iii FRG method
We use the functional renormalization group approach to describe the evolution from the bare action to an effective action at low energy as a function of the energy scale in an unbiased way, i.e. the structure of the effective low energy theory is not anticipated. This method accounts for effects beyond mean field and RPA as it also includes the interplay between different ordering tendencies. We use the fRG approach for the one-particle-irreducible (1PI) vertices with a momentum cutoff. For a recent review see Ref. metzner2012, . The 1PI vertices are generated by the effective action , which is the Legendre transform of the generating functional for the connected Green’s functions. In the effective action, the bare propagator of the system is modified by an infrared regulator
with the single-particle energy in band and the chemical potential . The regulator is chosen as a step function, which suppresses the modes with energy less than and reads . For numerical stability, we have slightly softened the step function in the actual implementation. With this modification of the bare propagator, we obtain a scale-dependent effective action and a variation of the scale provides a renormalization group (RG) flow. Integration of the RG flow from an initial scale (typically of the order of the bandwidth) down to low energies provides a smooth interpolation between the bare and the effective action of the system. During the lowering of the energy scale, some components of the effective two-particle interaction typically grow large and diverge at a critical scale indicating an instability towards an ordered state. The pronounced momentum structure of the near-critical interaction vertex then allows to extract an effective Hamiltonian for the low-energy degrees of freedom and determines the leading order parameter.
In the following, several approximations are employed to numerically integrate the resulting RG flow equations for the 1PI vertices efficiently: We truncate after the two-particle interaction vertex , which means that the results are of second order in the interaction. The flow of provides the essential information about the leading instabilities and we further neglect the effect of self-energy corrections as they couple to the flow of the interaction vertex only at third order, see e.g. Ref. SalmhoferHonerkamp, . In addition, we do not account for the frequency dependence of the vertex and set the external frequencies to zero to single out the most singular contribution of the flow for the determination of the ground state properties of the system. This strategy has proven to provide reliable results in a large number of different two-dimensional fermionic solid state systems, cf. Refs. metzner2012, ; thomale2013, . Despite these approximations, one obtains an infinite-order summation of second order diagrams which, importantly, accounts for the competition between different channels.
The flow equation for the coupling function reads
with labeling the Matsubara frequencies, the wave vectors and the bands. The particle-particle channel is given as
and . The direct and the crossed particle-hole channels are given by
respectively, and we define the wave vectors , and . denotes the area of the first Brillouin zone. The loop kernel reads
with the free propagator due to the neglect of the self-energy.
To solve the flow equations numerically, the wave vector dependence of the interaction vertex is discretized with a patching scheme that divides the Brillouin zone into patches, in which the vertex is approximated to be constant. This procedure has been established in Refs. SalmhoferHonerkamp, and successfully applied in a large number of works, cf. Ref. metzner2012, . The representative momentum for each patch is placed close to the Fermi level. Thereby the angular dependence of the interaction on the wavevectors is taken into account. We choose or patches as shown in Fig. 1 to check that our findings do not depend on this choice. In the actual numerical evaluation of a diverging interaction vertex, we stop the flow at a scale where the largest component of is of the order of ten times the bandwith and use this as an estimate for the critical scale . The flow to strong coupling signals that the self-energy negligence is no longer valid. Self-energy corrections would alter the low-energy spectrum, e.g. by the appearance of a gap, and without these corrections, an emerging order appears as a divergence. The critical scale can be interpreted as an estimate for the temperature below which ordering occurs or correlations of the order parameter become important.
Iv Instabilities and phase diagram
In this section, we analyze the phases of the honeycomb lattice system at temperature as a function of the bare interaction parameters and and the electron-phonon couplings that determine the impact of the phonons with wavevector close to and , respectively.
iv.1 Purely phonon-mediated interaction
We start with the study of the isolated effect of the phonon-mediated electron-electron coupling. As mentioned in the beginning, we focus on the phonons at and . This means that we set if is close to and corresponds to the optical branches. For in the vicinity of the Dirac points and labeling the three highest-energy phonons . The analogous approximation is used for the polarization vectors and all other modes are neglected. This ansatz accounts for the phonons that have been identified in the Raman spectrum of graphene and in DFT calculations to give the strongest electron-phonon coupling park2008 (); piscanec2004 (). They are often referred to as and or phonons, respectively. It has been shown in Refs. kharitonov2012, ; nomura2009, that the latter modes can give rise to a an instability corresponding to Kekulé ordering. We also find this to be the dominating instability for equally strong coupling of both modes, whereas for an enhanced coupling of the phonons, a nematic bond order is induced.
With these preliminaries and comments, we choose to parametrize the phonon-mediated contribution to the electron-electron interaction as
with due to momentum conservation and fixed polarization in as motivated before. labels the value of the phonon energy, so that the interaction parameter has units of energy111 and should not be confused with from the expansion of the hopping amplitude . With the ab-initio values for from above, the interaction parameter is of the order ..
As is well known and also discussed in the context of RG e.g. in Ref. fu2006, , the phonon-mediated interaction is suppressed for frequencies larger than the phonon frequency . Resolving this frequency dependence of the interactions in the RG flow with frequency-independent interactions requires some physical insights, at least if one wants to reduce the numerical effort222Note that in some recent fRG works honerkampfulee2007 (); uebelackerhonerkamp2012 (); giering2012 (), the frequency dependence of the interactions have been taken into account.. Usually one tries to replace the frequency dependence with a dependence on the electronic excitation energy. One reasonable choice for studying the phonon-mediated interaction case separately would then be to only include interactions of electrons with excitation energies below . This would correspond to start the RG flow only at RG scale . In this work, however, we chose to start the flow already at the bandwidth . This can be viewed as ignoring the retardation and artificially enhances the impact of electron-phonon interactions and makes its effects clearly visible.
We address several scenarios for the electron-phonon coupling by different choices for the interaction parameters and . We study the cases where only phonons from the vicinity of the point (), or only phonons close to the Dirac points (), contribute. Moreover, we include their mutual influence on each other by tuning through different ratios of with the most physical case around . The investigated parameter range spans from to . In these flows, divergences develop only if the interaction parameter is large enough, otherwise the system is a stable semimetal. This phenomenon is clearly related to the vanishing density of states at the Fermi level and has been seen for many other types of interaction-driven instabilities for fermions with this spectrum before (e.g. Refs. herbut2006, ; honerkamp2008, ; raghu2008, ). First, we discuss the results for . The discretized, bare interaction, which is the initial value of the flow equation, is shown in Fig. 2(a). Here, the critical parameter value for an instability to occur at half filling is , cf. Fig. 3(a). The momentum structure of the effective interaction close to the critical scale is presented in Fig. 2(b). It has the same structure as the bare phonon-mediated interaction, however, only for momentum transfer with and . Using this relation in the coupling function gives the effective low-energy Hamiltonian
with . For this expression, we can perform a mean-field decoupling with the molecular field yielding
ignoring the constant term.
If we compare this expression to the coupling Hamiltonian, Eq. (II.1), we recover the contribution of the -phonons from the beginning with the identification . But in order to get a diverging phonon-mediated interaction parameter , the phonon frequency must tend to zero. Thus the observed instability results in a static lattice distortion formed by the modes from , which are responsible for the Kekulé distortion shown in the inset of Fig. 3(b). Correspondingly, close to the Dirac points, the eigenenergies extracted from Eq. (IV.1) coincide with previous investigations of the Kekulé phaseChamon1999 (); Hou2007 (). In this state, the hopping is non-uniform in the three bond directions resulting in a tripled unit cell and the opening of a gap.
(b) Critical scale as function of the ratio at fixed . Insets show the nematic () and the Kekulé () state, respectively.
Now we tune the interaction to both extreme cases where one of the parameters is zero. For , we again find the Kekulé distortion as leading instability. The only difference is that the critical needed to induce the ordering is slightly decreased to , as visible in Fig. 3(a). This already shows that the small-wavector phonons controlled by have a destructive influence on the large-wavevector phonons controlled by , i.e. there is some degree of phonon-phonon interaction.
For , the behavior is qualitatively different. As shown in Fig. 3(a), the instability does not occur until a threshold value is reached and cannot be due to -phonons because they are not included in this case. Instead we find the momentum structure at low energies of Fig. 2(c), which mirrors the bare phonon-mediated interaction for zero momentum transfer . The extracted, low-energy Hamiltonian is
with and the abbreviation . The corresponding mean-field Hamiltonian results in
where is the molecular field . Comparison to the coupling Hamiltonian, Eq. (II.1) now gives the static lattice distortion due to the zone center phonons with . For , neighboring sites have different signs . This means that the two sub-lattices are moved in opposite directions in this state. The sixfold symmetry of the original lattice is reduced to a twofold one, the translational symmetry of the underlying Bravais lattice, however, is maintained, corresponding to a nematic ordering pattern. The best energy gain is a distortion along the bonds between two sites. As a result, we obtain the configuration shown in the inset of Fig. 3(b), where the hopping along one bond direction is enhanced and the Dirac points are shifted away from the Brillouin zone corners. Such a state was studied in Ref. Chamon1999, .
In Fig. 3(b), we also show the evolution of the critical scale for fixed supercritical when is increased. We clearly see that the ’Kekulé phonons’ with large wavevector transfer reduce the scale for the nematic instability and push it to zero already for . This exhibits a clear ’anharmonic’ interaction between phonon modes with different wavevectors that is revealed by the fRG treatment. When is increased further, one reaches the Kekulé-ordered phase again. Note that the rather high critical scales found here are not to be taken literally, due to the mentioned overestimation of the phonon effects when the retardation is ignored.
Usually, in more than one dimension, an important property of the electron-phonon interaction is to induce Cooper pairing, which seems to be suppressed here. We can indeed recover a conventional phonon-mediated superconducting state, but only if the RG flow equation for the interaction is reduced to the particle-particle term and all particle-hole terms are switched off. Through this the integration of the fRG equations is identical to a ladder summation in the particle-particle channel. However, the critical interaction strength needed to observe a flow to strong coupling for such an undisturbed Cooper instability is larger than in the bond-odering case. Here, without the particle-hole term, we find for . This shows that phonon-mediated superconductivity arises in the particle-particle channel as one would expect, however only when the competing contributions from the particle-hole channels are completely neglected.
iv.2 Inclusion of density-density interactions
We now also include the Coulomb-induced repulsive density-density interactions , and as given in Eq. (3). First, we include, in addition to the phonon-mediated interaction, each one of the three short-ranged interactions , and separately. This shows if the phonon-mediated interaction amplifies or weakens the effect of the respective electronic interaction. The results are compared to the case without the consideration of phonons.
Running the fRG flow with a fixed, supercritical on-site interaction for different phonon-mediated interaction strengths leads to an antiferromagnetic spin density wave (SDW) as in the case without phonons. But with increasing phonon-mediated interactions, the critical scale of the flow is enhanced. This amplifying tendency is also observed if we determine the critical on-site interaction needed to induce an instability. It reduces from to if a phonon-mediated interaction of, e.g. , is turned on, cf. Fig. 4. For the nearest-neighbor interaction, we obtain qualitatively the same behavior. As without phonons, the nearest-neighbor interaction triggers a charge density wave (CDW) whose critical scale is increased with increasing phonon mediated interaction. However, this effect is not as large as in the case of the on-site interaction. Nevertheless, the critical changes from for to for . The situation is different if we consider only a next-nearest-neighbor interaction which induces a quantum spin Hall state (QSH). Including an electron-phonon coupling suppresses the tendency for the formation of a QSH state as shown in the lower panel of Fig. 4.
These tendencies are confirmed when we run the fRG flow with all interactions, i.e. density-density repulsion up to the second nearest neighbor and the phonon-mediated interaction, included. The parameter range which we account for extends from zero through the cRPA values from Ref. wehling2011, . They are taken as upper bounds because the fRG tends to overestimate the critical scales. A summary of this investigation is given by Fig. 5.
In this work we have analyzed the impact of in-plane phonons on possible ground state orderings in a simple theoretical model for monolayer graphene. We focused on phonon eigenmodes arising from the modulation of nearest-neighbor bonds between the carbon -orbitals on the honeycomb lattice with wave vectors near and . These modes, classified as for small wavevector and , , for large wavevector transfer, are known from DFT calculations to couple most strongly to the electrons.
The electron-phonon coupling and the phonon dynamics for these modes was transformed into an effective electron-electron interaction, with some idealizations. Studying the effects of the phonon-mediated interaction without including Coulomb interactions between the electrons we find the following results. Near charge neutrality (i.e. for the undoped system) the dominant instability is in the particle-hole channel, and not in the pairing channel as is usually the case in non-nested systems. While our study may not be fully quantitative, the picture we find is that the phonons with large momentum transfer dominate in the low-energy effective interactions and that the predominant instability is towards Kekulé bond order, where the unit cell is tripled by a pattern of strengthened and weakened bonds. This state opens a gap at the Dirac points, i.e. is an insulator. Various works have argued for the existence of this state due to Coulomb interactionsnomura2009 (); kharitonov2012 (). Previous RG studies of the same model herbut2006 (); honerkamp2008 (); Scherers2012 () did not find the Kekulé order for Coulomb interactions of density-density type in the effective model, but it now occurs due to the bond-bond interactions mediated by the phonons. We can also weaken the influence of the large-wavevector phonons in the effective electron-electron interaction, emphasizing the small- phonons. In this case, a nematic instability becomes dominant where one of the three bond directions is enhanced with respect to the other two directions. The resulting spectrum features shifted Dirac points. Considering the competition between the different phonon channels, we found that the large-wavevector Kekulé phonons considerably weaken the tendencies towards nematic order driven by the small-wavevector phonons. This means that there is a significant amount of non-RPA or anharmonic physics at low energy scales, where phonons with different wavevectors interact destructively.
More realistically, the phonon-mediated interaction should be considered together with the Coulomb interaction between the electrons. The Coulomb interactions alone have been studied with RG and many other methods on honeycomb lattices in a number of works (e.g. Refs. tchougreff1992, ; herbut2006, ; honerkamp2008, ; raghu2008, ). In particular, quantum Monte Carlo calculationssorella1992 (); meng (); sorella2012 () have firmly established that the ground state for pure onsite repulsions becomes an antiferromagnetic spin-density-wave (SDW) state when the Hubbard- exceeds a threshold value. For interactions that extend further in space, only less controlled techniques are applicable. RG and saddle-point calculationsraghu2008 () found that charge-density wave states and interaction-induced quantum spin Hall (QSH) states are relevant competitors, depending on the profile of the effective interaction. Adding phonons to this interplay of the electronic ordering tendencies shifts the balance toward the SDW, while the competing QSH channel gets weakened. Interestingly, the bond phonons considered in the work actually increase the SDW and also potential CDW ordering tendencies. This can be seen most clearly in the lowering of the threshold value for the Hubbard interaction to change the semi-metal into the SDW state when the electron-phonon interaction is turned on. For the QSH state we found the reversed trend, indicating a destructive interplay with the phonons.
Hence, one important upshot of our study is the identification of the most relevant phonon-mediated effects on the ground state. Based on our study, we do not expect that the nature of potential ground state ordering or, more realistically according to the experimental state in single-layer graphene, the nature of the leading correlations if the overall interaction strength is insufficient to gap the semimetal, is determined by a phonon-mediated instability. However, we have shown that the phonon sector may actually shift the phase boundaries between different electronically-driven ordering tendencies. Hence, phonon effects may yet play an important role in deciding which of these channels wins. We can also try to extrapolate our results for the single-layer honeycomb lattice to multi-layer graphene, where experiments indeed show that the semimetal gives way to a gapped state at lower temperaturesbao (). From our theoretical experience with multi-layer honeycomb systemsScherers2012 () we can state that the main ordering tendencies are still the ones of the single layer, while the stacking only modifies the available density of states at low scales. Then, we should expect that in the multi-layer system, the phonon degrees of freedom have a similar influence in modifying the interplay of the electronic ordering tendencies as we find here. This makes the SDW state a more robust candidate to explain the observed gaps. Notably, even though the SDW state does not feature spin-resolved edge states like the QSH state, it may still be useful resource for nano-spintronic devices when the multilayer system is slightly doped and gatedyuan ().
We acknowledge discussions with L. Boeri and S.A. Maier. M.M.S. is supported by the grant ERC-AdG-290623 and DFG grant FOR 723. C.H. acknowledges support from DFG FOR 912 and SPP 1459.
- (1) A.K. Geim and K.S. Novoselov, Nature Materials 6, 183 (2007).
- (2) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
- (3) Novoselov, K. S., Geim, A. K., Morosov, S. V., Jiang, D., Katsnelson, M. I., Grigorieva, I. V., Dubonos, S. V., and Firsov, A. A., Nature 438, 197 (2005); Zhang, Y., Tan, Y.-W., Stormer, H. L., and Kim, P., Nature 438, 201 (2005).
- (4) Katsnelson, M. I., K. S. Novoselov, and A. K. Geim, Nat. Phys. 2, 620 (2006).
- (5) S. Sorella and E. Tosatti, Europhys. Lett. 19, 699 (1992).
- (6) D. V. Khveshchenko, Phys. Rev. Lett. 87, 246802 (2001).
- (7) I. F. Herbut, Phys. Rev. Lett. 97, 146401 (2006).
- (8) S. Raghu, X. L. Qi, C. Honerkamp, and S. C. Zhang, Phys. Rev. Lett. 100, 156401 (2008).
- (9) Maria Daghofer, Martin Hohenadler, Phys. Rev. B 89, 035103 (2014).
- (10) Z. Y. Meng, T. C. Lang, S. Wessel, F. F. Assaad and A. Muramatsu, Nature 464, 847 (2010).
- (11) Sandro Sorella, Yuichi Otsuka, Seiji Yunoki, Scientific Reports 2, 992 (2012).
- (12) A.M. Black-Schaffer, S. Doniach, Phys. Rev. B 75 134512, (2007).
- (13) C. Honerkamp, Phys. Rev. Lett. 100,146404 (2008)
- (14) H. B. Heersche et al., Nature(London) 446, 56 (2007).
- (15) N. B. Kopnin and E. B. Sonin, Phys. Rev. Lett 100, 246808 (2008).
- (16) B. Uchoa and A. H. Castro Neto, Phys. Rev. Lett 98, 146801 (2007).
- (17) D. M. Basko and I. L. Aleiner, Phys. Rev. B 77, 041409(R) (2008)
- (18) M. Einenkel and K. B. Efetov, Phys. Rev. B 84, 214508 (2011).
- (19) S. Piscanec et al., Phys. Rev. Lett. 93, 185503 (2004).
- (20) K. Nomura, S. Ryu, and D.-H. Lee, Phys. Rev. Lett. 103, 216801 (2009).
- (21) M. Kharitonov, Phys. Rev. B 85, 155439 (2012).
- (22) W. Metzner et al.., Rev. Mod. Phys. 84,299 (2012)
- (23) C. Platt, W. Hanke, R. Thomale, arXiv:1310.6191.
- (24) T. Wehling et al.., Phys. Rev. Lett. 106, 236805 (2011)
- (25) M. Lazzeri et al., Phys. Rev. B 78, 081406 (2008).
- (26) E. Cappelluti and G. Profeta, Phys. Rev. B 85, 205436 (2012).
- (27) J. Jiang et al., Phys. Rev. B 72, 235408 (2005)
- (28) C-H. Park et al.., Nano Letters 8,4229 (2008)
- (29) M. Salmhofer and C. Honerkamp, Progress of Theoretical Physics 105, 1 (2001)
- (30) H. C. Fu et al., Europhys. Lett. 75,146 (2006)
- (31) C. Honerkamp, H. C. Fu, and D.-H. Lee, Phys. Rev. B 75, 014503
- (32) S. Uebelacker and C. Honerkamp, Phys. Rev. B 86, 235140 (2012)
- (33) K.-U. Giering and M. Salmhofer, Phys. Rev. B 86 245122 (2012)
- (34) C. Chamon, Phys. Rev. B 62, 2806 (1999)
- (35) C.-Y. Hou et al., Phys. Rev. Lett. 98, 186809 (2007)
- (36) W. Bao et al., Nature Physics 7, 948 (2011); W. Bao et al., arxiv;1202.3212 (2012); J. Valesco et al., Nature Nanotechnology 7, 156 (2012)
- (37) M. M. Scherer et al., Phys. Rev. B 85, 235408 (2012); M. M. Scherer et al.., Phys. Rev. B 86, 155415 (2012)
- (38) A. L. Tchougreff and R. Hoffmann, J. Phys. Chem. 96, 8993 (1992).
- (39) Jie Yuan, Dong-Hui Xu, Hao Wang, Yi Zhou, Jin-Hua Gao, Fu-Chun Zhang, Phys. Rev. B 88, 201109(R) (2013)