Modeling of gate controlled Kondo effect at carbon point-defects in graphene

# Modeling of gate controlled Kondo effect at carbon point-defects in graphene

## Abstract

We study the magnetic properties in the vicinity of a single carbon defect in a monolayer of graphene. We include the unbound orbital and the vacancy induced bound state in an effective two-orbital single impurity model. The local magnetic moments are stabilized by the Coulomb interaction as well as a significant ferromagnetic Hund’s rule coupling between the orbitals predicted by a density functional theory calculation. A hybridization between the orbitals and the Dirac fermions is generated by the curvature of the graphene sheet in the vicinity of the vacancy. We present results for the local spectral function calculated using Wilson’s numerical renormalization group approach for a realistic graphene band structure and find three different regimes depending on the filling, the controlling chemical potential, and the hybridization strength. These different regions are characterized by different magnetic properties. The calculated spectral functions qualitatively agree with recent scanning tunneling spectra on graphene vacancies.

## I Introduction

The two orbital sub-bands of a honeycomb lattice in pristine graphene realize perfect Dirac fermions Novoselov et al. (2005); Castro Neto et al. (2009) with a linearly vanishing density of states (DOS) at the Dirac point. As a result, graphene is a semimetal with a half-filled conduction band Wallace (1947); Castro Neto et al. (2009). While vacancies in metals typically act as non-magnetic scattering centres on the conduction electrons, it has been shown Wang et al. (2009); Nair et al. (2012, 2013) that removing a carbon atom induces a local magnetic moment that can undergo a Kondo screening at low temperatures Chen et al. (2011) with substantial Kondo temperatures of order . A tunable Kondo resonance in combination with the possibility to manipulate the spin degrees of freedom Kane and Mele (2005); Cho et al. (2007) may be of fundamental use for practical graphene based spintronics.

Resistivity measurements Chen et al. (2011) or the magnetic response of graphene Nair et al. (2013) with vacancies induced by irradiation provide information on ensemble averaged impurity properties while a scanning tunneling microscopy (STM) gives direct access to single vacancy qualities. Since the graphene DOS vanishes linearly at the Dirac point such a system appears to be an ideal system for observing quantum criticality Vojta et al. (2010); Fritz and Vojta (2013). While the Kondo effect is absent for arbitrary coupling at the Dirac point and particle-hole (PH) symmetry in a linearly vanishing pseudo-gap DOS Withoff and Fradkin (1990); Gonzalez-Buxton and Ingersent (1998), finite doping controlled by an external gate voltage or strong PH asymmetry can lead to Kondo screening at low enough temperatures. In addition, the Kondo physics for adatoms such as Co on graphene Vojta (2006); Rudenko et al. (2012); Fritz and Vojta (2013); Mitchell and Fritz (2013); Ren et al. (2014) has drawn some attention in recent years.

While the transition metal Co carries an intrinsic magnetic moment, the generation of spin moments at carbon vacancies in graphene is less clear. Locally, two-orbitals of the three broken sp bonds form a spin singlet bond leaving one radical with free moment Nair et al. (2013). Since these orbitals are orthogonal to the orbital subsystem responsible for the Dirac fermions, they do not couple to the conduction band prohibiting the Kondo effect. Furthermore, the carbon vacancy induces a bound state in the vicinity of the defect that is energetically located at or close to the Dirac point Pereira et al. (2007); Nanda et al. (2012); Cazalilla et al. (2012). In a tight-binding formulation for the itinerant electrons, the exact position of the bound state depends on the nearest and next nearest neighbor hopping and . A vanishing results in a symmetrical DOS and a Dirac Point that coincides with the bound state. Those bound states are orthogonal to the itinerant states forming the conduction electron continuum and it is unclear Fritz and Vojta (2013) whether these states are responsible for the experimentally observed magnetic moments Nair et al. (2013). Due to the local graphene curvature in the vicinity of the carbon vacancy, the orbitals and the neighboring orbitals start to hybridize and the possibility of a Kondo effect emerges.

In this paper, we investigate the two-orbital model for single vacancies in graphene originally proposed by Cazalilla et al. Cazalilla et al. (2012) using Wilson’s numerical renormalization group (NRG) approach Wilson (1975); Bulla et al. (2008). We calculate local spectra as functions of the gate voltage controlled chemical potential as well as the hybridization strength which is connected to graphene’s curvature in the vicinity of the vacancy. The model comprises one unbound orbital, the locally bound state and the coupling to the remaining continuum.

Mitchell and Fritz used this model Mitchell and Fritz (2013) as starting point and constructed an effective Kondo model in the local moment regime comprising of a spin-1/2 coupled to a logarithmic divergent effective conduction band density of states Cazalilla et al. (2012); Peres et al. (2009) for a non-interacting subsystem. The authors include the vacancy induced bound state in the DOS which is singular at the DP. Thus, they neglect the Coulomb interaction between and orbital and treat the problem as an effective single impurity Kondo problem. However, the divergent DOS at the DP is irrelevant for doping away from charge neutrality. Since in this paper we are interested in local fluctuations in particular and treat Coulomb interactions explicitly, we incorporate the bound state in the impurity as an effective orbital and use an effective DOS comprising only of the itinerant graphene states.

The scanning tunneling spectra (STS) Mao et al. (2017) indicate that the system is located closely to local charge degeneracy points at which the local orbital occupation is changed by one indicating that the external gate voltage not only alters the filling in the graphene layer but also the local charge configuration. Charge fluctuations become important for understanding the STS questioning the applicability of the Schrieffer Wolff transformation in the experimentally relevant parameter space. Since a density functional theory (LDA) calculation Nanda et al. (2012) predicted a substantial ferromagnetic Hund’s rule coupling between the and the orbital, combining the subsystem into a single locally projected density of states (PDOS) of non-interacting degrees of freedom appears to be an oversimplification.

Previously, single orbital Kondo and Anderson models with a graphene type pseudo-gap DOS have been studied using the NRG see also the reviews Vojta (2006); Fritz and Vojta (2013). Miranda et al. Miranda et al. (2014) focused on a disorder-mediated Kondo effect. Ruiz-Tijerina et al. Ruiz-Tijerina and Dias da Silva (2017) concentrated on the investigation of a single orbital Anderson impurity model for different geometries in the context of a dilute ensemble of atomic impurities in graphene.

The inclusion of the ferromagnetic coupling Nanda et al. (2012) between both orbitals, the orbital and the zero-mode (ZM), can favor a local triplet formation in the n doped region leading to distinct local spin configurations in different parameter spaces, and consequently to different STS responses as function of the local graphene curvature that match the recent experimental findings Mao et al. (2017) as we will show below. We have also explicitly taken into account the hybridization between the orbital and the locally bound state as included in the starting point of Ref. Mitchell and Fritz (2013). This induces a level repulsion in the two-orbital model as function of the hybridization strength that turned out to be crucial for the development of different local charge and spin configurations. We investigate the competition between these different charge and spin states as function of the chemical potential and the hybridization caused by the graphene curvature. We are able to qualitatively explain the different experimental STS spectra obtained on different carbon vacancy locations using the same two-orbital model.

The paper is organized as follows. After giving a brief overview on the experimental motivation and the STS data in Ref. Mao et al. (2017) in Sec. II, we summarize the literature on the electronic configuration around carbon vacancies in graphene in Sec. III.1 and introduce a two-orbital model including our parameterization of the local curvature. First, we focus on the charge and spin configurations of the isolated impurity in Sec. III.2 in order to set the stage for the detailed NRG calculation presented in Sec. IV. Sec. III.3 is devoted to the basics of the Kondo effect in pseudo-gap systems such as graphene. We begin the result section IV by examining each of the three relevant hybridization parameter regimes independently - Sec. IV.1, IV.2, and IV.3 - before we combine a full parameter scan of hybridization and chemical potential into a finite temperature diagram in Sec. IV.5. A detailed analysis of the Kondo temperatures is presented in Sec. IV.7. We discuss the impact of parameter changes within the model on the spectral functions as well as possible extensions of the model in Sec. IV.8 and end with a summary (Sec. V) of our findings.

## Ii Experimental Motivation

In a recent STM study Mao et al. (2017), Mao et al. have shown that carbon vacancies in graphene exhibit Kondo physics adjustable by an external gate voltage . They identified two different classes of vacancies clearly distinguishable by their characteristic scanning tunneling spectra (STS) as a function of . The gate voltage can be directly converted into a chemical potential separating the system by either p () or n () doping.

The STS of the first class of vacancies show a Kondo resonance for electron and hole doping with a vanishing close to charge neutrality as expected from the pseudo-gap DOS in graphene. For the second class of vacancies a Kondo peak is only visible in the hole doped regime and disappears upon approaching charge neutrality Mao et al. (2017).

We see the same characteristic behavior in our NRG calculations where we use the hybridization strength to distinguish both scenarios. Thus, we will refer to first scenario as the ’intermediate hybridization’ and the second as the ’strong hybridization’ regime throughout this paper. In our calculations, we also identified a third regime, where a Kondo peak is only found for . This so-called ’weak hybridization’ regime has been experimentally observed as well Mao et al. (2017). Note, that we are dealing with broad crossovers between these regimes and not sharp phase boundaries. Therefore, we will do without giving absolute values for the boundaries of the different regimes.

The goal of this paper is to give a physical explanation of all three regimes using the same two-orbital impurity model and calculate the spectral function for those classes by small parameter changes related to the local curvature of the graphene sheet at the location of the impurity.

## Iii Theory

### iii.1 Modeling vacancies by a two-orbital Anderson model

Modeling the electronic degrees of freedom in the vicinity of a graphene vacancy has a long history Pereira et al. (2008). Removing a carbon atom in a 2D monolayer of graphene leads to three dangling orbitals in the neighboring atoms. Two of them will form a new chemical bond with a doubly occupied orbital causing a Jahn-Teller like distortion of the local lattice. Electronically active remains only the third unpaired orbital Cazalilla et al. (2012) (see also Fig. 1 (left)). Charge neutrality and a weakly screened local Coulomb interaction causes a free local moment formed in this chemical radical which can undergo a Kondo screening at low temperature.

There is a long-standing debate about the possibility of Kondo screening of this moment by coupling to the conduction electrons since the Wannier states are orthogonal to the orbital in a flat 2D monolayer of graphene and, therefore, do not hybridize. This situation however changes in the presence of a local curvature which is naturally induced by either suspending the sample or by supporting it on a corrugated substrate such as SiO. Due to the local curvature, the neighboring orbitals acquire a finite overlap with the localized orbital and start to hybridize with a hybridization strength that depends on the local curvature: the larger the curvature, the stronger the coupling.

In addition, removing a carbon atom from the lattice also influences the band in the vicinity of the impurity. Treating this as a single-particle scattering problem Pereira et al. (2007) has demonstrated that now the subsystem comprises itinerant states responsible for the typical graphene density of states and one additional bound state weakly localized in the vicinity of the vacancy being orthogonal to the itinerant states. This bound state is often referred to as zero-mode since it is located at the Dirac point in the tight-binding description Pereira et al. (2007). Furthermore, density functional theory predicts a rather significant ferromagnetic coupling between the local state and the orbital of eV Nanda et al. (2012).

When setting up an effective model for the vacancy one has to carefully review the local density of states of the effective non-interacting conduction band, particularly when hunting for very subtle changes of many-body wave function due to the Kondo effect.

A first guess would be to consider the local density of state (LDOS) of the neighboring orbitals of the dangling orbital. In this LDOS Pereira et al. (2007, 2008) the bound state causes a singularity at the Dirac point. Using such LDOS is somehow misleading, since (i) the zero-mode stems for a -distribution that (ii) is associated with an interacting orbital with a significant Hund’s rule coupling Nanda et al. (2012) which is neglected in the LDOS description.

Alternatively, this bound state could be ignored by setting up a simplified single-orbital pseudo-gap Anderson model. Such a model may be sufficient to qualitatively explore the limit of strong hybridization where only occupation fluctuations between a singly and a doubly occupied orbital is relevant – see below. However, it fails to capture the physics of the weak hybridization regime without major modifications and assumptions (e.g. hybridization dependent Coulomb interaction) which have to be added by hand into the model. This simplified model also does not address the issue of the bound state: the resonance resides at yielding a singular contribution to the DOS that must be properly dealt with. Furthermore, the ferromagnetic interaction between the bound state and dangling orbital is expected to be quite strong and should still be taken into account explicitly.

Cazalilla et al. Cazalilla et al. (2012) proposed the following two-orbital Anderson model

 H=Hloc+Hπ−band+Hhyb (1)

for the description of a carbon vacancy in graphene that will be the starting point of this paper. The local two-orbital Hamiltonian is given by

 Hloc=∑σ(ϵdndσ+ϵπnπσ)+Uddnd↑nd↓+Uππnπ↑nπ↓+Udπ∑σσ′ndσnπσ−JH→Sπ→Sd−JH(d†d↑d†d↓dπ↓dπ↑+h.c.) (2)

where and denotes the single-particle energy of the orbitals with the corresponding creation operators creating an electron in the (respectively ) orbital with spin (we used the subscript for the orbital to avoid interference with the spin index) and label the intra-orbital Coulomb repulsion and the inter-orbital Coulomb interaction (see Fig. 1 (right)). In order to ensure rotational invariance in the spin space Oleś (1983) the Hund’s rule coupling has to be augmented with a pair-hopping term with the same coupling strength.

Since we have included the bound state of the system in the effective impurity, we treat the continuum by a simple tight-binding model for a bipartite honeycomb lattice as

 Hπ−band = −t∑,σ(a†jσbiσ+b†iσajσ) (3) −t′∑≪i,j≫,σ(a†jσaiσ+b†iσbjσ+h.c.)

where denotes the summation over a next-nearest neighbor hopping Wallace (1947). The operators and destroy (create) an electron on the respective sublattice A or B. Fitting of the hopping integrals to DFT calculations for the band structure yields the values and Nanda and Satpathy (2009) where the inclusion of a next-nearest neighbor hopping leads to an asymmetrical DOS and a shift of the Dirac point.

Throughout the paper, refers to charge neutrality of the system where the Dirac point coincides with the Fermi energy. In STM experiments is taken as the zero of energy and its variation with the gate induced doping is monitored by the motion of the Dirac point relative to . Likewise, the energy spectrum is symmetrically discretized around the chemical potential in the NRG. Therefore the spectral functions are also calculated with respect to as zero of energy.

The tight-binding Hamiltonian (3) can be easily diagonalized which yields two energy bands (cf. Wallace (1947); Castro Neto et al. (2009))

 ϵ±(→k) =−t′f(→k)±t√3+f(→k), (4) f(→k) =2cos(√3kya)+4cos(√32kya)cos(32kxa) (5)

where denotes the band index.

The hybridization energy takes the standard form

 Hhyb=∑→kστVd→k(c†→kστdσ+d†σc→kστ) (6)

where the influence of the band on the impurity dynamics can be fully encoded into the complex hybridization function

 Δ(z)=∑→k,τ|Vd→k|2z−ετ(→k). (7)

This defines a local Wannier orbital created by and an effective hybridization strength via

 Vc†0σ = ∑→kστVd→kc†→kστ. (8)

Taking the imaginary part yields the coupling function

 Γ(ω)=π∑→k,τ|Vd→k|2δ(ω−ϵτ(→k)), (9)

where is the one-particle tight-binding solution Eq. (4). Since , is determined by the equation

 V2 = ∑→kσ|Vd→k|2=1π∫∞−∞Γ(ω)dω (10)

that is related to the effective hybridization strength by the integral

 ∫∞−∞Γ(ω)dω=2DΓ0=πV2, (11)

where we choose for the bandwidth of graphene.

In order to capture the essence of the pseudo-gap density of states, we also used the parameterization of by the analytical expression ( pins the van Hove singularities in the DOS while determines the edges)

 Extra open brace or missing close brace (12)

The calculations presented below are performed for both hybridization functions , i.e. (i) a direct -space summation of the single-particle dispersion and (ii) the approximation (12). The different are depicted in Fig. 2.

is used as a free parameter to adjust the vacancy specific hybridization introduced above. Thus, it serves as a direct parameterization of the local curvature in the vicinity of the impurity.

### iii.2 Local scenarios and choice of parameters

In order to set the stage for the full solution of the correlated two-orbital impurity problem, we will discuss the different local configurations in a model which will become relevant for the explanation of the three different regimes found in the experiments.

In this section, we restrict ourselves to defined in Eq. (2). At charge neutrality, the orbital and the bound state should be singly occupied. In the p doped region, the orbital remains half-filled while the bound state is unoccupied and thus the total local occupation is typically close to for .

Three distinct local configurations are of particular interest for positive : A second electron can populate either the or orbital depending on the local parameters. In addition, a third electron will occupy the bound state in some parameter regimes if the additional energy gain is larger than the differences in the impurity energies. This results in a doublet state where the orbital is fully occupied and forms a local singlet whereas the electron’s spin is providing a local moment for a spin- Kondo effect.

First, consider the case . The total energy for the doubly occupied orbital (singlet state) is given by

 ES=2ϵd+Udd, (13)

which competes with the triplet state where each orbital is singly occupied and whose energy reads

 ET=ϵd+ϵπ+Udπ−JH. (14)

A local singlet ground state formed by two electrons on the orbital prohibits the Kondo effect once conduction electrons are coupled to the orbitals. On the other hand, the local triplet state still leaves the possibility for the underscreened Kondo regime with a twofold degenerate ground state from the orbital spin.

For a total orbital occupation of three electrons, the energy of the doublet state reads

 ED=ϵd+2ϵπ+Uππ. (15)

Here, the spin of the electron may also be screened by the conduction electrons showing Kondo physics. Since the state contains three electrons and, therefore, its occupation is governed by in a decoupled impurity. This state becomes favored with increasing over the singlet and the triplet configuration with and can push the system close to a charge instability.

The local scenario depends sensitively on the parameters of the local two-orbital Hamiltonian Eq. (2). The exchange interacting is crucial in order to stabilize the local triplet state. Ab-initio calculations predict a substantial ferromagnetic Hund’s rules coupling of Nanda et al. (2012) between the orthogonal state and orbital. We propose that the difference between the two classes of impurities detected experimentally and labeled as intermediate and strong hybridization regime are related to a hybridization controlled level crossing between the local singlet state with energy and the local triplet state with .

DFT calculations Yazyev and Helm (2007); Padmanabhan and Nanda (2016) and recent experimental studies Mao et al. (2017) suggest that Miranda et al. (2016). We will adopt this value but address the consequence of possible smaller in Sec. IV.8.1.

For the remaining Coulomb matrix elements, we comply with the parameters stated in the supplementary material of Ref. Cazalilla et al. (2012) where . In order to stabilize the triplet state in the strong hybridization regime, we slightly increase the onsite Coulomb repulsion on the orbital, .

Our goal is to switch between the different scenarios by just a small change in parameters that are solely caused by the variations in the local curvature of graphene in the vicinity of the different impurities. We chose the impurity single-particle energies such that the system is located close to the local crossover between singlet and triplet state outlined above.

Without a microscopic ab-initio theory for each impurity configuration at hand, we assume that in the experiment only two parameters, and , are varied: The former by the vacancy location, i. e. the local curvature, the later by the external control parameter .

The effect of the hybridization between the orbitals on the neighboring carbon sites Kanao et al. (2012) and the local orbital is twofold: (i) It couples the orbital to the band continuum and (ii) generates a hopping term between the bound state and the orbital.

A single-particle Green function approach Nanda et al. (2012) has been employed for the tight-binding model with one carbon vacancy to determine the fraction of the local orbital contributing to the bound state. The single particle wave function of this bound state decays only as . Due to the extended nature of the wave function, the factor is relatively small and has been calculated as Nanda et al. (2012). Using the hybridization parameter and diagonalizing the impurity single-particle matrix

 Hsploc = U†(ϵ′d00ϵ′π)U=(ϵd√zV√zVϵπ) (16)

with the appropriate transformation matrix determines the new eigenenergies and of the effective orbitals. Combining an estimate for and with this hybridization dependent level repulsion induced by the local curvature reduces the parameter space to and as desired.

We also have considered a simplified linear shift interpolating the level energies by the equation

 ϵ′d,π(Γ0)=ϵ′d,π(Γi0)+ϵ′d,π(Γs0)−ϵ′d,π(Γi0)Γs0−Γi0(Γ0−Γi0). (17)

Here, stands for the hybridization at which the system is either in the intermediate or strong hybridization regime. The disadvantage is obvious: There is no microscopic justification and we need to determine four reference parameters for this interpolation to arbitrary .

In the following, we omit the prime and implicitly use the shifted -dependent orbital energies. We also refer to the extended bound state as orbital in the following, since the physical orbitals have been decomposed in their contribution to and to the bound state. Since the Coulomb matrix elements are not determined by an ab-initio approach, we assume that they remain constant as function of and are always defined with respect to the final orbital-basis to keep the modeling simple.

In addition to the previous considerations, the level positions are subject to a dynamical, -dependent shift stemming from the interaction with the conduction band continuum that is explicitly included in our numerical renormalization group approach. As a consequence, the local estimations for and can only be regarded as an approximate guideline for the level positions and the charge instability.

For p doping (), the lower lying orbital is singly occupied in all cases eventually resulting in Kondo screening. The exact value of strongly depends on and . Fig. 3 summarises the different situations for positive . For small , the double occupation of the bound states becomes favorable, leaving a local moment in a ground state that can undergo a Kondo screening at large enough . Increasing will favor the spin-triplet formation in the local moment regime with that can exhibit an underscreened Kondo regime at very low temperatures. Increasing further reduces the single-particle energy enough that the doubly occupied orbital singlet state is locally favored, and the possibility of a Kondo effect is excluded.

### iii.3 Review of Kondo effect in systems with a pseudo-gap density of states

The Kondo problem has been a major point of interest in sold-state theory for multiple decades. It is a prime example of a correlated many-body problem naturally arising in metallic hosts exhibiting small magnetic dilution. At sufficient low temperature the antiferromagnetic interaction between the impurity spin and the electrons’ spin of the host gives rise to the formation of a singlet state which screens the magnetic moment of the impurity. The standard Kondo Hamiltonian Wilson (1975); Hewson (1993) takes the form

 H=∑→k,σϵ→kc†→kσc→kσ+V0∑→k,→k′,σc†→kσc→k′σ+J0→S⋅→s0 (18)

where is the Kondo coupling between the impurity spin and the local conduction electron spin density , and accounts for an additional potential scattering term arising from particle-hole symmetry breaking in the local moment (LM) regime. For a constant metallic DOS the Kondo temperature shows an exponential dependence on Wilson (1975); Hewson (1993)

 TK=√DJ0e−1/(J0ρ0). (19)

For a pseudo-gap density of states, i.e. , Withoff and Fradkin were the first to point out the existence of a critical coupling Withoff and Fradkin (1990) using a perturbative renormalization group argument. This model exhibits a wide variety of different phases for . These phases are characterized by different fixed points (FP) whose properties and occurrence depend on the bath exponent as well as on particle-hole symmetry or the absence of it. For , there exists a critical coupling strength Chen and Jayaprakash (1995); Ingersent (1996); Bulla et al. (1997); Gonzalez-Buxton and Ingersent (1998) governing the transition between a LM phase for weak coupling and a strong-coupling (SC) phase for large coupling to the metallic host.

As discussed in section III.1 and depicted in Fig. 2, the case is relevant for graphene. In the following, we briefly summarized the results of the literature Gonzalez-Buxton and Ingersent (1998); Vojta et al. (2010); Fritz and Vojta (2013). For and PH-symmetry the system will always flow to the stable LM fixed point independent of the Kondo coupling . For , Kondo screening will arise at any coupling due to the finite density of states and we expect an exponentially suppressed Kondo scale Vojta et al. (2010).

This picture changes in case of a broken particle-hole symmetry. At , screening is possible above a critical coupling and for strong enough asymmetry . Thus, there is an unstable quantum critical point separating an unscreened LM and screened ASC phase. For finite chemical potential and near the system will eventually arrive at a strong coupling fixed point but in a different way with regards to particle or hole doping Vojta et al. (2010). As a result, close to the quantum critical point (QCP) but for p doping, the Kondo temperature scales as , while for it is proportional to , where and is a prefactor of order that depends only on the bath exponent (see Fig. 4 in Ref. Vojta et al. (2010)).

### iii.4 NRG approach to quantum impurity systems

In general, we are interested in the dynamics and the thermodynamics of an interacting quantum-mechanical many-body system with an approximate continuous quasi-particle excitation spectrum. The importance of a wide range of different energy scales, from several eV on the scale of the bandwidth to arbitrary small excitations between degenerate states, gives rise to infrared divergences in the perturbation theory for the antiferromagnetic Kondo problem or Anderson impurity models. One successful approach is the renormalization group (RG) that allows for a non-perturbative treatment of all energy scales. The RG forms the foundation of Wilson’s Numerical Renormalization Group (NRG) Bulla et al. (2008) approach which has been first applied to solve the Kondo problem Wilson (1975); Krishna-murthy et al. (1980a, b).

The Hamiltonian of such a quantum impurity problem takes the three-part form of Eq. (1) where the effect of the impurity (band) is completely encoded in the term. The bath is divided into intervals on a logarithmic mesh, characterized by the discretization parameter , around the chemical potential where reconstructs analytically the original problem. One then proceeds by mapping this star geometry via a Householder transformation onto a half-infinite tight-binding chain, the so-called Wilson chain. The hopping parameter between neighboring sites and decay exponentially owing to the logarithmic discretization. The problem is now solved in an iterative fashion where the starting point is the easily diagonalizable bare impurity devoid of any bath degree of freedoms. In each iteration one additional site of the Wilson chain is included, the new eigenvalue problem is solved numerically, and high-energy excitations are discarded in order to control the otherwise exponentially growth of the many-body Fockspace. Each iteration provides access to a smaller temperature scale due to the exponentially decreasing hopping parameters. The iterative procedure is then continued until the desired temperature is reached. In the end, the set of all eigenstates and -energies make up an approximate solution to the original many-body spectrum up until the final temperature.

The NRG algorithm is not limited to equilibrium calculations and has been adopted successfully to non-equilibrium problems Anders and Schiller (2005, 2006). The identification of a complete basis set is the foundation for a sum rule conserving formulation of the impurity Green’s function Weichselbaum and von Delft (2007); Peters et al. (2006) and eventually for present non-equilibrium extensions. Equilibrium NRG results provide an almost exact agreement with analytical results if available and often serve as a benchmark for other methods. For a more in-depth examination refer to state-of-the-art reviews such as Bulla et al. (2008).

## Iv Results

This section starts out with presenting our results for the three different generic scenarios that are characterized by before we compile everything into a finite temperature as well as a regime diagram. These diagrams will be presented in Sec. IV.5. Fig. 3 which summarizes the local configurations for positive chemical potential and three different generic values of serves as a road map. The scenarios represent the different classes of vacancies identified by STM experiments Mao et al. (2017).

We used a discretization parameter and kept states after each iteration in all our NRG calculations. The calculations were done for finite temperature unless stated otherwise to be comparable to experimental results.

We employed three different approximations for and the single-particle orbital energies: (i) the approximate stated in Eq. (12) and a simple linear interpolation, (ii) the realistic generated from the dispersion and a linear interpolation of , and (iii) the realistic as well as a re-diagonalization based on the factor for the bound state. In all cases, the single particle energies and are functions of the selected .

The parameters used in the NRG calculations are mainly extracted from different sources in the literature, as discussed in Sec. III.2 and as summarized in table 1. The question whether different sets of parameters, especially smaller values for , results in similar experimental findings is addressed in Sec. IV.8.1. In order to reach the desired temperature, we choose and the number of NRG iterations is given by .

If not stated otherwise, the spectral functions are for the orbital only, i.e. . The effect of the contribution to the spectrum is discussed in Sec. IV.4.

### iv.1 Weak hybridization regime

Fig. 4 shows the calculated spectral functions Weichselbaum and von Delft (2007); Peters et al. (2006) for the orbital for a weak hybridization (). The results are depicted in panel Fig. 4(a) for and those for in panel Fig. 4(b) and cover a range from to . For clarity, the spectral functions are shifted by a constant for each and normalized by being the maximum of all spectral functions. While remains featureless in the p doped regime and reflects the pseudo-gap DOS of the conduction band, Kondo resonances are visible in the spectra for .

In order to relate the NRG results to the scenarios presented in Fig. 3, we plot the local orbital occupations as function of in Fig. 5(a). In the weak hybridization regime, the lower orbital remains half-filled and nearly independent of while the occupation of the orbital depends on the chemical potential. For the orbital is also half-filled due to the energy gain by the Hund rule coupling forming a local triplet. Increasing to adds another electron to the orbital since the interaction is weak and a local moment with remains on the orbital. Therefore, we identify an underscreened (undercompensated Cox and Zawadowski (1998)) Kondo problem for , while we find a conventional pseudo-gap Kondo problem for since . The observation of a lower for a underscreened Kondo problem in recent NRG calculations (see Ref. Roch et al. (2009)) is compatible with our findings where for .

Additionally, the influence of the chemical potential on the Kondo temperature is asymmetric even for a fixed Kondo problem Vojta et al. (2010); Fritz and Vojta (2013) with respect to the sign change of .

### iv.2 Intermediate hybridization regime

Increasing moves the system into the intermediate hybridization regime. Here, a Kondo peak is found in for both strong p doping or n doping as clearly depicted in Fig. 6. The spectral functions are normalized as before for clarity. For small bias voltage, the Kondo effect breaks down since vanishes linearly at suppressing screening close to the Dirac point. The Kondo temperature as function of is extracted using a Goldhaber-Gordon fit Goldhaber-Gordon et al. (1998) of the zero bias conductivity

 G(T)=G0[1+(21/s−1)(T/TK)2]s, (20)

where and . shows an exponential dependence on the chemical potential close to , albeit with different slope for p and n doping (not shown).

We have chosen the impurity parameters in the intermediate hybridization regime in such a way that the occupation of the orbital changes around Nanda et al. (2012) as depicted in Fig. 5(b). The orbital remains roughly half-filled as function of while the occupancy of the orbital changes rapidly from zero to single occupation. There is only a very small difference in the results obtained from the approximated , Eq. (12), and the realistic hybridization function emphasizing the generic character of the power law approximation of .

Since the local occupancy changes from to , the underlying Kondo effect is fundamentally different for both types of doping. For the local spin moment is screened by the bath forming a conventional Kondo singlet. For positive , however, the two spins form a local triplet state due to the strong ferromagnetic Hund’s rule coupling . Only a fraction of the resulting moment can be screened by the bath resulting in an underscreened Kondo effect and a dangling moment Cox and Zawadowski (1998).

### iv.3 Strong hybridization regime

When is increased further compared to the intermediate regime, eventually the strong hybridization regime is reached.

For , it essentially resembles the intermediate regime. The single occupied orbital provides an unpaired spin that is screened to a Kondo singlet, and shows a Kondo resonance (Fig. 7(a)). The orbital, however, remains mostly empty for all chemical potentials considered here and the physics is entirely governed by the orbital; this regime could also be described by a single orbital Anderson model Mao et al. (2017).

Starting around and increasing , the orbital gets slightly filled such that the total occupation approaches – see also Fig. 5(c). Therefore, the system approaches the intermediate valence (IV) FP governed by local charge fluctuations Krishna-murthy et al. (1980b); Lo et al. (2014). A resonance develops in close to the DP that moves towards lower energies as increases as shown in Fig. 7(b). This coincides with an increase of the orbital occupation. The Kondo resonance for evolves continuously to a charge fluctuation resonance close to since the local and charge configurations become energetically degenerate for . While a doubly occupied orbital is predicted in a purely local picture, the substantial hybridization favors the itinerancy of the impurity electrons. The local occupation of the orbital is shown in Fig. 5(c) and approaches indicating the IV FP.

### iv.4 Zero-mode peak

We only dealt with the spectral functions for the orbital in the previous sections. However, the tunneling STS current measured in the experiments Mao et al. (2017) results from a superposition of three parts: The orbital, the orbital as well as the substrate.

The orbital is responsible for the Kondo physics but the orbital contributes to the total experimental curves in form of the zero-mode. The orbital occupation Fig. 5 gives a rough estimate for the position of the zero-mode peak relative to the Fermi energy in the different regimes: a vanishing occupation implies a zero-mode that rests above . We used the same parameters as in the previous sections. The spectral functions depicted in Fig. 8 are calculated by Lehmann representation of the NRG data Weichselbaum and von Delft (2007); Peters et al. (2006) using a logarithmic Gaussian broadening – see the NRG review for details Bulla et al. (2008).

The -spectral functions in the intermediate hybridization regime, , are shown for varying the broadening parameter Bulla et al. (2008) and in Fig. 8(a). This demonstrates the very sharp nature of the ZM excitation whose width is artificially broadened by . In addition, a Hund’s rule mediated excitation at slightly higher energies emerges for small which results from the exchange coupling between the orbital and the orbital.

The zero-mode peak exhibits a simple dependence as shown in Fig. 8(b): A higher chemical potential shifts the peak towards . This is the exact same behavior observed in the experiments (Fig. 2(a) in Mao et al. (2017)).

Since we are mainly interested in the Kondo physics which is only contained in the orbital spectral function we will not discuss the spectral function in the following.

### iv.5 Mapping out the parameter space

#### Regimes at T=4.2K

After presenting the orbital spectral functions for three different characteristic values for representing the three different classes of vacancies found in the STM experiments, we combine the results into a more elaborate scan of from weak to strong hybridization in a single regime diagram where and are the only free parameters.

We used the realistic hybridization calculated from the DOS and the single-particle orbitals energies calculated for a given via Eq. (16) in all NRG calculations in this section. For each data point , we obtained the impurity entropy Bulla et al. (2008) and total orbital occupation from an individual NRG run. The calculations were performed at a fixed finite temperature to make contact to the STM experiments. The diagram depicted in Fig. 9 shows the color-coded product of the residual entropy times the total impurity occupation.

In Fig. 9 three different areas are separated by two naturally appearing lines that define a sharp, but temperature broadened crossover region characterized by . In the left part of the figure, rapidly drops from three to two, while simultaneously the entropy increases from to . The second line occurs where and .

We identify seven different regimes in Fig. 9. Increasing in the p doped region (), we found three local moment and one strong coupling regime: (LM-a), (LM-b), (LM-c), and (SC-a). For positive , we find two strong coupling FPs ( SC-b and SC-c) as well as a frozen impurity regime (FI). Note that the LM regimes mentioned above also extend partially to due to the vanishing at the Dirac point.

For the smallest hybridization in the range of , the orbital is doubly occupied for positive , and we find a half-filled orbital. This regime is labeled as the LM-a since the entropy tends towards and the square effective magnetic moment Wilson (1975) is given by . This regime extends into p doping and crosses over into the SC-c regime upon simultaneous increase of and .

Increasing at a fixed enhances the level splitting between the two local orbitals: the doubly occupied orbital becomes less favorable and the system crosses over to the (LM-b) regime. The thick red line marks the crossover from a local occupation of to which also distinguishes LM-a and LM-b regimes. Thus, the entropy in the LM-b regime tends to while the effective moment takes the value . The ground state is a triplet state formed by both electrons which are coupled due to the strong ferromagnetic interaction .

Increasing the hybridization further delocalizes another impurity electron, and the system crosses over to the LM-c region where only a single orbital spin is locally present. The spectral function already shows the onset of a Kondo peak (see Fig. 6) even though the entropy and the effective moment still signifies a LM FP: and . The Kondo temperature is significantly increased above at very larger and , and the crossover from LM-c to the SC-a regime is observed. The entropy as well as the square effective magnetic moment tends to zero due to the formation of a Kondo singlet state in this regime.

Now we discuss the diagram for starting from the upper right area of Fig. 9. For strong enough hybridization and large , a doubly occupied orbital and an empty orbital are locally favored. Although that regime is labeled as frozen impurity (FI) regime, it is closer related the IV FP since . Interestingly enough, the FI ground state calculated by the NRG corresponds to a double occupied state while the impurity occupation is closer to . Since the ground state is a local singlet state, the impurity decouples effectively from the conduction band. The coupling to the conduction band, however, leads to a gain of kinetic energy and partially delocalizes the orbital electrons. The spectral function shows a pronounced excitation peak which shifts linearly with the chemical potential as depicted in Fig. 7(b).

Upon the reduction of , the local spin triplet state becomes energetically favorable. As in the p doped region, we see a broad crossover regime where entropy and magnetic moment still take their LM fixed point values although an underdeveloped Kondo peak is already visible – see Fig. 6(b) at . Both local electrons align ferromagnetically, and the hybridization is sufficient to initiate Kondo screening of the electron spin. For sufficiently large , the underscreened Kondo fixed point (SC-b) is found with a dangling electron spin. The SC-b regime is characterized by and impurity entropy of and square effective moment of .

At small – upper left corner of Fig. 9 – we observe another crossover from the LM-a to the strong coupling regime SC-c. Both regimes are characterized by a fully occupied orbital. Therefore, we end up with an effective isolated half-filled orbital which may now be screened by the conduction band’s electrons when is increased with increasing .

Additionally, we added three areas to the color contour plot Fig. 9 to indicate the three regimes with their distinct spectral evolution as function of : (i) the weak, (i) intermediate, and (iii) strong hybridization regime.

#### Fixed points T→0

Away from the line, all LM regimes are associated with unstable RG FPs. The diagram in Fig. 10 is calculated for the same parameters as Fig. 9 but for . Clearly, the LM-c regime characterized by a single orbital spin is almost completely suppressed and reduced to a narrow region in the vicinity of where a rapidly vanishing exponentially reduces . This regime has been replaced by the stable SC-a FP for and the FI FP for . The area of two other stable SC FPs, SC-b and SC-c, have also increased but the two other LM FP, LM-a and LM-b still cover a large area of the diagram. Further reduction to – not shown here – reduce these regions in favor of the corresponding SC FPs.

### iv.6 Charge crossover from single to double occupation and phase transition at T→0

The question arises whether the crossover from LM-c to LM-b in the intermediate regime upon increasing will develop to a quantum phase transition at . Fig. 11 shows the spectral functions, panel (a), as well as the orbital occupation, panel (b), in the crossover region.

Already at finite the crossover takes place in a narrow range between and giving rise to the clearly visible distinction in Fig. 9. The spectral functions for are characterized by a peak at whose spectral weight is shifted towards slightly higher energies for .

The orbital occupation of the impurity illustrates the change from a crossover to a quantum phase transition at . The low lying orbital remains at around half-filling regardless of temperature and . On the other hand, the orbital changes its occupation from empty to half-filled in a smooth manner at finite . This develops into a sharp transition for indicating a real phase transition and a change of ground states.

We see two lines of quantum critical points for our parameters, where the impurity changes its total occupation by an integral number.

### iv.7 The Kondo temperature

In section IV.5 we have shown that the system has not yet reached its stable fixed points for a large part of the parameter space at finite . The broad crossover regimes show qualities which can be ascribed to either a LM or SC phase, such as an entropy while simultaneously exhibiting an onset of a Kondo peak. Therefore, estimating the Kondo temperature by using the full width half maximum (FWHM) of the peak may result in values that differ significantly from the true low energy scale for the given parameter sets. In addition, the PH asymmetry and the fitting procedure does not uniquely determine Esat et al. (2015).

We adopted a twofold strategy: (i) we estimated using a Fano line shape for the Kondo peak in the spectral function in order to connect to the STS approach for extracting from spectra experimental obtained over a limited temperature range, and (ii) we calculated from the temperature evolution of the zero-bias conductance (ZBC) originally proposed by Goldhaber-Gordon Goldhaber-Gordon et al. (1998) in the context to quantum dots that turns out Esat et al. (2015) to coincide very accurately with Wilson’s definition Wilson (1975). However, this method for extracting might be questionable in the case of STS which shows Fano resonances Schiller and Hershfield (2000).

Fig. 12 shows the calculated values using both approaches. Note that at finite in general the ZBC has not yet reached its plateau value rendering the Goldhaber-Gordon approach useless, so that we iterated until the fixed point is reached. For all calculations we use that stems from the realistic DOS as well as local orbital energies obtained from , Eq. (16).

The first two curves – from top to bottom – are obtained for hybridization and can be identified with the strong hybridization regime, the next two () show the occurrence of a Kondo peak characteristic for the intermediate regime, and the last one belongs to the weak hybridization regime.

Note that the values of determined by the ZBC fit correspond to the crossover temperature for the entropy change towards the SC depicted in Fig. 9. The agreement between extracted from the ZBC fit and true thermodynamic energy scale characterizing the crossover to the stable FP in the NRG has been previously observed in the analysis of STS for Au-PTCDA complexes on Au surface Esat et al. (2015).

### iv.8 Modification of the model

All features observed experimentally and reported in Fig. 4(a) by Mao et al. Mao et al. (2017) are in qualitative agreement with Fig. 12. However, we also note some differences to the STS data present in Ref. Mao et al. (2017): (i) calculated by the two-orbital model is about a factor of too small and (ii) the Kondo resonance in the weak and intermediate regime sets in at around at the earliest for finite . Since the Kondo temperature is exponentially sensitive to the parameters of the model and the determination by a FWHM fit in the experiment is not reliable – see discussion in Ref. Esat et al. (2015) – the difference between the experimental and the theoretical Fano fit is clearly not as significant as the different onset of the Kondo effect for .

#### Effect of a smaller on-site U

The exact value of the intra-orbital Coulomb interaction is not agreed upon in the literature with predictions spanning roughly one order of magnitude Yazyev and Helm (2007); Nair et al. (2013); Wehling et al. (2011); Vergés et al. (2010). In the light of a recent study, we investigate the effect of a reduced Coulomb matrix element as proposed by Nair et al. Nair et al. (2013). A smaller Coulomb interaction enhances local charge fluctuations and presumably can enhance using the Schrieffer-Wolff transformation Schrieffer and Wolff (1966) as a guideline. A larger on the other hand suppresses those charge fluctuations and, therefore, requires an even larger, possible non-realistic hybridization .

Note that is particularly small in comparison with the Hund’s rule coupling obtained by density functional theory Nanda et al. (2012) and the inter-orbital interaction . In order to guarantee the occupation necessary for developing a stable local moment, we slightly increased to .

The small on-site had to be combined with a value of in the intermediate (strong) hybridization regime for the system to form a stable moment in order to find a Kondo resonance at the chemical potential. Furthermore, a smaller tended to result in a vacant impurity ground state abruptly destroying the Kondo effect.

We restricted ourselves to the most approximated (Eq. (12)) for the calculations in order to focus an the qualitative differences to . Fig. 13 shows for the modified parameter set at . The intermediate (solid lines) and strong hybridization (dashed lines) regimes are clearly visible with a pronounced Kondo peak. However, the width of the Kondo resonance is extremely small corresponding to a Kondo temperature not comparable to recent experimental STM data (Ref. Mao et al. (2017)).

Our calculations show that both regimes appear to be a generic feature of our two-orbital pseudo-gap model. Counterintuitively, higher Coulomb repulsion yields a significant higher Kondo temperature and better agreement with STM experiments but requires an increase of . We could increase and find adequate parameters for , but it is not clear whether the corresponding large hybridization matrix element can be justified by a curvature induced overlap between the local orbital and the neighboring tilted orbital, when the orbital matrix elements entering the band structure are given by .

#### Possible modifications of the model by neglected interactions

In the literature as well as in our calculations, the starting point is always a perfect ideal flat graphene sheet. The local modification by a rippled and curved graphene is only included in (i) the shift of the local Dirac point and by the finite orbital overlap resulting in a finite . We recall that the local curvature is a consequence of an energetically favored 3D over a purely 2D graphene sheet according to the Mermin-Wagner theorem and the underlying SiO substrate. This will also modify the phonon modes. Long-wave length Goldstone phonons with vanishing phonon energies will be present in the material but can be neglected for our problem. Locally, however, there will be breathing modes of longitudinal vibrations of the curved graphene sheet that might expand and contract slightly the distance of the graphene and the STM tip. Clearly such a local vibration will have a finite frequency similar to a molecular vibration.

If we assume the existence of a single local vibrational breathing mode it will have two effects: (i) it will modify the Hamiltonian (1) describing the impurity dynamics and (ii) it will modify the tunnel theory employed to interpret the STS data. The definition of in Eq. (1) has to be augmented by the two additional terms in

 Hph = Ω0b†b+λph(b+b†)Hhyb (21)

where annihilates a phonon with energy . The first term accounts for free phonons while the second term result from the expansion of the local hybridization in Eq. (8) by

 V(X) = V(1+λphX) (22)

with the local displacement operator and the dimensionless electron-phonon coupling strength .

Such an additional term has been investigated within the NRG in the context of the quantum transport through a deformable molecular transistor Cornaglia et al. (2005). In the limit of weak electron-phonon coupling, it was analytically shown that the Kondo temperature will increase due to an enhancement of the Kondo coupling. An investigation of the impact onto in the pseudo-gap two-orbital model as well as an estimation of inelastic tunnel processes is desirable but beyond the scope of this paper.

In the context of charge transport through a two-orbital molecule the influence of an electron-phonon coupling onto an additional spin anisotropy term in the spin-1 sector was investigated Cornaglia et al. (2011); Ruiz˘Tijerina et al. (2012). This spin anisotropy requires a spin-orbit coupling which is relevant in Co complexes but weak in carbon. The effect of a small spin anisotropy can be enhanced by a linear spin-lattice coupling in the anti-adiabatic limit resulting in a reduction of . Due to the weak spin-orbit coupling in carbon, we do not consider this effect to be of relevance for the spin dynamics of graphene vacancies. Furthermore, the zero-mode bound state is extended Nanda et al. (2012) and can only couple very weakly to a local phonon.

#### Anisotropic hybridization V

As indicated in Fig. 1, the dangling sp orbital hybridizes with both orbitals located at the base of the distorted triangle Kanao et al. (2012). The hybridization part of the Hamiltonian can be written as

 HHyb=∑σ(V1a†1σdσ+V2a†2σdσ+h.c.) (23)

where creates an electron in the orbital at the nearest-neighbor position relative to the vacancy. The localized orbital resides at . We consider the general case where , i.e. we allow the hybridization to differ in phase and absolute value. A point-defect distorts the non-planar geometry of graphene and may lead to a symmetry breaking hybridization.

Combing the conduction band operators in (23) to the effective operator leads to the general hybridization function

 ~Γ(ω) = 1N∑→kτ(|V1|2+|V2|2+2R(V1V∗2ei→k(→δ1−→δ2))) (24) ×δ(ω−ϵτ→k) .

For , we can follow the approaches of Refs. Jones and Varma (1987); Affleck et al. (1995); Kanao et al. (2012) and combine both orbitals that are situated opposite the free sp (or ) orbital to arise at a pair of new orthogonal operators (even and odd parity). The fermionic commutator determines the new effective DOS and hybridization function for even () and odd () parity respectively. The corresponding hybridization functions take the form

 Γ±(ω)=4|V1|2N∑→k,τδ(ω−ϵτ→k)⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩cos2(→k(→δ1−→δ2)+φ2−φ12)sin2(→k(→δ1−→δ2)+φ2−φ12). (25)

Here, is the dispersion for the energy band given by Eq. (4). In the case where both orbitals hybridize equally to the free sp orbital, only the even linear combination couples to the impurity Kanao et al. (2012) while the odd parity linear combination decouples from the sp orbital. The authors of Ref. Kanao et al. (2012) restricted themselves to a linear dispersion around both Dirac points and a fully symmetric hybridization in which case Eq. (25) is analytically solvable for and becomes proportional to the zeroth order Bessel function.

Here, we make use of the full tight-binding dispersion including second nearest neighbors, Eq. (4), and evaluate the summation over the first Brillouin zone explicitly. In the presence of a vacancy the difference can be approximated by the lattice constant Kanao et al. (2012). Both hybridizations, , are shown in Fig. 14. The characteristic linear dependence close to the Dirac Point is preserved for both linear combinations. However, the even linear combination would suppress the Kondo effect significantly for a fixed due to the reduced slope. In order to retain Kondo temperatures close to the experimental ones, one could simply increase which would eventually lead to questionable high values for . A relative phase shift of would result in a swap between even and odd hybridization and would increase . In order for , either the configuration belongs to a local odd parity, or the parity is broken in the vicinity of the vacancy: both orbitals are slanted differently in relation to the -plane.

In order to reveal the effect of the different energy dependence of hybridization functions for the same hybridization strength – see Eq. (11) — we present a direct comparison of a spectral functions from the intermediate regime at for based on a realistic DOS used throughout this paper and the even/odd hybridization as defined in Eq. (25) in Fig. 15 using identical local impurity parameters. While the Kondo temperature is strongly suppressed for the even , we observe an enhanced Kondo temperature for illustrating the effect of the significantly large slope of the region in the vicinity of the Dirac point. Note that, in order to compare the consequence of the changed hybridization throughout, the parameters of the impurity ought to be adjusted.

In a flat graphene sheet with parity conservation at the location of the orbital, the even hybridization function should be relevant Kanao et al. (2012). Due to the orthogonality of the Wannier orbitals, however, a hybridization is only generated for curved graphene in the vicinity of the carbon vacancy. This locally distorted environment breaks parity in general. The matrix elements will be different in modulus and phase, and the most general hybridization function is given by

 ~Γ(ω) = |¯V|2N∑→kτδ(ω−ϵτ→k) (26) ×(1+2|V1||V2||V|2cos(→k(→δ1−→δ2)+δφ))

defining . If , close to the Dirac point. Since the parity breaking should be taken into account and a detailed microscopic theory of the specific single-particle properties in the vicinity of a carbon vacancy is missing, we restricted ourselves to the defined in Eq. (9).

## V Summary and conclusion

We investigated an effective two-orbital quantum impurity model for a carbon vacancy in locally curved monolayer of graphene. Combining the estimations for the local impurity parameters Yazyev and Helm (2007); Cazalilla et al. (2012); Nanda et al. (2012) with the experimental STM data Mao et al. (2017) we identified three different regimes, namely the weak, intermediate, and strong hybridization regime. Our NRG calculations are carried out for matching the experimental temperature Mao et al. (2017). The intermediate regime is characterized by SC fixed points for both n and p doping which crosses over into unstable LM regimes for a constant temperature when approaches the Dirac point as the vanishing DOS of graphene exponentially suppresses screening of the impurity spins. The other two regimes only show a Kondo effect in the p doped region for large hybridization (curvature) or in the n doped region for small hybridization (curvature). For positive and significantly strong hybridization the systems is driven into a FI regime whose striking feature is a single peak in the spectral function shifting away from with increasing .

All three regimes are connected by changing the graphene curvature that also changes the effective single particle orbital energies: Not only the conduction band hybridizes with the orbital but also the vacancy induced bound state. This physical mechanism driving the transition from weak to intermediate and finally to the strong hybridization regime is related to the local rippling of graphene that has been verified experimentally Mao et al. (2017). All three regimes can be understood in terms of local spin and charge configurations of a two-orbital impurity model.

Our findings are in qualitative agreement with recent experimental studies by Mao et al. Mao et al. (2017) where scanning tunneling microscope (STM) measurements on irradiated graphene sheets found all three scenarios albeit with higher Kondo temperatures. We have shown that all three regimes are an intrinsic feature of our realistic two-orbital model and are not necessarily dependent on the exact form of the DOS as long as the pseudo-gap slope for is present. Even for particularly small Coulomb interactions we can still reproduce the intermediate and strong hybridization regime although is diminished substantially.

The only drawback of our approach is the relatively low Kondo temperature compared to the experiment. In order to predict larger s, a substantial local hybridization is required Vojta et al. (2010); Mitchell and Fritz (2013) whose microscopic origin is unclear: Note that the nearest neighbor hopping matrix element in pristine graphene has been approximate , and hybridization matrix elements of the same order are required to achieve s of the order of . This might be caused by several aspects not included in our model. Firstly, in all calculations we assume that the Dirac Fermion DOS prevails in the vicinity of the vacancy up to the bound state. However, Kondo temperatures depend sensitively on the local close to the Dirac point. The true energy dependence of in a curved environment – a result of the vacancy, the substrate, and the high gate voltages – might differ slightly from that of pristine graphene. Secondly, we assumed that the STS is directly proportional to the spectral function of the -orbital while in reality, the STM tip can locally couple to several different orbitals Schiller and Hershfield (2000). The additional asymmetries of such a resulting spectrum stemming from the interference between multiple transport channels, one to the local orbital and the others to the graphene p band, might lead to a larger value of by a two-parameter Fano fit. Thirdly, local rippling of graphene is caused by an instability of the ideal 2D graphene sheet. Local curvatures will modify the phonon spectrum, and it is expected that there is a significant electron phonon coupling contributing to the hybridization. This electron phonon interaction has been proven to enhance the effective hybridization Cornaglia et al. (2005) and therefore . In addition, one expects inelastic contributions to the tunnel current if such a breathing mode leads to an oscillation of the distance between the vacancy and the STM tip. Such an additional contribution could modify the overall shape of STS and therefore also change the FWHM used to estimate .

###### Acknowledgements.
We thank R. Bulla and M. Vojta for fruitful discussions. We acknowledge support from the Deutsche Forschungsgemeinschaft via project AN-275/8-1 (F.B.A. and D.M.), DOE-FG02-99ER45742 (E.Y.A. and J.M.), NSF DMR 1708158 (Y.J.), Ministry of Science and Technology and also Academia Sinica of Taiwan (G.Y.G. and P.W.L.).

### References

1. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos,  and A. A. Firsov, Nature 438, 197 (2005).
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. P. R. Wallace, Phys. Rev. 71, 622 (1947).
4. Y. Wang, Y. Huang, Y. Song, X. Zhang, Y. Ma, J. Liang,  and Y. Chen, Nano Letters 9, 220 (2009).
5. R. R. Nair, M. Sepioni, I.-L. Tsai, O. Lehtinen, J. Keinonen, A. V. Krasheninnikov, T. Thomson, A. K. Geim,  and I. V. Grigorieva, Nat Phys 8, 199 (2012).
6. R. R. Nair, I.-L. Tsai, M. Sepioni, O. Lehtinen, J. Keinonen, A. V. Krasheninnikov, A. H. Castro Neto, M. I. Katsnelson, A. K. Geim,  and I. V. Grigorieva, Nature Communications 4, 2010 EP (2013), article.
7. J.-H. Chen, L. Li, W. G. Cullen, E. D. Williams,  and M. S. Fuhrer, Nat Phys 7, 535 (2011).
8. C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
9. S. Cho, Y.-F. Chen,  and M. S. Fuhrer, Applied Physics Letters 91, 123105 (2007).
10. M. Vojta, L. Fritz,  and R. Bulla, EPL (Europhysics Letters) 90, 27006 (2010).
11. L. Fritz and M. Vojta, Reports on Progress in Physics 76, 032501 (2013).
12. D. Withoff and E. Fradkin, Phys. Rev. Lett. 64, 1835 (1990).
13. C. Gonzalez-Buxton and K. Ingersent, Phys. Rev. B 57, 14254 (1998).
14. A. N. Rudenko, F. J. Keil, M. I. Katsnelson,  and A. I. Lichtenstein, Phys. Rev. B 86, 075422 (2012).
15. A. K. Mitchell and L. Fritz, Phys. Rev. B 88, 075104 (2013).
16. J. Ren, H. Guo, J. Pan, Y. Y. Zhang, X. Wu, H.-G. Luo, S. Du, S. T. Pantelides,  and H.-J. Gao, Nano Letters 14, 4011 (2014), pMID: 24905855, http://dx.doi.org/10.1021/nl501425n .
17. V. M. Pereira, J. Nilsson,  and A. H. Castro Neto, Phys. Rev. Lett. 99, 166802 (2007).
18. B. R. K. Nanda, M. Sherafati, Z. S. PopoviÄ,  and S. Satpathy, New Journal of Physics 14, 083004 (2012).
19. M. A. Cazalilla, A. Iucci, F. Guinea,  and A. H. Castro Neto, ArXiv e-prints  (2012), arXiv:1207.3135 [cond-mat.str-el] .
20. K. G. Wilson, Rev. Mod. Phys. 47, 773 (1975).
21. R. Bulla, T. A. Costi,  and T. Pruschke, Rev. Mod. Phys. 80, 395 (2008).
22. N. M. R. Peres, S.-W. Tsai, J. E. Santos,  and R. M. Ribeiro, Phys. Rev. B 79, 155442 (2009).
23. J. Mao, Y. Jiang, P.-W. Lo, D. May, G. Li, G.-Y. Guo, F. Anders, T. Taniguchi, K. Watanabe,  and E. Y. Andrei, arXiv:1711.06942  (2017).
24. V. G. Miranda, L. G. G. V. Dias da Silva,  and C. H. Lewenkopf, Phys. Rev. B 90, 201101 (2014).
25. D. A. Ruiz-Tijerina and L. G. G. V. Dias da Silva, Phys. Rev. B 95, 115408 (2017).
26. V. M. Pereira, J. M. B. Lopes dos Santos,  and A. H. Castro Neto, Phys. Rev. B 77, 115109 (2008).
27. A. M. Oleś, Phys. Rev. B 28, 327 (1983).
28. Reich, S. and Maultzsch, J. and Thomsen, C. and Ordejón, P., Phys. Rev. B 66, 035412 (2002).
29. B. R. K. Nanda and S. Satpathy, Phys. Rev. B 80, 165430 (2009).
30. O. V. Yazyev and L. Helm, Phys. Rev. B 75, 125408 (2007).
31. H. Padmanabhan and B. R. K. Nanda, Phys. Rev. B 93, 165403 (2016).
32. V. G. Miranda, L. G. G. V. Dias da Silva,  and C. H. Lewenkopf, Phys. Rev. B 94, 075114 (2016).
33. T. Kanao, H. Matsuura,  and M. Ogata, Journal of the Physical Society of Japan 81, 063709 (2012).
34. A. C. Hewson, The Kondo Problem to Heavy Fermions, Cambridge Studies in Magnetism (Cambridge University Press, 1993).
35. K. Chen and C. Jayaprakash, Journal of Physics: Condensed Matter 7, L491 (1995).
36. K. Ingersent, Phys. Rev. B 54, 11936 (1996).
37. R. Bulla, T. Pruschke,  and A. C. Hewson, Journal of Physics: Condensed Matter 9, 10463 (1997).
38. H. R. Krishna-murthy, J. W. Wilkins,  and K. G. Wilson, Phys. Rev. B 21, 1003 (1980a).
39. H. R. Krishna-murthy, J. W. Wilkins,  and K. G. Wilson, Phys. Rev. B 21, 1044 (1980b).
40. F. B. Anders and A. Schiller, Phys. Rev. Lett. 95, 196801 (2005).
41. F. B. Anders and A. Schiller, Phys. Rev. B 74, 245113 (2006).
42. A. Weichselbaum and J. von Delft, Phys. Rev. Lett. 99, 076402 (2007).
43. R. Peters, T. Pruschke,  and F. B. Anders, Phys. Rev. B 74, 245114 (2006).
44. D. L. Cox and A. Zawadowski, Advances in Physics 47, 599 (1998), for a review on the multi-channel models.
45. N. Roch, S. Florens, T. A. Costi, W. Wernsdorfer,  and F. Balestro, Phys. Rev. Lett. 103, 197202 (2009).
46. D. Goldhaber-Gordon, J. Göres, M. A. Kastner, H. Shtrikman, D. Mahalu,  and U. Meirav, Phys. Rev. Lett. 81, 5225 (1998).
47. P.-W. Lo, G.-Y. Guo,  and F. B. Anders, Phys. Rev. B 89, 195424 (2014).
48. T. Esat, T. Deilmann, B. Lechtenberg, C. Wagner, P. Krüger, R. Temirov, F. B. Anders, M. Rohlfing,  and F. S. Tautz, Phys. Rev. B 91, 144415 (2015).
49. A. Schiller and S. Hershfield, Phys. Rev. B 61, 9036 (2000).
50. T. O. Wehling, E. Şaşıoğlu, C. Friedrich, A. I. Lichtenstein, M. I. Katsnelson,  and S. Blügel, Phys. Rev. Lett. 106, 236805 (2011).
51. J. A. Vergés, E. SanFabián, G. Chiappe,  and E. Louis, Phys. Rev. B 81, 085120 (2010).
52. J. R. Schrieffer and P. A. Wolff, Phys. Rev. 149, 491 (1966).
53. P. S. Cornaglia, D. R. Grempel,  and H. Ness, Phys. Rev. B 71, 075320 (2005).
54. P. S. Cornaglia, P. R. Bas, A. A. Aligia,  and C. A. Balseiro, EPL (Europhysics Letters) 93, 47005 (2011).
55. D. A. Ruiz˘Tijerina, P. S. Cornaglia, C. A. Balseiro,  and S. E. Ulloa, Phys. Rev. B 86, 035437 (2012).
56. B. A. Jones and C. M. Varma, Phys. Rev. Lett. 58, 843 (1987).
57. I. Affleck, A. W. W. Ludwig,  and B. A. Jones, Phys. Rev. B 52, 9528 (1995).