Instabilities on graphene’s honeycomb lattice with electronphonon interactions
Abstract
We study the impact of electronphonon interactions on the manybody instabilities of electrons on the honeycomb lattice and their interplay with repulsive local and nonlocal Coulomb interactions at charge neutrality. To that end, we consider inplane optical phonon modes with wavevectors close to the point as well as to the points and calculate the effective phononmediated electronelectron interaction by integrating out the phonon modes. Ordering tendencies are studied by means of a momentumresolved functional renormalization group approach allowing for an unbiased investigation of the appearing instabilities. In the case of an exclusive and supercritical phononmediated interaction, we find a Kekulé and a nematic bond ordering tendency being favored over the wave superconducting state. The competition between the different phononinduced orderings clearly shows a repulsive interaction between phonons at small and large wavevector transfers. We further discuss the influence of phononmediated interactions on electronicallydriven instabilities induced by onsite, nearest neighbor and nexttonearest neighbor densitydensity 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.
I Introduction
Electrons in graphene feature many unusual properties that can be captured within relatively simple theoretical frameworks based on a singleparticle description of the electrons close to the Fermi levelnovoselov (); castroneto2009 (). The remarkable success of the interplay between experiment and singleparticle theory for phenomena such as the halfinteger quantum Hall effectnovoselov2005a () or the Klein paradoxkatsnelson2006a () leads to the conclusion that electronelectron 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 nonvanishing 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 electronphonon interactions for different types of superconducting states has been investigateduchoa2007 (); basko2008 (); kopnin2008 (); einenkel2011 () with a focus on the effects of inplane phonons that were identified in the Raman spectra of graphenepiscanec2004 (). Further ordering patterns as, e.g., a Kekulé order due to the electronphonon 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 electronphonon mediated electronic interactions from inplane optical phonons as well as shortranged Coulomb interactions are present. Therefore, we employ a functional renormalization group approach in a momentumresolved 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 twodimensional solid state systems with stronglycorrelated 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 tightbinding Hamiltonian with nearestneighbor hopping and densitydensity interactions. Phononmodes are included upon expansion of the hopping amplitude in the displacements and integrated out to give a contribution to the electronelectron 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 inplane optical phonons to analyze their isolated effect, see Sec. IV.1. Depending on the ratio of the strengths of the electronphonon couplings for phonon modes close to or , we find that the leading instability is of Kekulé or nematic type. The interplay with shortranged Coulomb interactions is studied and the impact of the electronphonon coupling is discussed in Sec. IV.2. We draw conclusions in Sec. V.
Ii Model Hamiltonian
We consider a tightbinding model of electrons on the bipartite twodimensional honeycomb lattice with nearestneighbor hopping
(1) 
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 nearestneighbor hopping amplitude which in graphene has been estimated to be eV. After Fourier transformation with and , the tightbinding Hamiltonian reads
(2) 
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 nextnearest neighbor interactions
(3) 
An estimate for the interaction parameters can be obtained from constrained random phase approximationswehling2011 (). Diagonalizing the singleparticle Hamiltonian provides an orbital makeup for the interaction terms, i.e. a momentumdependent 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.
(4) 
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,
(5) 
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 electronphonon coupling in orbital momentum space
(6) 
where . Similar couplings were obtained for Carbon nanotubesjiang2005 () and by a fit to abinitio values in graphenepiscanec2004 ().
Integrating out the phonons in the functional integral representation gives an effective electronic interaction
(7) 
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 phononwavector . 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 highestenergy modes close to give the main contributions to the electronphonon coupling strengthpark2008 (). This is why, we will concentrate on them and use only their energy and polarization in the phononmediated 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 oneparticleirreducible (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
(8) 
with the singleparticle 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 scaledependent 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 twoparticle interaction typically grow large and diverge at a critical scale indicating an instability towards an ordered state. The pronounced momentum structure of the nearcritical interaction vertex then allows to extract an effective Hamiltonian for the lowenergy 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 twoparticle 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 selfenergy 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 twodimensional fermionic solid state systems, cf. Refs. metzner2012, ; thomale2013, . Despite these approximations, one obtains an infiniteorder summation of second order diagrams which, importantly, accounts for the competition between different channels.
The flow equation for the coupling function reads
(9) 
with labeling the Matsubara frequencies, the wave vectors and the bands. The particleparticle channel is given as
and . The direct and the crossed particlehole channels are given by
and
respectively, and we define the wave vectors , and . denotes the area of the first Brillouin zone. The loop kernel reads
(10) 
with the free propagator due to the neglect of the selfenergy.
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 selfenergy negligence is no longer valid. Selfenergy corrections would alter the lowenergy 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 electronphonon couplings that determine the impact of the phonons with wavevector close to and , respectively.
iv.1 Purely phononmediated interaction
We start with the study of the isolated effect of the phononmediated electronelectron 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 highestenergy 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 electronphonon 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 phononmediated contribution to the electronelectron interaction as
(11)  
(12)  
(14) 
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 energy^{1}^{1}1 and should not be confused with from the expansion of the hopping amplitude . With the abinitio 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 phononmediated interaction is suppressed for frequencies larger than the phonon frequency . Resolving this frequency dependence of the interactions in the RG flow with frequencyindependent interactions requires some physical insights, at least if one wants to reduce the numerical effort^{2}^{2}2Note 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 phononmediated 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 electronphonon interactions and makes its effects clearly visible.
We address several scenarios for the electronphonon 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 interactiondriven 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 phononmediated interaction, however, only for momentum transfer with and . Using this relation in the coupling function gives the effective lowenergy Hamiltonian
(15) 
with . For this expression, we can perform a meanfield decoupling with the molecular field yielding
(16) 
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 phononmediated 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 nonuniform 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 smallwavector phonons controlled by have a destructive influence on the largewavevector phonons controlled by , i.e. there is some degree of phononphonon 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 phononmediated interaction for zero momentum transfer . The extracted, lowenergy Hamiltonian is
(17) 
with and the abbreviation . The corresponding meanfield Hamiltonian results in
(18) 
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 sublattices 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 electronphonon interaction is to induce Cooper pairing, which seems to be suppressed here. We can indeed recover a conventional phononmediated superconducting state, but only if the RG flow equation for the interaction is reduced to the particleparticle term and all particlehole terms are switched off. Through this the integration of the fRG equations is identical to a ladder summation in the particleparticle 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 bondodering case. Here, without the particlehole term, we find for . This shows that phononmediated superconductivity arises in the particleparticle channel as one would expect, however only when the competing contributions from the particlehole channels are completely neglected.
iv.2 Inclusion of densitydensity interactions
We now also include the Coulombinduced repulsive densitydensity interactions , and as given in Eq. (3). First, we include, in addition to the phononmediated interaction, each one of the three shortranged interactions , and separately. This shows if the phononmediated 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 onsite interaction for different phononmediated interaction strengths leads to an antiferromagnetic spin density wave (SDW) as in the case without phonons. But with increasing phononmediated interactions, the critical scale of the flow is enhanced. This amplifying tendency is also observed if we determine the critical onsite interaction needed to induce an instability. It reduces from to if a phononmediated interaction of, e.g. , is turned on, cf. Fig. 4. For the nearestneighbor interaction, we obtain qualitatively the same behavior. As without phonons, the nearestneighbor 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 onsite interaction. Nevertheless, the critical changes from for to for . The situation is different if we consider only a nextnearestneighbor interaction which induces a quantum spin Hall state (QSH). Including an electronphonon 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. densitydensity repulsion up to the second nearest neighbor and the phononmediated 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.
V Conclusions
In this work we have analyzed the impact of inplane phonons on possible ground state orderings in a simple theoretical model for monolayer graphene. We focused on phonon eigenmodes arising from the modulation of nearestneighbor 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 electronphonon coupling and the phonon dynamics for these modes was transformed into an effective electronelectron interaction, with some idealizations. Studying the effects of the phononmediated 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 particlehole channel, and not in the pairing channel as is usually the case in nonnested systems. While our study may not be fully quantitative, the picture we find is that the phonons with large momentum transfer dominate in the lowenergy 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 densitydensity type in the effective model, but it now occurs due to the bondbond interactions mediated by the phonons. We can also weaken the influence of the largewavevector phonons in the effective electronelectron 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 largewavevector Kekulé phonons considerably weaken the tendencies towards nematic order driven by the smallwavevector phonons. This means that there is a significant amount of nonRPA or anharmonic physics at low energy scales, where phonons with different wavevectors interact destructively.
More realistically, the phononmediated 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 spindensitywave (SDW) state when the Hubbard exceeds a threshold value. For interactions that extend further in space, only less controlled techniques are applicable. RG and saddlepoint calculationsraghu2008 () found that chargedensity wave states and interactioninduced 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 semimetal into the SDW state when the electronphonon 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 phononmediated 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 singlelayer graphene, the nature of the leading correlations if the overall interaction strength is insufficient to gap the semimetal, is determined by a phononmediated instability. However, we have shown that the phonon sector may actually shift the phase boundaries between different electronicallydriven 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 singlelayer honeycomb lattice to multilayer graphene, where experiments indeed show that the semimetal gives way to a gapped state at lower temperaturesbao (). From our theoretical experience with multilayer 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 multilayer 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 spinresolved edge states like the QSH state, it may still be useful resource for nanospintronic devices when the multilayer system is slightly doped and gatedyuan ().
Acknowledgements
We acknowledge discussions with L. Boeri and S.A. Maier. M.M.S. is supported by the grant ERCAdG290623 and DFG grant FOR 723. C.H. acknowledges support from DFG FOR 912 and SPP 1459.
References
 (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. BlackSchaffer, 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) CH. 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, DongHui Xu, Hao Wang, Yi Zhou, JinHua Gao, FuChun Zhang, Phys. Rev. B 88, 201109(R) (2013)