Kondo effect due to a hydrogen impurity in graphene: a multichannel Kondo problem with diverging hybridization

TL;DR

This paper investigates the Kondo effect caused by a hydrogen impurity in graphene, revealing a complex multichannel Kondo problem with diverging hybridization and analyzing its ground states using numerical renormalization group methods.

Contribution

It introduces a detailed model of hydrogen impurity-induced Kondo effect in graphene, highlighting the diverging hybridization and analyzing the resulting multi-channel Kondo physics.

Findings

Ground state is a spin-1/2 doublet with ferromagnetic coupling at weak scattering.

Moderate scattering leads to a spin singlet with antiferromagnetic coupling.

Diverging local density of states influences the Kondo behavior.

Abstract

We consider the Kondo effect arising from a hydrogen impurity in graphene. As a first approximation, the strong covalent bond to a carbon atom removes that carbon atom without breaking the $C_{3}$ rotation symmetry, and we only retain the Hubbard interaction on the three nearest neighbors of the removed carbon atom which then behave as magnetic impurities. These three impurity spins are coupled to three conduction channels with definite helicity, two of which support a diverging local density of states (LDOS) $\propto 1/ [ | \omega \ | \ln ^{2}( \Lambda /| \omega \ | \ ) \ ] $ near the Dirac point $\omega \rightarrow 0$ even though the bulk density of states vanishes linearly. We study the resulting 3-impurity multi-channel Kondo model using the numerical renormalization group method. For weak potential scattering, the ground state of the Kondo model is a particle-hole symmetric spin-$1/2$ doublet, with ferromagnetic coupling between the three impurity spins; for moderate potential scattering, the ground state becomes a particle-hole asymmetric spin singlet, with antiferromagnetic coupling between the three impurity spins. This behavior is inherited by the Anderson model containing the hydrogen impurity and all four carbon atoms in its vicinity.

Tables

Table 1. Table 1: Abbreviations used in this paper.

Abbreviation	Meaning
LDOS	local density of states
RG	renormalization group
NRG	numerical renormalization group
p-h	particle-hole
RKKY	Ruderman-Kittel-Kasuya-Yosida
ALM	p-h asymmetric local moment fixed point of the 2-channel spin-$S$ Kondo problem with a logarithmically divergent LDOS
SSC	p-h symmetric strong-coupling fixed point of the 2-channel spin-$S$ Kondo problem with a logarithmically divergent LDOS
LM	free-spin p-h symmetric local moment fixed point
free-ALM	free-spin p-h asymmetric local moment fixed point
K-S	ferromagnetic p-h symmetric Kondo fixed point
F-ALM	ferromagnetic p-h asymmetric local moment fixed point
K-I	antiferromagnetic p-h symmetric isospin Kondo fixed point
AF-ASC	antiferromagnetic p-h asymmetric strong-coupling fixed point
AF-ALM	antiferromagnetic p-h asymmetric local moment fixed point
F-ASC	ferromagnetic p-h asymmetric strong-coupling fixed point

Table 2. Table 2: Properties of various fixed points of the 3-impurity 3-channel Kondo model Eq. ( 9 ). The charge number is measured relative to half-filling; we assume a negative charge if the p-h symmetry is explicitly broken by potential scattering. Pristine graphene is chosen as the reference system for S imp subscript 𝑆 imp S_{\text{imp}} and χ imp subscript 𝜒 imp \chi_{\text{imp}} . Spin ( 1 / 2 ) 3 subscript 1 2 3 \left(1/2\right)_{3} refers to three independent spin- 1 / 2 1 2 1/2 impurities; 0 < ⟨ 𝑺 1 ⋅ 𝑺 2 ⟩ ≤ 1 / 4 0 delimited-⟨⟩ ⋅ subscript 𝑺 1 subscript 𝑺 2 1 4 0<\left<\bm{S}_{1}\cdot\bm{S}_{2}\right>\leq 1/4 in the K-S phase. The results for K-S and AF-ASC also apply to the 5-atom-cluster Anderson model (see Sec.  IV ), with the exception that − 1 / 4 ≤ ⟨ 𝑺 1 ⋅ 𝑺 2 ⟩ < 0 1 4 delimited-⟨⟩ ⋅ subscript 𝑺 1 subscript 𝑺 2 0 -1/4\leq\left<\bm{S}_{1}\cdot\bm{S}_{2}\right><0 in the AF-ASC phase of the Anderson model.

Fixed point	Stability	Non-normalizable zero modes	Spin	Charge	Helicity degeneracy	$S_{\text{imp}}$	$T{\chi }_{\text{imp}}$
free-spin symmetric local moment (LM)	unstable	$2$	${\left(\frac{1}{2}\right)}_3$	$0$	-	$7\ln 2$	$1$	$0$
free-spin asymmetric local moment (free-ALM)	unstable	$0$	${\left(\frac{1}{2}\right)}_3$	$-2$	-	$3\ln 2$	$\frac{3}{4}$	$0$
ferromagnetic symmetric Kondo (K-S)	stable	$0$	$\frac{1}{2}$	$0$	-	$\ln 2$	$\frac{1}{4}$	$\in (0,\frac{1}{4}]$
ferromagnetic asymmetric local moment (F-ALM)	stable	$0$	$\frac{3}{2}$	$-2$	-	$\ln 4$	$\frac{5}{4}$	$\frac{1}{4}$
antiferromagnetic symmetric isospin Kondo (K-I)	unstable	$0$	$0$	$\pm 1$	-	$\ln 2$	$0$	$-\frac{1}{4}$
antiferromagnetic asymmetric strong-coupling (AF-ASC)	stable	$0$	$0$	$-1$	-	$0$	$0$	$-\frac{1}{4}$
antiferromagnetic asymmetric local moment (AF-ALM)	stable	$0$	$\frac{1}{2}$	$-2$	$2_{\text{spin}}$	$\ln 4$	$\frac{1}{4}$	$-\frac{1}{4}$
ferromagnetic asymmetric strong-coupling (F-ASC)	stable	$0$	$1$	$-1$	$2_{\text{channel}}$	$\ln 6$	$\frac{2}{3}$	$\in (0,\frac{1}{4})$

Equations

$$H=H_{\text{vac}}+H_{\text{hyb}}+H_{\text{imp}}.$$

$$H_{\text{vac}}=-t{\sum }_{\vec{R}}^{\prime }\left\{\left[b^{†}\left(\vec{R}\right)+b^{†}\left(\vec{R}-{\vec{a}}_2\right)+b^{†}\left(\vec{R}-{\vec{a}}_1\right)\right]a\left(\vec{R}\right)+\text{h.c.}\right\},$$

$${\vec{a}}_1=\frac{a}{2}\left(\sqrt{3},1\right),{\vec{a}}_2=\frac{a}{2}\left(\sqrt{3},-1\right),$$

$$H_{\text{hyb}}=-\sqrt{2}t\sum _{j=1}^3{a}_j^{†}{b}_j+\text{h.c.},$$

$$b(\vec{0})=b_1,b\left(-{\vec{a}}_2\right)=b_2,b\left(-{\vec{a}}_1\right)=b_3,$$

$$a_1^{†}$$

$$a_2^{†}$$

$$a_3^{†}$$

$$H_{\text{imp}}=\sum _{j=1}^3\left({ϵ}_b{n}_{b,j}+U{n}_{b,j\uparrow }{n}_{b,j\downarrow }\right)$$

$$c_{h=0}^{†}$$

$$c_{h=1}^{†}$$

$$c_{h=\bar{1}}^{†}$$

$$H=$$

$$+J_{11}\left({\mathbf{s}}_{11}+{\mathbf{s}}_{\bar{1}\bar{1}}\right)\cdot {\mathcal{S}}_0+J_{1\bar{1}}\left({\mathbf{s}}_{1\bar{1}}\cdot {\mathcal{S}}_1+{\mathbf{s}}_{\bar{1}1}\cdot {\mathcal{S}}_{\bar{1}}\right)$$

$$\displaystyle+J_{01}\bigg{[}\left(\bm{s}_{01}+\bm{s}_{\bar{1}0}\right)\cdot\mathscr{S}_{1}+\left(\bm{s}_{10}+\bm{s}_{0\bar{1}}\right)\cdot\mathscr{S}_{\bar{1}}\bigg{]},$$

$${\mathcal{S}}_h=\sum _{j=1}^3{e}^{-ih2\pi \left(j-1\right)/3}{\mathbf{S}}_j,{\mathbf{S}}_j=\frac{1}{2}\sum _{\alpha \beta }{b}_{j\alpha }^{†}{\mathbf{\sigma }}_{\alpha \beta }{b}_{j\beta },n_{b,j}=1.$$

$${\mathbf{s}}_{h{h}^{\prime }}=\frac{1}{2}\sum _{\alpha \beta }{c}_{h\alpha }^{†}{\mathbf{\sigma }}_{\alpha \beta }{c}_{h^{\prime }\beta }$$

$$J_{h{h}^{\prime }}\approx \frac{4}{3}{t}^2\left(\frac{1}{U+{ϵ}_b}+\frac{1}{-{ϵ}_b}\right),V_h\approx t^2\left(\frac{1}{-{ϵ}_b}-\frac{1}{U+{ϵ}_b}\right)$$

$$c_0^{†}$$

$$c_1^{†}$$

$$c_{\bar{1}}^{†}$$

$${\rho }_{h=0}\left(\omega \right)=\frac{3\sqrt{3}{a}^2}{8\pi v_F^2}\left|\omega \right|\text{,}$$

$${\rho }_{h=1,\bar{1}}\left(\omega \right)\approx \frac{\pi }{\sqrt{3}\left|\omega \right|{\ln }^2\frac{{\Lambda }^2}{{\omega }^2}}\text{.}$$

$$v_h\left(D\right)\equiv {\rho }_h\left(D\right)V_h\left(D\right)\text{,}$$

$$j_{h{h}^{\prime }}\left(D\right)\equiv \sqrt{{\rho }_h\left(D\right){\rho }_{h^{\prime }}\left(D\right)}{J}_{h{h}^{\prime }}\left(D\right)\text{.}$$

$$-\frac{d{j}_{11}}{d\ln D}=\left(1-\frac{2}{\ln \frac{\Lambda }{D}}\right)j_{11}+\left(j_{11}^2+j_{1\bar{1}}^2+j_{01}^2\right)\text{,}$$

$$-\frac{d{j}_{1\bar{1}}}{d\ln D}=\left(1-\frac{2}{\ln \frac{\Lambda }{D}}\right)j_{1\bar{1}}+\left(2j_{11}{j}_{1\bar{1}}+j_{01}^2\right)\text{,}$$

$$-\frac{d{j}_{00}}{d\ln D}=-j_{00}+j_{00}^2+2j_{01}^2\text{,}$$

$$-\frac{d{j}_{01}}{d\ln D}=-\frac{1}{\ln \frac{\Lambda }{D}}{j}_{01}+j_{01}\left(j_{11}+j_{1\bar{1}}+j_{00}\right)\text{,}$$

$$-\frac{d{v}_0}{d\ln D}=-v_0\text{,}$$

$$-\frac{d{v}_1}{d\ln D}=\left(1-\frac{2}{\ln \frac{\Lambda }{D}}\right)v_1\text{.}$$

$$T_K\propto J_K/{\ln }^2\left(\Lambda /J_K\right)\text{,}$$

$$T_P\propto V_1\left(D_0\right)/{\ln }^2\left(\Lambda /V_1\left(D_0\right)\right)\text{,}$$

$$H_{\text{RKKY}}=I\sum _{i<j}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

Full text

I Introduction

Interest in graphene magnetism and its potential applications in spintronics started to grow soon after the isolation and characterization of this two-dimensional materialNovoselov et al. (2004, 2005); Castro Neto et al. (2009); Yazyev (2010); Tombros et al. (2007); Yazyev and Katsnelson (2008); Soriano et al. (2010); Ortiz et al. (2018); Han et al. (2014). One particularly intensively explored approach to making graphene magnetic is through point defectsYazyev and Helm (2007); Kumazaki and Hirashima (2007); Yazyev (2008); Singh and Kroll (2009); Haase et al. (2011); Mashkoori and Jafari (2015); Charlebois et al. (2015), such as adsorbing adatomsWu et al. (2008); Fujimoto and Saito (2011); Kim and Cho (2013); Uchoa et al. (2008); Santos et al. (2010a, b); Wu et al. (2010); Pike and Stroud (2014) and vacanciesPalacios et al. (2008); Dharma-wardana and Zgierski (2008); Wang and Pantelides (2012); Chen et al. (2014); Lee et al. (2014); Li et al. (2014); Zhang et al. (2016); Valencia and Caldas (2017); Ronchi et al. (2017); Hadipour (2019). While transition metal adatoms with $d$ or $f$ electrons constitute an obvious option, rather amazingly, it has been shown both theoretically and experimentally that hydrogen impurities in graphene are also capable of inducing local magnetic momentsBoukhvalov et al. (2008); Ferro et al. (2008); McCreary et al. (2012); Sofo et al. (2012); Soriano et al. (2015); González-Herrero et al. (2016); Şaşıoğlu et al. (2017); García-Martínez et al. (2017), of the order of one Bohr magneton per defect. An intuitive explanation follows from the strong coupling between the hydrogen impurity and the carbon atom directly below itRobinson et al. (2008); Wehling et al. (2010). In the limit of this coupling going to infinity, the large energy cost of transferring an electron from or to the hydrogen-carbon pair effectively removes the $p_{z}$ orbital of the carbon atom from the graphene sheet. If we approximate graphene as a Hubbard model defined on the bipartite honeycomb lattice, then Lieb’s theorem predicts that the total spin of the ground state should be $1/2$ after the removal of the $p_{z}$ orbitalLieb (1989).

In metals with dilute magnetic impurities, conduction electrons screen the impurity magnetic moments at low temperatures, forming many-body singlets in the famous Kondo effectKondo (1964); Hewson (1997). Now understood in great detail, the Kondo effect is frequently employed in various unconventional materials in order to locally probe the bulk properties of the conducting host. From a theoretical perspective, graphene is predicted to support many exotic variants of the Kondo effect[][; andreferencestherein.]RepProgPhys.76.032501, thanks to the Dirac cones in its electronic structure and a diversity of possible impurity locations in the unit cellDell’Anna (2010); Zhu et al. (2010).

On the experimental side, it has been found early on that irradiation induced carbon vacancies in single-layer graphene produce a resistivity minimum versus temperatureChen et al. (2011), which is consistent with a high Kondo temperature ($\sim 70$K) if attributed to Kondo effect. However, subsequent magnetization measurements in irradiated thick graphite laminates suggest paramagnetism down to the lowest accessed temperaturesNair et al. (2012). While alternative explanations of the resistivity minimum based on weak localization or electron-electron interactions have been proposedJobst et al. (2013); Takehana et al. (2017); Araujo et al. (2018), the apparent contradiction between resistivity and magnetization measurements is eventually resolved by scanning tunneling spectroscopyJiang et al. (2018); May et al. (2018). The Kondo screening of vacancy-induced moments generally takes place when the graphene layer is not locally perfectly flat, and the Kondo temperature depends sensitively on the local curvature. In corrugated graphene samples, vacancies come with different local curvatures and are subjected to varying degrees of screening, but the distribution of Kondo temperatures is not fully captured by either resistivity or magnetization measurements, with the former probing screened moments and the latter probing unscreened ones.

Unlike a carbon vacancy which is subject to the Jahn-Teller distortion, a hydrogen impurity preserves the $C_{3}$ rotational symmetry of the graphene lattice around the carbon atom directly below it (henceforth referred to as the “central site”). The induced magnetic moment predominantly resides on the carbon sublattice where the impurity is not adsorbed. Despite extending many lattice constants, the magnetization is the strongest on the three nearest neighbors of the central siteFerro et al. (2008); Haase et al. (2011); González-Herrero et al. (2016); Sofo et al. (2012); García-Martínez et al. (2017). This has motivated Ref. Sofo et al., 2012 to examine, among other models, a reduced Hamiltonian where the Hubbard interactions are taken into account only on the 5-atom cluster including the hydrogen impurity, the central site and its three nearest neighbors. To study the Kondo effect, Ref. Sofo et al., 2012 replaces the rest of the system by a non-interacting bath of Dirac electrons; the cluster hybridizes with the bath via a local density of states (LDOS) that vanishes linearly as a function of energy near the Dirac points, as is characteristic of the density of states in bulk grapheneCastro Neto et al. (2009).

Upon closer inspection, however, one realizes that the system with the 5-atom cluster removed cannot be equivalent to bulk graphene, but rather comes with a 4-site vacancy: one site is removed from one of the sublattices and three sites are removed from the other. In the nearest-neighbor tight-binding model, it is known that a sublattice site number imbalance produces the same number of zero energy states (or “zero modes”) dwelling exclusively on the sublattice with more sitesPereira et al. (2008). In the simpler case of a single-site vacancy, the single zero mode produced by the vacancy cannot be normalized, because its wave function decays as $1/r$ away from the vacancyPereira et al. (2006); Peres et al. (2009); Nanda et al. (2012); Ducastelle (2013). This is intimately related to a strongly enhanced LDOS around the vacancyPereira et al. (2008), which diverges as $1/\left[\left|\omega\right|\ln^{2}\left(\Lambda/\left|\omega\right|\right)\right]$ ($\Lambda$ is a high-energy cutoff) near the Dirac point $\omega\rightarrow 0$. (We refer to this as a “logarithmic divergence” in the following.) When an impurity magnetic moment located at the vacancy is coupled to the rest of the graphene sheet, the divergent LDOS has a profound impact on the ensuing single-channel Kondo effect at half fillingMitchell and Fritz (2013); Cazalilla et al. (2012); Shirakawa and Yunoki (2014); Kanao et al. (2015); Kogan et al. (2017): in stark contrast to the linear-LDOS caseWithoff and Fradkin (1990); Cassanello and Fradkin (1996); Gonzalez-Buxton and Ingersent (1998); Fritz and Vojta (2004); Kogan (2018), both potential scattering and Kondo scattering perturbations become strongly relevant in the renormalization group (RG) sense, leading to a high Kondo temperature, and the low-energy behavior of the system is always controlled by a strong-coupling fixed point. Similar impurity-related LDOS enhancement mechanisms and their effects on Kondo screening have been discussed in the context of $d$-wave superconductorsBalatsky et al. (1995); Zhu and Ting (2000); *PhysRevB.64.060501.

In this paper, we apply the above considerations to the hydrogen impurity problem under the 5-atom cluster approximation. The non-interacting bath with the 4-site vacancy allows $3-1=2$ zero modes, both of which are non-normalizable. Correspondingly, we show that two conduction channels with a diverging hybridization appear in the Kondo problem. There is a third conduction channel with a vanishing LDOS at low energies, but the perturbations associated with it are strongly irrelevant; its importance is therefore diminished by the other two channels.

As a first approximation, we consider the limit of infinite coupling between the hydrogen impurity and the central site. The central site is essentially eliminated from the low-energy theory in this limit. Consequently, only the three nearest neighbors of the central site are left to host the magnetic moment, mapping to a 3-impurity Kondo problem where the impurities are symmetric under $Z_{3}$ permutations. In a metallic host, this problem is known to yield a rich phase diagramPaul and Ingersent (1996); Ingersent et al. (2005), and many insights are carried over to our case of a diverging hybridization.

To tackle the 3-impurity problem, we first construct an auxiliary problem: a single impurity spin of arbitrary size $S$ coupled to two conduction channels through the same diverging hybridization. This is analyzed with the aid of the numerical renormalization group (NRG) algorithmWilson (1975); Bulla et al. (2008); Z̆itko (2019). Systematic studies have been performed on the pseudogap case where the hybridization vanishes at zero energySchneider et al. (2011). However, as in the single-channel case, the diverging hybridization makes a qualitative difference. We find that all low-energy fixed points are strong-coupling and Fermi-liquid like; in particular, in the presence of particle-hole (p-h) symmetry, the low-energy fixed point involves the impurity spin screened by conduction electrons from both channels, forming a residual spin of size $\left|S-1\right|$ for any $S$, including $S=1/2$ together with two phase-shifted conduction channels. This is very different from the 2-channel spin-$1/2$ Kondo problem with a constant hybridization, whose ground state is a non-Fermi liquidAffleck and Ludwig (1993, 1991); Andrei and Destri (1984); Tsvelick and Wiegmann (1985).

Returning to the 3-impurity Kondo model, we map out the phase diagram by NRG and study the thermodynamics and the impurity spin correlations in each phase. The most important coupling constants are the relevant ones associated with potential and Kondo scattering in the two conduction channels with a diverging hybridization. For weak and intermediate potential scattering, we find two stable low-energy fixed points: a p-h symmetric spin-$1/2$ fixed point with the Kondo effect taking place in the spin sector which we label as K-S, and a p-h asymmetric spin-singlet fixed point which we label as AF-ASC. The latter is connected to an unstable fixed point of the Kondo effect taking place in the isospin sector. K-S is favored by weak potential scattering at the impurities and ferromagnetic Ruderman-Kittel-Kasuya-Yosida (RKKY) coupling between the impurities, and exhibits ferromagnetic impurity spin correlations, whereas AF-ASC exhibits antiferromagnetic impurity spin correlations; the transition between the K-S phase and the AF-ASC phase is shown to be a simple level crossing. On the other hand, we find that very strong potential scattering can overwhelm the Kondo scattering and suppress the Kondo effect, in which case the impurity spins couple to form a magnetic moment decoupled from the conduction electrons. Finally, the divergence of the hybridization is inevitably cut off at low energies in more realistic models of graphenePereira et al. (2008), and we examine the consequences of such a cutoff on a phenomenological level.

We then proceed to analyze with NRG the Anderson model with Hubbard interactions on the 5-atom cluster. The system flows to the p-h symmetric spin-$1/2$ K-S fixed point when p-h symmetry breaking perturbation is weak, and to the p-h asymmetric spin-singlet AF-ASC fixed point otherwise; the two phases are again respectively characterized by ferromagnetic and antiferromagnetic impurity spin correlations. We conclude that the Kondo effect occurs at the hydrogen impurity both in the 3-impurity Kondo model and the 5-atom cluster Anderson model.

The rest of this paper is organized as follows. In Sec. II, we introduce the 3-impurity Kondo model of a hydrogen impurity in the infinite hydrogen-carbon coupling limit, highlighting the diverging hybridization between the magnetic impurities and the two conduction channels. Sec. III is devoted to the 3-impurity Kondo model: we give the scaling behavior of various perturbations and the RKKY interactions in the weak-coupling limit in Sec. III.1, then discuss the NRG results on the 2-channel spin-$S$ Kondo model with a diverging hybridization in Sec. III.2. The numerical results on the 3-impurity Kondo model with a diverging hybridization are analyzed in depth in Sec. III.3. Sec. III.4 closes our discussion of the 3-impurity Kondo model by demonstrating the effects of a low-energy cutoff on the divergent LDOS. In Sec. IV we interpret our numerical results for the 5-atom cluster Anderson model with a diverging hybridization. Sec. V concludes the paper and discusses some open problems. Appendix A contains a derivation of the divergent LDOS from the 4-site vacancy lattice model, and the corresponding zero mode solutions of this model are discussed in Appendix B. Finally, in Appendix C, we calculate the RKKY interaction in the 3-impurity Kondo model to the second order in Kondo couplings. The abbreviations used in this paper are summarized in Table 1.

II Model

We consider the nearest-neighbor tight-binding model of a graphene layer defined on a bipartite honeycomb lattice with two sublattices A and B. The “central site”, to which the hydrogen impurity is coupled, is assumed to be on the A sublattice. Throughout the paper, we follow Ref. Sofo et al., 2012 and only retain the local Hubbard interactions on the 5-atom cluster, composed of the hydrogen impurity, the central A site and its three nearest-neighbor B sites. To highlight the Kondo physics, in this section we furthermore approximate the hydrogen-carbon coupling as infinity, which leaves both the hydrogen and the central site decoupled from the rest of the system. In the limit of strong Hubbard interactions, we effectively localize the electrons on the three nearest-neighboring B sites. The effect of the hopping from the 3 B sites to the 6 outer A sites can be considered to leading order in perturbation theory. The associated Anderson model corresponds to 3 interacting impurities hybridizing with a non-interacting bath through the outer A electrons:

[TABLE]

The non-interacting bath part of the effective 3-impurity Hamiltonian describes a graphene sheet with a 4-site vacancy

[TABLE]

where the summation excludes the central A site and its three nearest neighbor B sites, as shown in Fig. 1. Here, $\vec{R}=n\vec{a}_{1}+m\vec{a}_{2}$ are the Bravais lattice vectors with

[TABLE]

where $a$ is the lattice constant, and the B site labeled by $\vec{R}$ is displaced by $(a/\sqrt{3},0)$ from the corresponding A site. The four sites removed from $H_{\text{vac}}$ are then labeled by $a(\vec{0})$, $b(\vec{0})$, $b\left(-\vec{a}_{1}\right)$ and $b\left(-\vec{a}_{2}\right)$.

The impurity-bath hybridization, given by

[TABLE]

is invariant under $C_{3}$ rotations. Here we have relabeled the 3 impurity B sites as

[TABLE]

and defined the symmetric linear combinations of pairs of the neighboring $a^{\dagger}(\vec{R})$ electrons as

[TABLE]

Finally, the Hamiltonian of the 3 nearest-neighbor B sites consists of a local Hubbard interaction term and a local on-site potential term:

[TABLE]

where $n_{b,j\alpha}=b_{j\alpha}^{{\dagger}}b_{j\alpha}$, and $n_{b,j}=n_{b,j\uparrow}+n_{b,j\downarrow}$ is the number operator for B electrons at site $j$. This model is p-h symmetric when $\epsilon_{b}=-U/2$ and the chemical potential $\mu=0$.

The lattice with a 4-site vacancy inherits the $C_{3}$ symmetry of the pristine lattice. Hence, we can construct helicity eigenstates $h=0,1,\bar{1}$ from the $a_{1,2,3}^{\dagger}$ states:

[TABLE]

A counter-clockwise $2\pi/3$ rotation about the central A site acts as a permutation of the three $a_{1,2,3}^{\dagger}$ states. Then, under this rotation we have $c_{h}^{\dagger}\rightarrow e^{i2\pi h/3}c_{h}^{\dagger}$.

In the limit $U\sim\left|\epsilon_{b}\right|\gg t$, by applying a Schrieffer-Wolff projection as in Ref. Paul and Ingersent, 1996, we obtain an effective Kondo model which includes potential scattering and Kondo interactions:

[TABLE]

where the local moment operators of definite helicities arePaul and Ingersent (1996); Ingersent et al. (2005)

[TABLE]

We note that $\mathscr{S}_{h=0}$ is simply the total impurity spin operator, and the total impurity spin quantum number can be $S=3/2$ (one quartet) or $S=1/2$ (two doublets); the two $S=1/2$ doublets can be distinguished by their behavior under the $Z_{3}$ permutation of impurity spins. We have also defined the particle number operators for $c_{h}$, $n_{h}=\sum_{\alpha}c_{h\alpha}^{\dagger}c_{h\alpha}$, and the spin operators

[TABLE]

which can involve conduction electrons of different helicities. For our particular microscopic model, the unrenormalized couplings are

[TABLE]

Generally $\left|V_{h}\right|/J_{hh^{\prime}}\leq 3/4$. Note that $V_{h}$ will be generated by breaking the p-h symmetry. This can occur not only when $\epsilon_{b}\neq-U/2$, but also when we move away from the charge neutrality point or take into consideration second neighbor hoppingMalecki and Affleck (2010).

The scaling dimensions of the various couplings in $H$ are determined by the LDOS of the $c_{h}$ conduction channels for the 4-site-vacancy graphene. A detailed solution of this non-interacting problem in Appendix A gives the following leading contributions to the $c_{h}$ channels:

[TABLE]

where $\Lambda\sim t$ is an ultraviolet energy cutoff, $v_{F}=\sqrt{3}ta/2$ is the Fermi velocity, and $\tilde{\phi}_{\vec{K}/\vec{K}^{\prime},m,k}^{\dagger}$ creates an electron in the eigenstate of $H_{\text{vac}}$ in valley $\vec{K}$ or $\vec{K}^{\prime}=\left(\sqrt{3},\pm 1\right)\left(2\pi/3a\right)$, with angular momentum $m$ and momentum amplitude $k$. (The low energy Dirac theory has full rotational symmetry so that eigenstates can be labeled by $m$, the 2-dimensional angular momentum quantum numberCassanello and Fradkin (1996).) Note that the low-energy spectrum is determined from $\epsilon_{\vec{K}+\vec{k}}\approx v_{F}k$. Due to the additional factor of $a$ in Eq. (13a), $c_{0}$ should become less and less important compared to $c_{1}$ and $c_{\bar{1}}$ at low energy scales, as will be confirmed in Sec. III. From this, is it straightforward to determine the leading contributions to the LDOS for $c_{h}$ in the low-energy limit:

[TABLE]

[TABLE]

While the helicity-0 channel has a behavior similar to pristine graphene, helicities $1$ and $\bar{1}$ show a logarithmic divergence in their LDOS. We attribute such a divergence to the presence of two non-normalizable zero modes in the 4-site-vacancy graphene, whose wave functions behave as $1/r$ when the distance to the vacancy $r$ is largePeres et al. (2009); Pereira et al. (2006); Ducastelle (2013); see Appendix B for details. Because $N_{A}=1$ A site and $N_{B}=3$ B sites are removed, there are $\left|N_{A}-N_{B}\right|=2$ zero modes (for each spin) living on the A sublattice; due to the $C_{3}$ symmetry in this case, these two zero modes can be chosen as helicity eigenstatesPereira et al. (2008). Even though these zero modes are not true eigenstates in an infinite system, they hybridize strongly with the low-energy itinerant states in pristine graphene, forming low-energy scattering states which are true eigenstates of $H_{\text{vac}}$ and contribute to the divergent LDOSDucastelle (2013).

III Kondo model

In this section, we establish the phase diagram of the 3-impurity 3-channel Kondo model Eq. (9), using a combination of analytical arguments and NRG.

III.1 Scaling and RKKY interactions at weak coupling

It is instructive to begin by analyzing the weak-coupling fixed point. To find the scaling behavior of various coupling constants, we first define dimensionless couplings at the running energy cutoff $D$:

[TABLE]

and

[TABLE]

The second-order weak-coupling RG equations then read

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

From these equations we can also see that the relation $J_{11}=J_{1\bar{1}}$, if true for the bare couplings, is preserved along the RG flow.

Due to the singular LDOS for helicities $1$ and $\bar{1}$, $J_{11}$ and $J_{1\bar{1}}$ are relevant at low energies ($D\ll\Lambda$), in analogy to the single-channel problem with a divergent LDOS discussed in Ref. Mitchell and Fritz, 2013. This leads to a greatly enhanced Kondo temperature for the corresponding Kondo couplings:

[TABLE]

where $J_{K}$ is either $J_{11}\left(D_{0}\right)$ or $J_{1\bar{1}}\left(D_{0}\right)$, $D_{0}$ being the initial semi-bandwidth of the Kondo model. The potential scattering term $V_{1}$ is likewise relevant, and has its own characteristic energy scale

[TABLE]

at which it flows to strong-coupling. On the other hand, $J_{01}$ is weakly irrelevant, becoming almost marginal only at very low energies, even though it generates relevant couplings $J_{11}$ and $J_{1\bar{1}}$ at the second order. Finally, the linear LDOS of the helicity-0 channel renders $V_{0}$ and $J_{00}$ strongly irrelevant.

It is also possible to consider the RKKY interactions between magnetic impurities at weak coupling mediated by conduction electronsJones and Varma (1987); *PhysRevLett.61.125; *PhysRevB.40.324; Paul and Ingersent (1996); this gives us some intuition on possible magnetic orders of the impurities. It should be clarified that these interactions are only introduced to help us understand the Kondo model; they are not part of the NRG input (as the Kondo couplings are), and we make no a priori assumptions about the associated RKKY energy scale in our NRG calculations. As in Refs. Paul and Ingersent, 1996; Ingersent et al., 2005, the RKKY interactions are of the form

[TABLE]

where we labeled a generic RKKY interaction by $I$ in order to avoid confusion with the valley-momenta $\vec{K},\vec{K}^{\prime}$. This expression can be re-cast as

[TABLE]

$H_{\text{RKKY}}$ takes the value $3I/4$ in the $S=3/2$ multiplet state, and $-3I/8$ otherwise. Hence, strong antiferromagnetic RKKY interactions ($I>0$) project onto the $S=1/2$ manifold, while strong ferromagnetic RKKY ($I<0$) prefers the $S=3/2$ configuration.

The RKKY coupling strength $I$ has been evaluated in bulk grapheneSaremi (2007); Bunder and Lin (2009); Black-Schaffer (2010); Sherafati and Satpathy (2011); Kogan (2011), under the assumption that each magnetic impurity interacts with one carbon atom but does not disrupt the graphene lattice (e.g. by introducing vacancies). In that case, it has been shown to be ferromagnetic between impurities on the same sublattice, and antiferromagnetic between impurities on different sublattices. Nevertheless, antiferromagnetic RKKY interactions have also been reported between same-sublattice magnetic impurities when other non-magnetic impurities are presentLee et al. (2012), or when a large on-site potential energy is associated with the magnetic impuritiesAgarwal and Mishchenko (2017).

In Appendix C, we analyze the RKKY interaction between the three effective magnetic impurities $b_{1,2,3}$ generated by Kondo couplings to $O\left(J_{hh^{\prime}}^{2}\right)$, carefully including the effects of the vacancy. At low temperatures, we find $I\propto J_{1\bar{1}}^{2}-2J_{11}^{2}$; thus $I$ is ferromagnetic in a model with $J_{11}$ only, and antiferromagnetic in a model with $J_{1\bar{1}}$ only. In the Kondo model obtained through the Schrieffer-Wolff transformation, where $J_{11}=J_{1\bar{1}}$, $I$ is expected to be ferromagnetic. The RG flow of the RKKY interaction is controlled by

[TABLE]

In other words, the RKKY interaction does not change sign along the RG flow near the weak coupling fixed point.

III.2 Auxiliary model: 2-channel spin-$S$ Kondo model with a logarithmically divergent LDOS

The weak-coupling analysis in Sec. III.1 shows that all coupling constants associated with the helicity-0 conduction channel are irrelevant, and the collective state of the impurity spins can be either an $S=3/2$, helicity-0 multiplet or an $S=1/2$, helicity-$\pm 1$ multiplet. We are therefore motivated to study the 2-channel spin-$S$ Kondo model where both conduction channels $1$ and $\bar{1}$ are characterized by the logarithmically divergent LDOS Eq. (15):

[TABLE]

with an antiferromagnetic Kondo coupling $J>0$. Although this model is interesting in its own right, it is quite different from Eq. (9) even after ignoring the helicity-0 conduction channel [see Eq. (25)]. In Eq. (9) we have spin operators of conduction electrons that are diagonal ($\bm{s}_{11}$, $\bm{s}_{\bar{1}\bar{1}}$) and off-diagonal ($\bm{s}_{1\bar{1}}$, $\bm{s}_{\bar{1}1}$) in the helicity basis, while in Eq. (24) we only have diagonal operators in the helicity basis; also, the impurity in Eq. (9) comprises three $S=1/2$ spins appearing in different total spin sectors and different helicity combinations, whereas the impurity in Eq. (24) is a single spin of rigid size $S$. Nevertheless, as we will show in Sec. III.3, the auxiliary model Eq. (24) provides helpful intuitions for understanding the limiting cases of the 3-impurity model, where the impurity can be effectively viewed as a rigid spin. The auxiliary model also serves as an excellent benchmark for the impurity contributions to thermodynamic quantities, which we will discuss later in this section.

Let us briefly review the case where both channels have a constant LDOS. In this case, at the weak-coupling fixed point, $V$ is exactly marginal and $J$ is marginally relevant due to the dynamics of the impurity spin. It is well known that the low-energy behavior of this problem depends on the size of the impurity spinNozières, Ph. and Blandin, A. (1980); Affleck (1995). When $S>1$ (underscreened) or $S=1$ (exactly screened), the low-energy fixed point is simply the strong-coupling one. This involves a decoupled residual impurity spin of size $S-1$ (when $S>1$) and a Fermi liquid theory for the conduction electrons, with a phase shift $\pi/2+\delta$ ($\delta$ being an odd function of $V$) for each conduction channel taking part in the screening of the impurity. On the other hand, in the overscreened case $S=1/2$, the strong-coupling fixed point is no longer stable. If we consider a naive $J\rightarrow+\infty$ theory on a nearest-neighbor tight-binding lattice with hopping amplitude $t\ll J$, then the electrons $c_{1}$, $c_{\bar{1}}$ are strongly bound to the impurity spin, forming an effective spin-1/2. The effective spin is in turn coupled to the lattice with $c_{1}$ and $c_{\bar{1}}$ removed (which we name as the “strong-coupling lattice”). This $O\left(t^{2}/J\right)$ coupling constant is antiferromagnetic, which we already know to be marginally relevant. In contrast, the underscreened strong-coupling fixed point is stable because the large-$J$ effective model is a Kondo model with a ferromagnetic Kondo coupling, which is marginally irrelevant (and turns the system into a singular Fermi liquidMehta et al. (2005)). The stable low-energy fixed point in the overscreened case is therefore an intermediate-coupling, non-Fermi-liquid one. Its properties can be obtained by a range of theoretical methods, including boundary conformal field theoryAffleck and Ludwig (1993), Bethe ansatzAndrei and Destri (1984); Tsvelick and Wiegmann (1985) and Abelian bosonizationEmery and Kivelson (1992); in particular an impurity entropy of $\frac{1}{2}\ln 2$ and an impurity magnetic susceptibility that depends logarithmically on temperature have been predicted and numerically verified.

We return to the 2-channel Kondo model Eq. (24) with a divergent LDOS Eq. (15). As discussed in Sec. III.1, at the weak-coupling, p-h symmetric local moment fixed point, both $J$ and $V$ are relevant on account of the energy dependence of the LDOS. It is therefore natural to investigate the corresponding strong-coupling theories on a lattice with hopping amplitude $t\ll J$. The $J\rightarrow+\infty$ lattice theory as before gives rise to an effective impurity spin of size $\left|S-1\right|$ formed by the impurity spin and the electrons $c_{1}$ and $c_{\bar{1}}$, which is coupled to the remaining “strong-coupling lattice” with an either ferromagnetic (when $S>1$) or antiferromagnetic (when $S=1/2$) coupling constant of $O\left(t^{2}/J\right)$. In the $V\rightarrow\pm\infty$ lattice theory, on the other hand, no Kondo physics takes place; the $c_{1}$ and $c_{\bar{1}}$ states become nevertheless inaccessible to other conduction electrons, being either empty or occupied depending on the sign of $V$, and the original spin $S$ is coupled to the remaining “strong-coupling lattice”.

At this point we should accentuate the crucial difference between the logarithmically divergent LDOS and the constant LDOS. In the constant LDOS case, removing $c_{1}$ and $c_{\bar{1}}$ from the original lattice does not change the LDOS at the sites which the effective impurity can couple to, because the original lattice and the “strong-coupling lattice” are both composed of one-dimensional chains. In graphene with a 4-site vacancy, however, projecting out $c_{1}$ and $c_{\bar{1}}$ will remove the logarithmic divergence in the LDOS, and the leading contribution to the LDOS becomes linear near the Dirac point as in the bulk. This can be shown explicitly using the method outlined in Appendix A. On a more intuitive level, we can also explain the disappearance of the logarithmic divergence through the removal of the two non-normalizable zero modes. In Appendix B, we see that the wave functions associated with $c_{1}$ and $c_{\bar{1}}$ cannot be simultaneously zero for a nontrivial zero-energy solution of the lattice Schrödinger equation; in other words, if we project out $c_{1}$ and $c_{\bar{1}}$ from the lattice, i.e. demand the wave functions associated with $c_{1}$ and $c_{\bar{1}}$ should vanish, then no zero mode exists.

Because of the linear LDOS in the “strong-coupling lattice”, any perturbation at the naive strong-coupling fixed point– Kondo coupling or potential scattering– is strongly irrelevantFritz and Vojta (2004). It is therefore reasonable to conclude that the large-$J$ and large-$V$ strong-coupling fixed points are stable in the 2-channel Kondo model Eq. (24) with a divergent LDOS Eq. (15), irrespective of the size of the impurity spin $S$.

This picture is verified by NRG calculations, which we perform with the “NRG Ljubljana” codeZ̆itko (2019). The schematic phase diagram is shown in Fig. 2. As conjectured, the two regimes are controlled by the large-$J$ p-h symmetric strong-coupling (SSC) fixed point and the large-$V$ p-h asymmetric local moment (ALM) fixed point, with residual spin sizes $\left|S-1\right|$ and $S$ respectively. At both fixed points it is possible to construct the entire finite-size spectrum from single-particle excitations, in contrast to the non-Fermi-liquid overscreened Kondo fixed point in the case of constant LDOS. As in the single-channel spin-1/2 case discussed in Ref. Mitchell and Fritz, 2013, the two phases are separated by a second-order transition.

To shed further light on the nature of the low-energy fixed points, let us examine their thermodynamic properties. The impurity contribution to any quantity $\Omega$ in a quantum impurity system is defined as $\Omega_{\text{imp}}=\Omega-\Omega_{0}$, the difference between this quantity evaluated in the entire impurity system and in the reference “clean” system without the impurity. In our case, there are two possibilities for the reference system: the pristine graphene, and the non-interacting bath of a 4-site-vacancy graphene lattice with a logarithmically divergent LDOS. Following Ref. Mitchell and Fritz, 2013, we choose the 4-site-vacancy graphene as the reference system while presenting our numerical results, but impurity quantities measured with respect to pristine graphene will also be discussed because they are experimentally directly accessible (see Table 2). We focus on the impurity entropy $S_{\text{imp}}=-\partial\mathscr{F_{\text{imp}}}/\partial T$ ($\mathscr{F}$ being the free energy) and the impurity magnetic susceptibility $\chi_{\text{imp}}=\left\langle S_{z}^{2}\right\rangle_{\text{imp}}/T$.

It is useful to first discuss the effect of non-normalizable zero modes. In short, with respect to pristine graphene, each non-normalizable zero mode contributes $\ln 4$ to the zero-temperature impurity entropy, and $1/\left(8T\right)$ to the magnetic susceptibilityMitchell and Fritz (2013). We can directly derive these results by viewing pristine graphene as a (non-interacting) resonant-level model with a logarithmically divergent LDOS, as has been done in Appendix B of Ref. Mitchell and Fritz, 2013. Alternatively, when we calculate the impurity-induced density of states (which appears in the Friedel sum rule)Hewson (1997) in the single-vacancy graphene lattice, we find a $\delta$-function peak at zero energyPeres et al. (2009); Pereira et al. (2006); Ducastelle (2013); this $\delta$-function is what contributes to the zero-temperature impurity entropy and the impurity susceptibility, as if it were a real spin-degenerate single-particle eigenstate of the system. When the non-normalizable zero mode is removed, for instance at strong coupling, the impurity contributions vanish correspondingly. This applies equally to the 4-site-vacancy graphene lattice, except there are now two channels with a logarithmically divergent LDOS, so that the strength of the $\delta$-function peak is doubled together with the impurity entropy and impurity magnetic susceptibility.

We can now determine the limiting behavior of the thermodynamic properties at the fixed points of the 2-channel spin-$S$ Kondo model Eq. (24) with a divergent LDOS Eq. (15). At the ALM fixed point, the non-normalizable zero modes are removed by potential scattering; the local moment then results in a zero-temperature impurity entropy of $S_{\text{imp}}\left(T=0\right)=\ln\left(2S+1\right)$ and an impurity susceptibility of $\chi_{\text{imp}}=S\left(S+1\right)/\left(3T\right)$ relative to pristine graphene. At the SSC fixed point, the zero modes are removed by strong Kondo screening, so the residual spin $\left|S-1\right|$ yields $S_{\text{imp}}\left(T=0\right)=\ln\left(2\left|S-1\right|+1\right)$ and $T\chi_{\text{imp}}=\left|S-1\right|\left(\left|S-1\right|+1\right)/3$. Finally, at the p-h symmetric local moment fixed point, the zero modes together with the local moment give $S_{\text{imp}}\left(T=0\right)=\ln\left(2S+1\right)+2\ln 4$ and $T\chi_{\text{imp}}=1/4+S\left(S+1\right)/3$.

The high- to low-temperature crossover of $S_{\text{imp}}$ and $T\chi_{\text{imp}}$ is plotted in Fig. 3, with the reference system chosen as the 4-site-vacancy graphene (where $S_{\text{imp}}\left(T=0\right)=2\ln 4$ and $T\chi_{\text{imp}}=2\cdot 1/8=1/4$). In the $S=1/2$ case, the crossover between the unscreened $S=1/2$ and the overscreened $S=1/2$ is clearly visible from the non-monotonicity of $T\chi_{\text{imp}}$. At low temperatures, $S_{\text{imp}}$ and $T\chi_{\text{imp}}$ have logarithmic corrections in the form of $1/\ln\left(\Lambda/T\right)$ near all fixed pointsMitchell and Fritz (2013). However, unlike in a singular Fermi liquid, this logarithmic behavior does not originate from a marginally irrelevant operator, because the allowed operators at any of these fixed points are strongly irrelevant irrespective of the spin size. We can verify that the logarithmic corrections vanish when the reference system is pristine graphene (in the form of a non-interacting 4-site cluster coupled to the 4-site-vacancy graphene); therefore, the logarithmic behavior in Fig. 3 can be fully attributed to the non-normalizable zero modes of the reference system.

III.3 Phases of the 3-impurity model

We are in the position to present the NRG results on the 3-impurity 3-channel Kondo model Eq. (9). Again, as a first approximation, we neglect the helicity-0 conduction channel entirely on the grounds that all couplings $J_{00}$, $J_{01}$ and $V_{0}$ associated with it are irrelevant; we shall see later that this approximation is usually justified. This leaves us with a 3-impurity 2-channel Kondo model,

[TABLE]

with two relevant Kondo couplings $J_{11}$ and $J_{1\bar{1}}$ as well as a relevant potential scattering $V_{1}$. As in Refs. Paul and Ingersent, 1996; Ingersent et al., 2005, one may measure the relative strengths of $J_{11}$ and $J_{1\bar{1}}$ with the dimensionless RKKY coupling strength $\tilde{I}\equiv\left(J_{1\bar{1}}^{2}-2J_{11}^{2}\right)/\left(2J_{11}^{2}+2J_{1\bar{1}}^{2}\right)$. $J_{11}=J_{1\bar{1}}=V_{1}=0$ marks the unstable p-h symmetric local moment (“LM”) fixed point, which has three decoupled impurity spins and two non-normalizable zero modes; therefore, at the LM fixed point, $S_{\text{imp}}\left(T=0\right)=3\ln 2+2\ln 4=7\ln 2$, and $T\chi_{\text{imp}}=3\cdot 1/4+2\cdot 1/8=1$ relative to pristine graphene.

In Fig. 4 we present the NRG phase diagram of the 3-impurity 2-channel Kondo model Eq. (25). Panel (a) is the phase diagram on the $\tilde{I}$-$V_{1}$ plane when $J_{11}^{2}+J_{1\bar{1}}^{2}=\left(0.1D_{0}\right)^{2}$, and panel (b) is the phase diagram on the $J_{11}$-$V_{1}$ plane for $J_{11}=J_{1\bar{1}}$. We also show the high- to low-temperature crossover of $S_{\text{imp}}$ and $T\chi_{\text{imp}}$ in different phases in Figs. 5 and 6.

In the following we discuss the different regimes in the phase diagram. We begin by considering each of the two relevant Kondo couplings $J_{11}$ and $J_{1\bar{1}}$ separately.

(i) $J_{11}>0,J_{1\bar{1}}=0$:

The RKKY interaction is ferromagnetic ($\tilde{I}=-1$), favoring $S=3/2$; therefore it is plausible that, at low energies, the model is reduced to a 2-channel problem with a logarithmically divergent LDOS and a single $S=3/2$ impurity. This picture for $J_{1\bar{1}}=0$ is confirmed by NRG, and the coupling between the two spin sectors $S=3/2$ and $S=1/2$ turns out to be irrelevant.

Applying our results from Sec. III.2, we find two different phases: a p-h symmetric strong-coupling phase and a p-h asymmetric local moment phase. In the symmetric phase, the effective $S=3/2$ impurity is strongly screened by both $c_{1}$ and $c_{\bar{1}}$, and a residual $S=1/2$ impurity emerges; thus we name it “K-S” after the Kondo screening taking place in the spin sector. In the asymmetric phase, the effective $S=3/2$ local moment is unscreened, while both $c_{1}$ and $c_{\bar{1}}$ are blocked locally by strong potential scattering; we call this phase “F-ALM” after the ferromagnetic RKKY interaction. K-S and F-ALM are both stable fixed points as discussed in Sec. III.2; they are separated by a critical boundary marked by $\left|V_{1}\right|/J_{11}\sim O\left(1\right)$, in analogy to the single-channel caseMitchell and Fritz (2013). The low-temperature thermodynamic properties of K-S (F-ALM) is identical to those of the SSC (ALM) fixed point in the 2-channel $S=3/2$ problem with a logarithmically divergent LDOS. Therefore, at K-S, $S_{\text{imp}}\left(T=0\right)=\ln 2$ and $T\chi_{\text{imp}}=1/4$ relative to pristine graphene; at F-ALM, $S_{\text{imp}}\left(T=0\right)=\ln 4$ and $T\chi_{\text{imp}}=5/4$ relative to pristine graphene.

While our heuristic picture correctly predicts the thermodynamic quantities and the finite-size spectrum, it would only be strictly true had we introduced by hand a strong ferromagnetic RKKY interaction Eq. (21) into the Hamiltonian, assuming that the RKKY strength $I$ were much greater in magnitude than any other energy scale in the problem. In reality, there might be no clear separation between the RKKY energy scale, at which the effective $S=3/2$ impurity forms, and the Kondo temperature $T_{K}$ (or $T_{P}$), at which $J_{11}$ (or $V_{1}$) flows to strong coupling. These energy scales must be extracted numerically, e.g. from thermodynamics (Fig. 5) and impurity spin correlations (Fig. 7 which we will discuss later).

(ii) $J_{1\bar{1}}>0,J_{11}=0$:

The RKKY interaction is now antiferromagnetic ($\tilde{I}=1/2$), and favors the $S=1/2$ state for the magnetic impurities. While one might assume that the low-energy physics would be captured by the 2-channel $S=1/2$ Kondo model with a logarithmically divergent LDOS, we must also take note of the additional helicity degeneracy of the $S=1/2$ subspace, and the fact that $J_{1\bar{1}}$ is not a conventional Kondo coupling.

We can again glean some insight from the constant-LDOS version of this problemPaul and Ingersent (1996); Ingersent et al. (2005). In the constant-LDOS 3-impurity 3-channel Kondo problem, when the RKKY interactions are antiferromagnetic and $J_{1\bar{1}}$ overwhelms $J_{01}$, it has been reported that the helicity-0 channel decouples, and the low-energy effective model is a 2-channel spin-1/2 Kondo model with spin and isospin sectors interchanged. (The isospin for helicities $1$ and $\bar{1}$ is defined as usual by $\hat{I}^{z}=\frac{1}{2}\sum_{h=\pm 1,\alpha}\left(c_{h\alpha}^{{\dagger}}c_{h\alpha}-1/2\right)$ and $\hat{I}^{+}=\frac{1}{2}\sum_{h=\pm 1,\alpha\beta}\epsilon_{\alpha\beta}c_{h\alpha}^{{\dagger}}c_{\bar{h}\beta}^{{\dagger}}$, where $\epsilon_{\alpha\beta}$ is the Levi-Civita symbol.) To be more concrete, the effective model involves an isospin-1/2 impurity screened by two conduction channels, one spin-up and the other spin-down. We emphasize that such a fictitious isospin-1/2 impurity is only invoked to describe the isospin state of the conduction electrons, since the impurities themselves always possess p-h symmetry and have isospin 0 by construction. The resulting low-energy non-Fermi-liquid fixed point, dubbed “isospin Kondo” in Ref. Paul and Ingersent, 1996, is unstable against an infinitesimal p-h symmetry breaking perturbation $V_{1}$, which plays the role of a magnetic field in the isospin sector; the system flows to a Fermi-liquid state in the presence of $V_{1}$.

In our problem with a logarithmically divergent LDOS, the condition of $J_{1\bar{1}}$ dominating over $J_{01}$ is always satisfied when the bare couplings are weak, since $J_{1\bar{1}}$ is relevant while $J_{01}$ is irrelevant. We find that, as in the constant-LDOS case, the low-energy physics is governed by a 2-channel spin-1/2 Kondo model with spin and isospin sectors interchanged; however, as shown in Sec. III.2, the logarithmically divergent LDOS dictates that the low-energy fixed point for $V_{1}=0$ should be located at strong coupling rather than intermediate coupling. We label this strong-coupling fixed point as “K-I” due to the Kondo effect taking place in the isospin sector. The ground state of K-I has spin zero and isospin $1/2$, i.e. one electron is either removed from or added to half filling. This yields $S_{\text{imp}}\left(T=0\right)=\ln 2$ and $T\chi_{\text{imp}}\rightarrow 0$ at the K-I fixed point relative to pristine graphene. As an intuitive picture, in the lattice version of K-I, the impurity spins form an effective spin-$1/2$, which is in turn strongly coupled to (i.e. screened by) one of the conduction channels; the other conduction channel can be either empty or doubly occupied at the lattice site closest to the impurity.

K-I is unstable against an infinitesimal $V_{1}$, which picks out a preferred isospin state from the doublet (i.e. one electron removed from or added to half filling). We call the resulting p-h asymmetric strong-coupling fixed point “AF-ASC”. In the lattice version of AF-ASC, one conduction channel forms a spin singlet with the impurities, and the other channel takes advantage of the local potential scattering to lower the ground state energy; the remaining conduction electrons are essentially free apart from constraints imposed by the Pauli principle. With the isospin symmetry broken and the ground state a spin singlet, we simply have $S_{\text{imp}}\left(T=0\right)=0$ and $T\chi_{\text{imp}}\rightarrow 0$ at AF-ASC relative to pristine graphene. Fig. 6 shows how the system flows from K-I to AF-ASC for $V_{1}/D_{0}$ ranging from from $10^{-7}$ to $0.1$.

Increasing $V_{1}$ further, for sufficiently large $\left|V_{1}\right|/J_{11}\sim O\left(1\right)$, we eventually encounter a second transition to the large-$V_{1}$ fixed point, where both $c_{1}$ and $c_{\bar{1}}$ become blocked by strong potential scattering and the ground state electric charge differs from half-filling by two. The antiferromagnetic RKKY interactions remain in effect even though the impurity spins are already decoupled from the conduction electrons, so we are left with an $S=1/2$ local moment with an additional helicity degeneracy $h=\pm 1$. Therefore, at this fixed point which we call “AF-ALM”, $S_{\text{imp}}\left(T=0\right)=\ln 4$ and $T\chi_{\text{imp}}=1/4$ relative to pristine graphene.

In addition to cases (i) and (ii), it is also worth mentioning that taking $J_{11}=J_{1\bar{1}}=0$ but $V_{1}\neq 0$ will lead to another free-spin p-h asymmetric local moment fixed point (“free-ALM”), where the impurity spins completely decouple and the non-normalizable zero modes vanish. This is an unstable fixed point, because even an infinitesimal RKKY interaction induced by $J_{11}$ or $J_{1\bar{1}}$ drives the impurity spins into the $S=3/2$ or the $S=1/2$ state. Obviously, due to the three impurity spins, the free-ALM fixed point has $S_{\text{imp}}\left(T=0\right)=3\ln 2$ and $T\chi_{\text{imp}}=3/4$ relative to pristine graphene. The three fixed points F-ALM, AF-ALM and free-ALM differ only in their impurity spin states.

Cases (i) and (ii) represent the limits of maximally ferromagnetic and maximally antiferromagnetic RKKY interaction respectively. As shown in Fig. 4, when $J_{11}^{2}+J_{1\bar{1}}^{2}=\left(0.1D_{0}\right)^{2}$ and $\Lambda/D_{0}=1.5$, the K-S fixed point controls a large region of the parameter space where $J_{11}/J_{1\bar{1}}$ is not too small and $V_{1}$ is not too large. The K-I phase only occurs when the p-h symmetry is preserved and $J_{11}/J_{1\bar{1}}$ is very small (i.e. $\tilde{I}>\tilde{I}_{c0}$ where $\tilde{I}_{c0}$ is close to $1/2$)111We have obtained the phase diagram Fig. 4 using the NRG parameters $\Lambda_{\text{NRG}}=4$ and $z=1$ in the conduction band discretization scheme of Ref. Frota and Oliveira, 1986; *PhysRevB.49.11986. With these parameters we find that the system flows to the K-I fixed point when $J_{11}/J_{1\bar{1}}<0.0060$, or $\tilde{I}>\tilde{I}_{c0}\approx 0.494$. Changing $\Lambda_{\text{NRG}}$ and $z$ shifts the phase boundaries only slightly; however, where (and indeed, whether) the K-I fixed point is numerically attainable depends sensitively on the value of $z$. In the continuum limit $\Lambda_{\text{NRG}}\to 1$, it is not clear to us whether $\tilde{I}_{c0}$ is exactly $1/2$, i.e. whether the K-I phase exists only in the limiting case $J_{11}=0$., but its direct descendant– the AF-ASC phase– becomes progressively more important at intermediate values of $V_{1}$ when the RKKY interactions are not too strongly ferromagnetic. Finally, at very large values of $V_{1}$, the F-ALM and AF-ALM phases come into play, separated approximately by the $\tilde{I}=0$ line.

While the phase boundaries between K-S, AF-ASC and F-ALM seemingly meet at a single tricritical point on the $\tilde{I}$-$V_{1}$ plane, a more careful survey of the parameter space reveals the presence of another phase in a small area separating these three phases. The ground state in this phase is again p-h asymmetric with one electron removed from or added to half filling. However, in contrast to the AF-ASC phase, the ground state has a residual $S=1$ impurity, consistent with the ferromagnetic RKKY interaction ($\tilde{I}<0$); thus we name this phase “F-ASC”. Because one of the two conduction channels couples to the impurity spins and the other is blocked by potential scattering, the $S=1$ residual impurity in the ground state can form in two distinct helicities. Combined with the threefold degeneracy of the spin state, this helicity degeneracy gives an impurity entropy of $S_{\text{imp}}\left(T=0\right)=\ln 6$ relative to pristine graphene. The impurity susceptibility from the residual impurity is simply $T\chi_{\text{imp}}=2/3$.

It is also enlightening to look at the correlation between the impurity spins. In Fig. 7 we plot the expectation value $\left\langle\bm{S}_{1}\cdot\bm{S}_{2}\right\rangle$ as a function of temperature at various points in the phase diagram. At high temperatures its sign is simply opposite to that of $\tilde{I}$, as perturbation theory predicts. At low temperatures, $\left\langle\bm{S}_{1}\cdot\bm{S}_{2}\right\rangle$ takes the minimum possible value $-1/4$ in the K-I, AF-ASC and AF-ALM phases, and the maximum possible value $1/4$ in the F-ALM phase; these results are consistent with our previous analysis that the impurity spins form a spin-$1/2$ in the K-I, AF-ASC and AF-ALM phases, and a spin-$3/2$ in the F-ALM phase. The low-temperature spin correlation in the K-S phase is more interesting: $\left\langle\bm{S}_{1}\cdot\bm{S}_{2}\right\rangle$ varies smoothly from $1/4$ to [math] as $\tilde{I}$ increases from $-1$ to $\tilde{I}_{c0}$ (recall that the K-I phase takes over for $\tilde{I}>\tilde{I}_{c0}$ and $V_{1}=0$), and is only weakly dependent on $V_{1}$ as long as the system remains in the K-S phase. Therefore, away from the strongly ferromagnetic limit $\tilde{I}=-1$, the impurity spins generally form a superposition of spin-$3/2$ and spin-$1/2$ states in the K-S phase, even though the residual spin is always $1/2$, with two electrons participating in screening. A similar statement can be made for the F-ASC phase: the impurity spins form a superposition of spin-$3/2$ and spin-$1/2$ states, which couples to one conduction electron to produce a residual spin-$1$.

In the Kondo model obtained by Schrieffer-Wolff transforming the 3-impurity Anderson model, $J_{11}=J_{1\bar{1}}$ and $\tilde{I}=-1/4$; in this case it is clear from Fig. 4 that K-S and F-ALM control small- and large-$V_{1}$ physics respectively as in the strongly ferromagnetic RKKY limit, whereas the AF-ASC phase sets in for intermediate values of $V_{1}$ as in the strongly antiferromagnetic RKKY limit. For $J_{11}=J_{1\bar{1}}=0.1D_{0}/\sqrt{2}$, the AF-ASC phase occurs for $0.126<\left|V_{1}\right|/D_{0}<0.152$; this cannot be realized by a Schrieffer-Wolff transformation which requires $\left|V_{1}\right|/J_{11}\lesssim 3/4$. However, panel (b) of Fig. 4 indicates that when $J_{11}=J_{1\bar{1}}\lesssim 10^{-4}D_{0}$, the critical value of $\left|V_{1}\right|/J_{11}$ at the K-S/AF-ASC transition can be reduced dramatically, well below $3/4$. This strongly suggests that both K-S and AF-ASC (or at least their generalizations) are accessible in an Anderson model, although F-ALM and AF-ALM may still be out of reach. We will confirm this picture in Sec. IV.

We now examine the K-S/AF-ASC transition, motivated by the observation that a transition of similar nature may exist in an Anderson model. For $J_{11}=J_{1\bar{1}}=0.1D_{0}/\sqrt{2}$, the K-S/AF-ASC transition takes place at $\left|V_{1}\right|/D_{0}=v_{c}\approx 0.126$. Fig. 8 shows the high- to low-temperature crossover of $S_{\text{imp}}$ and $T\chi_{\text{imp}}$ as $V_{1}/D_{0}-v_{c}$ sweeps from $-0.01$ to $0.1$, and Fig. 9 shows the corresponding behavior of $\left\langle\bm{S}_{1}\cdot\bm{S}_{2}\right\rangle$. We can explain the unstable critical point separating the two phases as a simple level crossing of the spin-$1/2$ doublet ground states in the K-S phase and the spin singlet ground state in the AF-ASC phase. The doublet and the singlet do not mix, as they belong in different sectors of the Hilbert space. At the critical point, $S_{\text{imp}}\left(T=0\right)=-\ln\left(16/3\right)$ relative to the 4-site-vacancy graphene, which is greater than the AF-ASC value by $\ln 3$, a signature of the accidental degeneracy. Moreover, the values of both $T\chi_{\text{imp}}$ and $\left\langle\bm{S}_{1}\cdot\bm{S}_{2}\right\rangle$ at the critical point can be obtained as weighted averages of the K-S value (with weight $2/3$) and the AF-ASC value (with weight $1/3$). We mention that this simple level crossing picture also applies to the AF-CR critical point in the single-channel caseMitchell and Fritz (2013).

Finally we briefly discuss the effect of the helicity-0 channel. We assume that the irrelevant couplings $J_{00}$, $J_{01}$ and $V_{0}$ are not too large compared to the relevant couplings, so that the intermediate-coupling phase transition identified in Ref. Withoff and Fradkin, 1990 does not take place. In most cases, these irrelevant couplings merely shift the phase boundaries without modifying the phase diagram qualitatively. A notable exception is the K-I fixed point. In the p-h symmetric strongly antiferromagnetic RKKY limit $J_{11}=V_{1}=0$ and $J_{1\bar{1}}\neq 0$, we find that $V_{0}$ by itself or the combination of $J_{00}$ and $J_{01}$ does not affect the low-energy K-I behavior. However, the combination $J_{00}\neq 0$, $J_{01}=0$ and $V_{0}\neq 0$ drives the system into the K-S phase, while the combination $J_{00}=0$, $J_{01}\neq 0$ and $V_{0}\neq 0$ drives the system into the AF-ASC phase. Therefore, as with its constant-LDOS analogPaul and Ingersent (1996), the K-I phase is highly fragile against p-h symmetry breaking, and unlikely to be experimentally observed.

We summarize our results on the fixed points of the 3-impurity 3-channel Kondo model Eq. (9) in Table 2. The results for $S_{\text{imp}}\left(T=0\right)$ and $T\chi_{\text{imp}}$ are in full agreement with Figs. 5 and 6 upon changing the reference system to graphene with a 4-site vacancy, i.e. subtracting $2\ln 4$ from $S_{\text{imp}}$ and $1/4$ from $T\chi_{\text{imp}}$.

III.4 Logarithmic LDOS with an infrared cutoff

In a more realistic model of the graphene sheet, the small next-nearest-neighbor hopping $t^{\prime}$ between carbon atoms replaces the zero mode associated with a vacancy by a number of quasi-localized states shifted slightly away from the Dirac point. While the LDOS remains strongly enhanced near the energies of these quasi-localized states, it is no longer logarithmically divergentPereira et al. (2008). Nevertheless, following Ref. Mitchell and Fritz, 2013, we can fine-tune the Fermi energy to the vacancy-induced peak of the LDOS, and heuristically model the effect of a next-nearest-neighbor hopping by imposing an infrared energy cutoff $X$ on the LDOS.

To be concrete, we replace $\rho\left(\omega\right)$ by a constant $\rho\left(X\right)$ for $\left|\omega\right|<X$ in Eq. (15), so that the LDOS becomes a large constant value at small energies. Whereas such a “hard” cutoff scheme is slightly different from the “soft” cutoff adopted in Ref. Mitchell and Fritz, 2013, we can verify that the two cutoff schemes do not lead to qualitatively different results. In the $X\rightarrow 0$ limit our LDOS recovers the logarithmic divergence in the $t^{\prime}=0$ case. For our choice of the ultraviolet energy cutoff in the LDOS $\Lambda=1.5D_{0}$, we find that $X/D_{0}\sim 0.01$ reproduces the LDOS peak height found by solving the tight-binding modelPereira et al. (2008) with the experimentally estimated value $t^{\prime}\approx 0.1t$Kretinin et al. (2013). Fig. 10 shows $S_{\text{imp}}$ and $T\chi_{\text{imp}}$ in the K-S phase of the 3-impurity 2-channel Kondo model Eq. (25) as we increase the infrared cutoff $X$ from [math] to $10^{-2}D_{0}$, and Fig. 11 shows the corresponding behavior of $\left\langle\bm{S}_{1}\cdot\bm{S}_{2}\right\rangle$.

At sufficiently low energies $T\ll X$, there is no longer any contribution to $S_{\text{imp}}$ and $T\chi_{\text{imp}}$ from the non-normalizable zero mode, so both $S_{\text{imp}}$ and $T\chi_{\text{imp}}$ recover their values in the constant-LDOS spin-$3/2$ 2-channel Kondo problem in this limit, namely $\ln 2$ and $1/4$. Also, when $X$ is far smaller than $T_{K}^{X=0}$ (the Kondo temperature at $X=0$), the RG flow is still towards the K-S fixed point in the energy range $X\ll T\ll T_{K}^{X=0}$. These features are also present in the single-channel caseMitchell and Fritz (2013). On the other hand, as $X$ increases, Figs. 10 and 11 both show an increase of the total impurity spin at low energies, which is particularly pronounced for larger $X$ ($X\gtrsim 10^{-4}D_{0}$). When $X=10^{-2}D_{0}$, $\left\langle\bm{S}_{1}\cdot\bm{S}_{2}\right\rangle$ is close to $1/4$, so that the effective spin is almost completely a spin-$3/2$; this effective spin controls the physics across a wide range of energies between its formation around $T\sim 10^{-3}D_{0}$ and the onset of screening below $T\sim 10^{-6}D_{0}$.

IV Anderson model

Having discussed the 3-impurity Kondo model in great detail, in this section we turn back to our initial approximation of ignoring the hydrogen impurity and the central A site. This approximation is based on an infinite hydrogen-carbon coupling strength. Realistic estimates put the hydrogen-carbon coupling around twice the nearest-neighbor hopping between carbon atomsRobinson et al. (2008); Wehling et al. (2010); Gmitra et al. (2013); it is therefore important to check whether our intuitions from the Kondo model carry over to the full 5-atom cluster Anderson model. This Anderson model is also represented by the Hamiltonian

[TABLE]

$H_{\text{vac}}$ and $H_{\text{hyb}}$ are already given in Eqs. (2) and (4). The impurity Hamiltonian has additional terms:

[TABLE]

Here we have labeled the hydrogen impurity as $g$, the central A site $a(\vec{0})$ as $a_{0}$, and defined $n_{a,0\alpha}\equiv a_{0\alpha}^{{\dagger}}a_{0\alpha}$, $n_{H\alpha}\equiv g_{\alpha}^{{\dagger}}g_{\alpha}$, $n_{a,0}=n_{a,0\uparrow}+n_{a,0\downarrow}$ and $n_{H}=n_{H\uparrow}+n_{H\downarrow}$. Compared to Eq. (1), Eq. (26) has a number of new coupling constants: the on-site chemical potentials $\epsilon_{a}$ and $\epsilon_{H}$, the Hubbard interaction on the hydrogen impurity $U_{H}$, the hydrogen-carbon coupling strength $t_{H}$, and the nearest-neighbor hopping between the central A site and its nearest neighbors $t_{0}$. $t_{0}$ is generally different from $t$ due to the presence of the hydrogen impuritySofo et al. (2012). We also note that, due to the two additional impurity sites, Eq. (26) cannot be mapped to the simple Kondo model Eq. (9) even in the limit $U\sim\left|\epsilon_{b}\right|\gg t$.

Because of the immense size of the parameter space, we now focus on the experimentally relevant regime where all parameters (including $U$ and $U_{H}$) are of comparable magnitudes. We also continue to neglect the helicity-0 channel with a linear LDOS and keep only the helicity-$\pm 1$ channels with a logarithmically divergent LDOS. Under these assumptions, quite generally, we find that the ground state of the system is a charge-neutral spin doublet state when the p-h symmetry breaking terms are weak, or a spin singlet state with charge $+1$ (or $-1$) when the p-h symmetry breaking terms are strong. For reasons that will become clear later we again call these two phases K-S and AF-ASC respectively.

To be more concrete, we choose $U=t$, $U_{H}=2.8t$, $\epsilon_{b}=-U/2$, $t_{H}=2t$ and $t_{0}=0.6t$. When $\epsilon_{a}=-U/2-0.7t$ and $\epsilon_{H}=-U_{H}/2$, a Hartree-Fock calculation of the LDOS in the full Anderson-Hubbard model (where the Hubbard interaction is also included in $H_{\text{vac}}$) has been reported to agree qualitatively with density-functional theory resultsSofo et al. (2012). However, we argue that the p-h symmetry breaking should be stronger on the hydrogen impurity than on the central A site. In the following we therefore let $\epsilon_{a}=-U/2$ and vary $\epsilon_{H}$ instead, placing the p-h symmetry breaking term on the hydrogen impurity. We nevertheless note that our results below are not qualitatively modified by the presence of additional p-h symmetry breaking terms on the central A site or its nearest-neighboring B sites, as long as these terms are not too large compared to $t$.

Fig. 12 shows the typical behavior of $S_{\text{imp}}$ and $T\chi_{\text{imp}}$ in the K-S and the AF-ASC phases and across the phase transition in between. The K-S/AF-ASC transition occurs at $\left(\epsilon_{H}+U_{H}/2\right)/t=\tilde{\epsilon}_{c}\approx-1.055$, and we tune $\left(\epsilon_{H}+U_{H}/2\right)/t$ from [math] to $-2$. The low-temperature behavior of $S_{\text{imp}}$ and $T\chi_{\text{imp}}$ is completely identical to that of their Kondo model counterparts, not only inside each phase but also at the transition (cf. Figs. 5 and 8). The K-S/AF-ASC transition can again be explained as a simple level crossing of the spin-$1/2$ doublet ground states in the K-S phase and the spin singlet ground state in the AF-ASC phase. It is also interesting to consider the equal-time spin correlations $\left\langle\bm{S}_{1}\cdot\bm{S}_{2}\right\rangle$ and $\left\langle\bm{S}_{1}\cdot\bm{S}_{H}\right\rangle$, where $\bm{S}_{1}$ is again the spin on the nearest-neighbor B site $b_{1}$, and $\bm{S}_{H}\equiv\frac{1}{2}\sum_{\alpha\beta}g_{j\alpha}^{\dagger}\bm{\sigma}_{\alpha\beta}g_{j\beta}$ is the spin on the hydrogen impurity; these are plotted in Fig. 13 for different values of $\epsilon_{H}$. We see that $\left\langle\bm{S}_{1}\cdot\bm{S}_{2}\right\rangle$ is ferromagnetic in the K-S phase, but becomes antiferromagnetic in the AF-ASC phase (although now much smaller in magnitude than $-1/4$), changing sign across the phase transition as in the Kondo model (cf. Figs. 7 and 9). A similar behavior is seen in $\left\langle\bm{S}_{1}\cdot\bm{S}_{H}\right\rangle$, i.e. the spin correlation between the hydrogen impurity and the nearest-neighbor B sites also goes from ferromagnetic to antiferromagnetic across the K-S/AF-ASC transition. On the other hand, the spin correlations involving the spin on the central A site, $\bm{S}_{0}\equiv\frac{1}{2}\sum_{\alpha\beta}a_{0\alpha}^{\dagger}\bm{\sigma}_{\alpha\beta}a_{0\beta}$, are almost unchanged at the transition.

Further information of the K-S and AF-ASC phases is given in Fig. 14, where we plot the low-energy expectation values of the following operators: the total spin on the nearest neighbor B sites $(\bm{S}_{1}+\bm{S}_{2}+\bm{S}_{3})^{2}$, the total spin on the hydrogen impurity and the central A site $(\bm{S}_{0}+\bm{S}_{H})^{2}$, the total spin on the 5-atom cluster $\bm{S}_{\text{tot}}^{2}\equiv(\bm{S}_{1}+\bm{S}_{2}+\bm{S}_{3}+\bm{S}_{0}+\bm{S}_{H})^{2}$, the occupancy of a nearest-neighbor B site $n_{b,1}$, the occupancy of the central A site $n_{a,0}$, and the occupancy of the hydrogen impurity $n_{H}$. $(\bm{S}_{0}+\bm{S}_{H})^{2}$ is almost unchanged across the phase transition and remains small (less than $1/4$) for $-2t\leq\epsilon_{H}+U_{H}/2\leq 0$, while both $(\bm{S}_{1}+\bm{S}_{2}+\bm{S}_{3})^{2}$ and $\bm{S}_{\text{tot}}^{2}$ fall dramatically when the system goes from the K-S phase to the AF-ASC phase. While $(\bm{S}_{1}+\bm{S}_{2}+\bm{S}_{3})^{2}$ decreases from $1.6$ to $1.0$, $\bm{S}_{\text{tot}}^{2}$ undergoes a sharper drop from $1.7$ to $0.8$, which is a natural result of $\left\langle\bm{S}_{1}\cdot\bm{S}_{H}\right\rangle$ changing from ferromagnetic in the K-S phase to antiferromagnetic in the AF-ASC phase. Recalling that the K-S phase has total spin $1/2$ and the AF-ASC phase has total spin [math], we conclude that the spin of the 5-atom cluster is Kondo-screened by the conducting channels in both the K-S phase and the AF-ASC phase. We stress that, in contrast to the linear-LDOS caseSofo et al. (2012), such a Kondo effect does not require a very strong coupling between the 5-atom cluster and the surrounding non-interacting bath. Meanwhile, $n_{b,1}$, $n_{a,0}$ and $n_{H}$ all experience a sudden increase across the transition as $\epsilon_{H}<-U_{H}/2$ decreases, which is consistent with the total charge in the ground state increasing by one; $n_{b,1}$ increases by the largest amount of the three, from $0.98$ to $1.2$.

Finally, we briefly discuss the effect of an infrared energy cutoff $X$ of the LDOS on the Anderson model. Fig. 15 shows $S_{\text{imp}}$ and $T\chi_{\text{imp}}$ for different values of $X$ in both the K-S phase and the AF-ASC phase. Not too surprisingly, the behavior of $S_{\text{imp}}$ and $T\chi_{\text{imp}}$ approaches the $X=0$ case at intermediate energy scales $X\ll T\ll T_{K}^{X=0}$, while at lower energies $T\ll X$, $S_{\text{imp}}$ and $T\chi_{\text{imp}}$ return to their values in the constant-LDOS version of the model where there is no contribution from non-normalizable zero modes. While increasing $X$ also leads to an increase in magnitude for the spin correlations $\left\langle\bm{S}_{1}\cdot\bm{S}_{2}\right\rangle$ and $\left\langle\bm{S}_{1}\cdot\bm{S}_{H}\right\rangle$, this is a tiny effect in comparison with the spin correlation of the Kondo model shown in Fig. 11.

To summarize this section, our results indicate that the behavior of the Anderson model which contains the hydrogen impurity and the four nearest carbon atoms is qualitatively captured by the 3-impurity Kondo model. There exist a p-h symmetric spin-$1/2$ K-S phase where the impurity spins align ferromagnetically, and a p-h asymmetric spin-singlet AF-ASC phase where the impurity spins align antiferromagnetically. It is possible that the K-S/AF-ASC transition picture is applicable to even more realistic models of the hydrogen impurity: weaker p-h symmetry breaking favors ferromagnetic spin correlation and leads to magnetic ground states, while stronger p-h symmetry breaking favors antiferromagnetic spin correlation and tends to suppress the ground state degeneracy.

V Conclusions and outlook

In this paper, we have studied the Kondo effect associated with a single hydrogen impurity on graphene. The hydrogen impurity is strongly coupled to the “central” carbon atom directly below it. First we consider the limit of infinite coupling, so that the hydrogen atom and the central carbon atom are effectively decoupled from the rest of the system, but the $C_{3}$ rotation symmetry of the system is preserved. To model the induced magnetization, we place a strong Hubbard interaction on the three nearest neighbor carbon atoms, creating three magnetic impurities. The remaining graphene sheet with four vacancy sites, approximated to be nearest-neighbor and non-interacting, supports two conduction channels which hybridizes with the three impurities with a local density of states diverging logarithmically as a function of energy near the Dirac point, in addition to a conduction channel whose LDOS vanishes linearly.

We study the resulting 3-impurity, 3-channel Kondo model with the numerical renormalization group method. Couplings to the conduction channel with a linear LDOS are irrelevant and usually negligible, and the phase diagram is controlled by the Kondo and potential scattering coupling constants associated with the two conduction channels with a logarithmically divergent LDOS. The regime where the potential scattering is not too strong sees the competition between a p-h symmetric Kondo phase (K-S) and a p-h asymmetric strong-coupling phase (AF-ASC). Ferromagnetic RKKY interactions between the magnetic impurities and weaker potential scattering favor the K-S phase, where the p-h symmetric ground state is a residual spin-$1/2$ after screening by the two conduction channels, and the impurity spins tend to align ferromagnetically. On the other hand, antiferromagnetic RKKY interactions and stronger potential scattering favor the AF-ASC phase, where the ground state is a spin singlet with one electron removed from or added to half filling, and the impurity spins align antiferromagnetically. In the strong potential scattering regime, the potential scattering coupling strength renormalizes to infinity, and the magnetic impurities decouple from the conduction channels, forming a local moment whose size depends on the RKKY interactions.

Relaxing the approximation of infinite hydrogen-carbon hybridization, we obtain an Anderson model with 5 impurity sites: the hydrogen atom, the central carbon atom and its three nearest neighbors in the tight-binding model. For realistic Hubbard interaction strengths on impurities, we find through NRG that the ground state is the p-h symmetric spin-$1/2$ K-S phase when the p-h symmetry breaking is not too strong, and the particle-hole asymmetric spin singlet AF-ASC phase otherwise. Kondo screening is shown to take place in both phases of the Anderson model. In the K-S phase, the spins of the nearest neighbor carbon atoms align ferromagnetically with each other and with the spin of the hydrogen atom, whereas in the AF-ASC phase they align antiferromagnetically with each other and with the hydrogen spin.These provide evidence that our 3-impurity Kondo model approximation is qualitatively reasonable.

Many open questions remain to be answered. First of all, we have assumed throughout this work that the bulk chemical potential is fine-tuned to the singularity of the vacancy-induced logarithmically divergent LDOS, which coincides with the zero point of the bulk density of states. Perturbations at the LM fixed point are thus strongly relevant for the helicity-1, 1̄ channels with a logarithmically divergent LDOS, and strongly irrelevant for the helicity-0 channel with a linear LDOS. If the bulk chemical potential is shifted by an applied gate voltage, it is interesting to check whether the helicity-0 channel will have a progressively more important influence on the Kondo temperature and the transport properties of the system, as one may expect from results on the single-channel Kondo problem in gated or doped grapheneCornaglia et al. (2009); Vojta et al. (2010); Lo et al. (2014); Ruiz-Tijerina and Dias da Silva (2017); May et al. (2018).

As a closely related point, in obtaining the 5-impurity Anderson model, we have neglected the next-nearest-neighbor hopping between carbon atoms. While imposing an infrared cutoff on the logarithmically divergent LDOS partially mimics its effectsPereira et al. (2008), the next-nearest-neighbor hopping will also change the wave functions of the bulk conduction electrons and their coupling to the hydrogen impurity. A more careful treatment of the next-nearest-neighbor hopping is thus necessary for a quantitative comparison with experiments.

The electron-electron interaction on carbon atoms farther away from the hydrogen impurity than the three nearest neighbors is another essential ingredient in a more realistic model, since experiments have shown that the spin-polarized state induced by the hydrogen impurity has a large spatial extensionGonzález-Herrero et al. (2016). Such interactions have been taken into account in previous studies within the Hartree-Fock approximationKumazaki and Hirashima (2007); Sofo et al. (2012) and using dynamical mean-field theoryHaase et al. (2011); Charlebois et al. (2015). To the best of our knowledge, it is not clear how the vacancy-induced non-normalizable zero modes behave in the presence of bulk electron-electron interactions when they are not strong enough to turn graphene into a Mott insulator. Addressing this issue will be useful for a theoretical understanding of the unusually long-ranged coupling between the magnetic moments induced by different hydrogen adatomsGonzález-Herrero et al. (2016); Buividovich et al. (2017).

Acknowledgements.

We thank Eva Andrei, Josh Folk and Eran Sela for helpful discussions. The research of IA was supported by NSERC Discovery Grant 04033-2016 and the Canadian Institute for Advanced Research.

Appendix A Conduction channels at low energies

In this appendix we derive the low-energy behavior of the three conduction channels of definite helicities, Eqs. (13a), (13b) and (13c). This is achieved by diagonalizing the non-interacting Hamiltonian with a 4-site vacancy $H_{\text{vac}}$ and finding the scattering state wave functions.

The vacancy can be implemented by strong potential scattering: $H_{\text{vac}}=H_{0}+V_{i}$, with $H_{0}$ describing the translationally invariant pristine graphene,

[TABLE]

and $V_{i}$ simulating the vacancy,

[TABLE]

The limit $U_{2}\rightarrow\pm\infty$ corresponds to the vacancy sites $a(\vec{0})$, $b(\vec{0})$, $b\left(-\vec{a}_{2}\right)$ and $b\left(-\vec{a}_{1}\right)$. Physically we do not expect the value of $U_{1}$ to affect the scattering state wave functions in the $U_{2}\rightarrow\pm\infty$ limit, because $a(\vec{0})$ will be isolated from the other sites.

We will work in the basis of the $H_{0}$ eigenstates. These are given by

[TABLE]

here we have performed the Fourier transform

[TABLE]

where $S_{1BZ}=8\pi^{2}/(\sqrt{3}a^{2})$ is the area of the hexagonal first Brillouin zone. It is straightforward to rewrite $H_{0}$ and $V_{i}$ in terms of $\psi_{\pm}$,

[TABLE]

with the dispersion

[TABLE]

and

[TABLE]

The scattering states are given by

[TABLE]

where for $\lambda$, $\lambda^{\prime}=\pm 1$, $G_{\lambda,\lambda^{\prime}}\left(\vec{p},\vec{p}^{\prime}\right)\left(\lambda\epsilon_{\vec{p}}-\lambda^{\prime}\epsilon_{\vec{p}^{\prime}}\right)$ is finite. By definition $\phi$ diagonalizes $H_{\text{vac}}$:

[TABLE]

The energy eigenvalues are those of the $H_{0}$ eigenstates because the impurity is localized. This equation is then solved for $G_{\lambda,\lambda^{\prime}}\left(\vec{p},\vec{p}^{\prime}\right)$.

In the limit of $U_{2}\rightarrow\pm\infty$, the scattering states are indeed independent of $U_{1}$:

[TABLE]

The negative-energy states are found by interchanging $+$ with $-$, and inverting the signs of all absolute values. We have introduced the shorthand

[TABLE]

It is useful to give a low-energy asymptotic formula for $L\left(\omega^{+}\right)\equiv L\left(\omega+i0\right)$, valid for $\left|\omega\right|\ll\Lambda\sim t$, obtained by only keeping contributions from near the two Dirac points:

[TABLE]

where $\Lambda_{0}=\Lambda e^{\frac{\pi}{6\sqrt{3}}}$ is another ultraviolet energy cutoff.

Using Eqs. (31) and (30), we can rewrite $a$ in terms of $\phi$:

[TABLE]

where the $\epsilon_{\vec{p}}\rightarrow-\epsilon_{\vec{p}}$ part is the contribution from the negative energy eigenstates $\phi_{-,\vec{p}}$, and we have further defined

[TABLE]

[TABLE]

It is useful to note that $L\left(\vec{R},z\right)$ have all the symmetries of the hexagonal lattice, and that $L\left(\vec{a}_{1},z\right)$ is related to $L\left(z\right)$ by

[TABLE]

According to Eq. (40), the symmetric linear combinations $a_{1,2,3}$ have the form

[TABLE]

[TABLE]

[TABLE]

For the helicity-$\pm 1$ combinations $c_{h=1}$ and $c_{h=\bar{1}}$, at low energies it is permissible to keep only the terms that are logarithmically divergent at the Dirac points:

[TABLE]

[TABLE]

Here $\phi_{\vec{K},\pm,\vec{k}}\equiv\phi_{\pm,\vec{K}+\vec{k}}$. On the other hand, for the helicity-[math] linear combination $c_{h=0}$, the divergent terms are suppressed by $O\left(k^{2}\ln k\right)$ at low energies, and the constant incident terms dominate instead:

[TABLE]

We now take advantage of the rotational invariance at low energies, and introduce angular momentum eigenmodes labeled by the quantum number $m$:

[TABLE]

These eigenmodes $\tilde{\phi}$ obey

[TABLE]

and in terms of $\tilde{\phi}$,

[TABLE]

Inserting Eq. (48) into Eqs. (45), (46) and (47) then yields Eqs. (13a), (13b) and (13c).

Appendix B Non-normalizable zero modes

This appendix elaborates on the zero modes of the infinite graphene sheet with 4 vacancy sites. As discussed in Sec. II, these non-normalizable zero modes are responsible for the logarithmic divergence in the LDOS of our impurity models. We will solve the lattice Schrödinger equation by generalizing the method of Ref. Pereira et al., 2006, give the long-distance asymptotics of the two solutions, and briefly discuss their fate in the strong-coupling regime of the impurity models.

It is convenient to relabel the lattice sites as in Fig. 16, with the two solid red lines dividing the plane into three parts: the left half-plane with a zigzag edge, the right half-plane with a “bearded” edge, and the middle strip that contains the four vacancy sites. The zero mode wave functions vanish on the entire B sublattice, so we focus on the wave function on the A sublattice, which we denote as $\phi_{l,j}$; here $l$ is an integer and $j$ is either an integer or a half-integer, but $l+2j$ is always even.

Away from the vacancy, the Schrödinger equation at zero energy reads

[TABLE]

this allows the expansion of the zero mode wave function on the left half-plane in edge states of the zigzag edgeWakabayashi et al. (1999),

[TABLE]

as well as the expansion on the right half-plane in edge states of the bearded edge,

[TABLE]

Inserting these expansions into the $l=0$ and $l=1$ equations and eliminating $\phi_{0,j}$, we have

[TABLE]

which is true for any integer $j$ as long as $j\neq 0$ and $j\neq-1$. Using the relation

[TABLE]

valid for nonzero integer $j$, we find two nontrivial solutions by inspection:

[TABLE]

and

[TABLE]

These solutions are linearly independent, and are therefore the only zero energy solutions allowedPereira et al. (2008).

We can show the long-distance asymptotic behavior of these solutions is given by

[TABLE]

and

[TABLE]

where $x=\sqrt{3}la/2$, $y=ja$, and $r=\sqrt{x^{2}+y^{2}}\rightarrow\infty$. For instance, on the right half-plane, the first solution

[TABLE]

is dominated by momenta near $k^{\prime}a=\pm 2\pi/3$ when $l\gg 1$:

[TABLE]

where $\tilde{\Lambda}$ is a momentum cutoff of $O\left(1/a\right)$. Performing the integrals and taking the $\tilde{\Lambda}\rightarrow\infty$ limit, we promptly obtain Eq. (58). Eqs. (58) and (59) suggest that the linear combinations $e^{-i\frac{\pi}{3}}\phi^{\left(1\right)}-e^{i\frac{\pi}{3}}\phi^{\left(2\right)}$ and $e^{i\frac{\pi}{3}}\phi^{\left(1\right)}-e^{-i\frac{\pi}{3}}\phi^{\left(2\right)}$ are eigenstates of the $C_{3}$ rotation, which is indeed the case: the former has helicity 1 and the latter has helicity 1̄.

Since $\phi^{\left(1\right)}$ and $\phi^{\left(2\right)}$ are linearly independent, we are unable to construct a normalizable zero mode whose wave function drops to zero faster than $1/r$ as $r=\sqrt{x^{2}+y^{2}}\rightarrow\infty$. However, if we consider removing even more sites from the graphene lattice, more than two zero modes may be allowed. (The simplest example of removing a site together with its three nearest neighbors and six next nearest neighbors produces $\left|1+6-3\right|=4$ zero modes.) In such a situation, at most two zero modes decaying as $1/r$ are linearly independent, which we can choose as $1/\left(x+iy\right)$ and $1/\left(x-iy\right)$; any other zero mode can be combined with these two non-normalizable modes to yield a wave function that decays faster than $1/r$. In other words, at most two conduction channels have a logarithmically divergent LDOS.

We conclude this appendix by explaining why the zero modes cease to exist in the “strong-coupling lattice”. The removal of $c_{1}$ and $c_{\bar{1}}$ from the original lattice amounts to the condition that the corresponding wave functions vanish. We can directly calculate these wave functions for the two solutions:

[TABLE]

It is easy to verify that $c_{\bar{h}}$ vanishes for the helicity-$h$ zero mode $e^{-ih\frac{\pi}{3}}\phi^{\left(1\right)}-e^{ih\frac{\pi}{3}}\phi^{\left(2\right)}$. However, $c_{h}$ does not vanish for the helicity-$h$ solution, which means the electronic states $c_{1}$ and $c_{\bar{1}}$ cannot be projected out without removing both zero modes.

Appendix C RKKY interaction in graphene with a 4-site vacancy

In this appendix we calculate the RKKY interaction between magnetic impurities in the Kondo model Eq. (9) to the second order in the Kondo couplings. We show that the RKKY interaction at low temperatures is dominated by the helicity $1$ and $\bar{1}$ channels, and remains ferromagnetic despite the presence of the vacancy.

Following Ref. Kogan, 2011, to the second order in Kondo couplings $J_{hh^{\prime}}$, we can write the RKKY interaction between $b_{1}$ and $b_{2}$ as

[TABLE]

where the (isothermal) static spin susceptibilities $\chi_{hh^{\prime}}$ are evaluated using Wick’s theorem for the non-interacting Hamiltonian with a 4-site vacancy $H_{\text{Vac}}$,

[TABLE]

The factor of $1/4$ comes from spin degrees of freedom, and $\beta=1/T$. The imaginary time Green’s function $G^{c}$ is defined by

[TABLE]

$G^{c}$ is diagonal in the helicity index, and may be expressed as linear combinations of the real space Green’s function $G_{aa}\left(\vec{R},\vec{R}^{\prime},\tau\right)\equiv-\left\langle T_{\tau}a\left(\vec{R},\tau\right)a^{{\dagger}}\left(\vec{R}^{\prime},0\right)\right\rangle$.

We proceed to find $G_{aa}$ by solving its equation of motion (coupled with that of $G_{ba}\left(\tau\right)\equiv-\left\langle T_{\tau}b\left(\tau\right)a^{{\dagger}}\left(0\right)\right\rangle$) in momentum spacePeres et al. (2009); Nanda et al. (2012). The result is

[TABLE]

[TABLE]

we find $G_{hh}^{c}$ in particularly simple forms:

[TABLE]

[TABLE]

The low-energy expressions of these Green’s functions can also be found from Eqs. (13a), (13b) and (13c). One can show, term by term, that $G_{hh}^{c}\left(z\right)$ is analytic everywhere except on the real axis.

We are ready to compute $\chi_{hh^{\prime}}$:

[TABLE]

where $n_{F}\left(\omega\right)=1/\left(e^{\beta\omega}+1\right)$, $\omega^{\pm}\equiv\omega\pm i0$, and we have deformed the contour of integration into two straight lines $\operatorname{Im}z=\pm 0^{+}$.

While $\chi_{00}$ and $\chi_{01}=\chi_{0\bar{1}}=\chi_{10}=\chi_{\bar{1}0}$ are finite, it turns out that $\chi_{11}=\chi_{1\bar{1}}=\chi_{\bar{1}1}=\chi_{\bar{1}\bar{1}}$ is divergent for temperatures $T\ll\Lambda$:

[TABLE]

Inserting this into Eq. (63), we find that at $T_{K}\ll T\ll\Lambda$, the RKKY interaction can be approximated as

[TABLE]

where $J_{hh^{\prime}}\left(T\right)$ are the renormalized Kondo couplings at energy scale $T$.