# Local multiplet formation around a single vacancy in graphene: an effective
Anderson model analysis based on the block-Lanczos DMRG method

Tomonori Shirakawa  Computational Materials Science Research Team, RIKEN
Center for Computational Science (R-CCS), Kobe, Hyogo 650-0047, Japan  Quantum
Computational Science Research Team, RIKEN Center for Quantum Computing (RQC),
Wako, Saitama 351-0198, Japan  Seiji Yunoki  Computational Materials Science
Research Team, RIKEN Center for Computational Science (R-CCS), Kobe, Hyogo
650-0047, Japan  Quantum Computational Science Research Team, RIKEN Center for
Quantum Computing (RQC), Wako, Saitama 351-0198, Japan  Computational Quantum
Matter Research Team, RIKEN Center for Emergent Matter Science (CEMS), Wako,
Saitama 351-0198, Japan  Computational Condensed Matter Physics Laboratory,
RIKEN Cluster for Pioneering Research (CPR), Saitama 351-0198, Japan

###### Abstract

To better understand the electronic structure of a single vacancy in graphene,
we study the ground state property of an effective Anderson model, consisting
of three dangling $sp^{2}$ orbitals of the surrounding carbon atoms around the
vacancy and the $\pi$ orbitals of carbon atoms that form the honeycomb lattice
with a single vacancy. This model possesses $\mathcal{C}_{3}$ point group
symmetry around the vacancy and thereby the local multiplets can be
characterized by their irreducible representations. Employing the block-
Lanczos density-matrix renormalization group (DMRG) scheme proposed by the
present authors [T. Shirakawa and S. Yunoki, Phys. Rev. B 90, 195109 (2014)],
we show that there are two phases in the relevant parameter space, i.e., a
nonmagnetic phase in the weak coupling region and a free magnetic moment phase
in the realistic parameter region. The systematic analysis finds that, in the
free magnetic moment phase, local multiplets of the doubly degenerate
irreducible representation with spin 1 become dominant in the ground state,
and approximately half of this local spin 1 is screened by electrons in the
surrounding $\pi$ orbitals, indicating the emergence of the residual spin-1/2
free magnetic moment. The symmetry of this local multiplets is compatible with
the occurrence of the in-plane Jahn-Teller distortion to lift the degeneracy,
found in the previous ab-initio calculations based on the density functional
theory. Furthermore, we find that the emergence of the free magnetic moment is
robust against carrier doping, which is in sharp contrast to the case of
graphene with an adatom, thus explaining the qualitative difference observed
experimentally in these two classes of systems. We also find the enhancement
of the spin correlation function between $\pi$ electrons around the vacancy
and those in the conduction band away from the vacancy in the undoped case, as
compared to that in the doped case, while the spin correlation function
between the $\sigma$ electrons in the $sp^{2}$ dangling orbitals around the
vacancy and the $\pi$ electrons in the conduction band remains large in both
undoped and doped cases. This implies that there is an additional contribution
for the free magnetic moment from the $\pi$ electrons, which is fragile
against the carrier doping, besides the free magnetic moment due to the
$\sigma$ electrons in the $sp^{2}$ dangling orbitals, which is robust against
the carrier doping. Our calculations thus support qualitatively the previous
experiment that suggests the emergence of free magnetic moment with two
distinct origins.

##  I Introducion

Graphene with vacancies has attracted a great deal of attention because it
provides a route for additional functionality of graphene, i.e., the emergence
of magnetism as predicted by ab-initio calculations based on the density-
functional theory (DFT) Yazyev and Helm (2007); Dharma-wardana and Zgierski
(2008); Dai  _et al._ (2011); Paz  _et al._ (2013); Padmanabhan and Nanda
(2016); Valencia and Caldas (2017), similar to the case of graphite El-Barbary
_et al._ (2003); Lehtinen  _et al._ (2004). Indeed, the magnetization
measurements have observed paramagnetism in the presence of vacancies Ney  _et
al._ (2011); Nair  _et al._ (2012, 2013), while the pristine graphene is
diamagnetic Sepioni  _et al._ (2010). Moreover, the tunnel scanning microscope
measurements have observed the spectra showing two spin-polarized peaks
located in the vicinity of vacancy Zhang  _et al._ (2016).

The introduction of a single vacancy in graphene induces dangling orbitals
that are localized around the vacancy. These orbitals would contribute to
reconstruct the electronic structure around the vacancy because the energy
levels of these dangling orbitals are located near the Fermi level. Therefore,
the correlation effect is expected to play an important role to determine the
electronic structure even in the carbon based system.

The previous ab-initio DFT calculations have shown that there occurs the in-
plane Jahn-Teller distortion around the vacancy Yazyev and Helm (2007);
Dharma-wardana and Zgierski (2008); Dai  _et al._ (2011); Paz  _et al._
(2013); Padmanabhan and Nanda (2016), in which two of three carbon atoms
around the vacancy become closer in distance. Let us briefly consider the
spatial symmetry of graphene in the presence of a single vacancy. If one of
carbon atoms forming the honeycomb lattice structure is removed from graphene,
the system has $\mathcal{C}_{3}$ point group symmetry around the vacancy.
Under the $\mathcal{C}_{3}$ rotational symmetry, the eigenstates of the system
are characterized with the non-degenerate symmetric irreducible representation
$\mathcal{A}$ or the doubly degenerate irreducible representation
$\mathcal{E}$. The distortion found in the ab-initio DFT calculations is
considered as the consequence of the coupling through the vibration of
$\mathcal{E}$-type, called $e$ vibration mode. The electronic state coupled to
this vibration mode is expected to be a state in $\mathcal{E}$ symmetry
Casartelli  _et al._ (2013). This can indeed be compatible, for example, with
a noninteracting state composed of the three $sp^{2}$-dangling orbitals around
the vacancy occupied by three electrons. However, if we consider a
noninteracting system composed of three $\pi$ orbitals around the vacancy in
addition to the three $sp^{2}$-dangling orbitals, there are six electrons and
the ground state can be in $\mathcal{A}$ symmetry. Furthermore, it is not
trivial how the correlation effect along with the surrounding other orbitals
can be incorporate in these pictures.

Another issue is the size of the free magnetic moment induced by the presence
of vacancy in graphene. It has been reported that the magnetization curve of
graphene in the presence of vacancy generated by the proton irradiation is
well described by the Brillouin function with $1.0\mu_{B}$ Nair  _et al._
(2012, 2013) ($\mu_{B}$ being the Bohr magneton), while the other experiment
for graphene nanoflakes irradiated with nitrogen ions suggests $2.0\mu_{B}$
paramagnetism induced per single vacancy Ney  _et al._ (2011). Since graphene
with vacancy is a complex system to analyze experimentally, the theoretical
investigation based on the numerical approaches is highly anticipated to be
useful. However, the estimated values by the ab-initio DFT calculations are
rather widely distributed from $1.0\mu_{B}$ to $2.0\mu_{B}$ Yazyev and Helm
(2007); Dharma-wardana and Zgierski (2008); Dai  _et al._ (2011); Paz  _et
al._ (2013); Padmanabhan and Nanda (2016), depending on the use of different
types of the exchange potential Valencia and Caldas (2017) as well as the
treatment of the further corrections for the electron correlation Casartelli
_et al._ (2013); Jr  _et al._ (2019). Note also that the size of the free
magnetic moment observed experimentally could not be simply estimated
theoretically by the static one-body approximation because it fails to
describe multiplet electronic structures and Kondo-like physics. Therefore, it
is valuable to address these issues based on the many-body model calculations.

Moreover, other experiments suggest that the emergence of the free magnetic
moment around the vacancy is robust against carrier doping Nair  _et al._
(2013); Zhang  _et al._ (2016). This is in sharp contrast to the magnetism
induced by hydrogen absorption, where the magnetic moment vanishes by carrier
doping Nair  _et al._ (2013); González-Herrero  _et al._ (2016). It is also
claimed in Ref. Nair  _et al._ (2013) that the vacancy magnetism in graphene
has a dual origin and approximately half of the free magnetic moment abruptly
diminishes with carrier doping. Therefore, the mechanism for the emergence of
the free magnetic moment in graphene would be different between these two
systems with vacancies and adatoms.

To better understand the origin of magnetism, here we examine the ground state
property of an effective Anderson model for a single vacancy in graphene.
Although there have been several models proposed for graphene in the presence
of vacancy Kanao  _et al._ (2012); Cazalilla  _et al._ (2012); Mitchell and
Fritz (2013), the essential difference of our model is to contain all three
$sp^{2}$ dangling orbitals around the vacancy, in addition to the $\pi$
orbitals of carbon atoms that form the honeycomb lattice with a single
vacancy. We also impose the $\mathcal{C}_{3}$ rotational symmetry around the
vacancy to discuss a possible mechanism for the Jahn-Teller distortion.

First, we analyze the local multiplet structures around a vacancy based on a
three-site two-orbital Hubbard model, composed of the three dangling $sp^{2}$
orbitals and three $\pi$ orbitals closest to the vacancy, which will be
incorporated into an effective Anderson model as the impurity sites. We show
that there are two phases in a reasonable parameter region: one is a weak
coupling phase in which the multiplet structures are obtained by the
perturbation theory to the open shell electron configuration realized in the
noninteracting limit, and the other is a strong coupling phase in which each
of carbon sites forms spin triplet due to the Hund’s coupling. The effective
model in the strong coupling region is thus well described by the spin-1
antiferromagnetic Heisenberg model whose ground state is often referred to as
a valence bond state Haldane (1983); Affleck  _et al._ (1987). In both phases,
the low-lying multiplet structures are identified by the irreducible
representations of the $\mathcal{C}_{3}$ rotational point group and the total
spin.

We then examine how the local multiplet ground state is affected via the
hybridization with the surrounding $\pi$ electrons in the conduction band. For
this purpose, we construct an effective Anderson model composed of the three-
site two-orbital Hubbard model, which describes the local electronic states
around the vacancy and thus serves as the impurity sites, and the surrounding
$\pi$ orbitals of carbon atoms forming the honeycomb lattice with the single
vacancy. The Anderson model is solved by employing the density matrix
renormalization group (DMRG) method White (1992, 1993) combining with the
block-Lanczos technique Shirakawa and Yunoki (2014); Allerdt  _et al._ (2015,
2017); Allerdt and Feiguin (2019), which gives the numerically accurate
solution within controlled errors. With this numerical analysis, we can fully
take into account the local multiplet structures as well as the coupling to
the surrounding $\pi$ electrons, which thus provides the valuable information
complementary to the one-body type approximation used in the ab-initio DFT
calculations.

Main results are summarized as follows. We find that (1) in the reasonable
parameter region for graphene, the local multiplets with $\mathcal{E}$
irreducible representation and local spin $S=1$ are dominant in the ground
state, due to the coupling with the surrounding $\pi$ electrons in the
conduction band. The symmetry of the local $\mathcal{E}$ multiplets is thus
compatible to the occurance of the in-plane Jahn-Teller distortion by the
coupling through the $\mathcal{E}$-type vibration mode Casartelli  _et al._
(2013). We also find that (2) the local magnetic moment is not completely
screened but partially screened by the surrounding $\pi$ electrons in the
conduction band, suggesting the existence of the residual free magnetic
moment. The estimated local free magnetic moment around the vacancy is as
large as $\sim 0.55$ for the undoped case and $\sim 0.45$ for the doped case,
where the g factor $g=2$ is assumed, implying the $S=1/2$ paramagnetism formed
locally around the vacancy. We also confirm that (3) this free magnetic moment
is robust against carrier doping, which is compatible to the experimental
results Nair  _et al._ (2013); Zhang  _et al._ (2016). Finally, analyzing the
spin correlation function between the electrons in the local impurity sites
and the surrounding $\pi$ electrons in the conduction band, we find that (4)
the spin correlation functions between the electrons in the dangling $sp^{2}$
orbitals and the surrounding $\pi$ electrons are insensitive to the carrier
doping, while those between the electrons in the local $\pi$ orbitals around
the vacancy and the surrounding $\pi$ electrons behave differently in the
undoped and doped cases. The latter is attributed to the nature of the
hybridization function that shows a pseudo-gap structure in the undoped case,
thus leading to the unscreening of the local magnetic moment, but does not in
the doped case.

The rest of this paper is organized as follows. We first introduce an
effective Anderson model for a single vacancy in graphene in Sec. II. We then
summarize, as the building blocks of the Anderson model, the low-lying
multiplet structures of the local impurity part of the Anderson model,
corresponding to an effective three-site two-orbital Hubbard model, and the
hybridization function that describes the energy-resolved coupling between the
local multiplets and the surrounding $\pi$ orbitals forming the conduction
band of the Anderson model. The numerical method used in this study is briefly
explained in Sec. III and the numerical results are shown in Sec. IV. The
paper is concluded in Sec. V by discussing the relevance of our results to the
experiments and also making the comparison with the previous theoretical
studies. The explicit form of the local multiplet states in the strong
coupling limit is provided in Appendix A. The spin correlation functions for
similar but much simpler models are calculated as the reference for the
comparison in Appendix B.

##  II Model

###  II.1 Anderson model for a single vacancy in graphene

Graphene is composed of carbon atoms forming a two-dimensional honeycomb
lattice structure with a single carbon atom locating at each site. The
$sp^{2}$ orbitals of a carbon atom in graphene are three local orbitals formed
by the linear combination of the 2$s$ atomic orbital and two 2$p$ atomic
orbitals (2$p_{x}$ and 2$p_{y}$ orbitals) in the plane. Each of these $sp^{2}$
orbitals has a large lobe pointing to one of three different directions and
the angle between any two of them is $120^{\circ}$, hence compatible to the
honeycomb lattice structure of graphene. In graphene, two of these $sp^{2}$
orbitals in neighboring sites point to each other (one from each site) and
hybridize strongly (i.e., $\sigma$ bonding) to form bonding and antibonding
orbitals between the neighboring sites. Because of the large energy level
difference between the bonding and antibonding orbitals, the bonding orbital
is fully occupied by two electrons and the antibonding orbital is completely
empty. Therefore, these orbitals formed by the $sp^{2}$ orbitals are usually
inactive in the pristine graphene. The electronic properties of graphene are
thereby determined by the remaining $\pi$ orbitals, i.e., $2p_{z}$ orbitals
perpendicular to the graphene plane, and the resulting energy band exhibits
the Dirac linear dispersion across the Fermi level Novoselov  _et al._ (2005).
The $sp^{2}$ orbitals are referred to simply as $\sigma$ orbitals hereafter.

These features can be depicted schematically in Fig. 1(a) with orange (blue)
spheres representing $\pi$ ($\sigma$) orbitals. Each carbon atom has six
electrons. Two of them occupy the 1s core orbital and the remaining four
electrons occupy the $\pi$ and $\sigma$ orbitals with one electron per
orbital. All of the $\sigma$ orbitals have their neighboring pairs to form the
covalent $\sigma$ bonding.

Now, let us consider to introduce a single vacancy of the carbon atom in
graphene. As shown in Fig. 1(b), in the presence of a single vacancy, there
appear three $\sigma$ orbitals around the vacancy, which have no neighboring
pairs. These unpaired $\sigma$ orbitals are called dangling orbitals and are
tightly localized around the vacancy. On the other hand, the removal of one
site generates the imbalance between sublattices in the honeycomb lattice Lieb
(1989), which induces a zero energy mode in the $\pi$ orbital system Pereira
_et al._ (2006); Peres  _et al._ (2009); Ugeda  _et al._ (2010). The spectral
weight (or equivalently the wave function) of the zero energy mode is mostly
found around the vacancy, although it should be more extended in space than
the dangling $\sigma$ orbitals.

Figure 1:  (a) Schematic depiction of the pristine graphene composed of $\pi$
orbitals (orange spheres) and $\sigma$ orbitals (blue spheres). Each carbon
site has one $\pi$ orbital and three $\sigma$ orbitals. Here, $\sigma$
orbitals refer to the $sp^{2}$ orbitals and their directional nature of the
shape of orbitals is indicated by locating them slightly away from the lattice
point toward the neighboring sites. (b) Schematic depiction of graphene in the
presence of a single vacancy. Three dangling $\sigma$ orbitals appear around
the vacancy. (c) Schematic depiction of the effective Anderson model
$\mathcal{H}_{\rm AM}$ studied here. The model consists of the three $\sigma$
orbitals around the vacancy and all $\pi$ orbitals of carbon atoms forming the
honeycomb lattice structure with a single vacancy. Three sites around the
vacancy are treated as an “impurity site” in the Anderson model (inside a blue
dashed circle denoted as $\mathbb{I}$) and the remaining sites are treated as
the conduction sites (outside a blue dashed circle denoted as $\mathbb{L}$).
Green solid (dashed) lines in (b) and (c) indicate the bonds inside region
$\mathbb{I}$ with a finite transfer integral $t_{\sigma}$ ($t_{\pi}$) between
the $\sigma$ orbitals ($\pi$ orbitals).

For the modeling of this system, i.e., a single vacancy in graphene, we
introduce the following simplifications. First, as shown in Fig. 1(c), we
neglect all paired $\sigma$ orbitals because they are inactive. Next, we
assume that the conduction electrons in the surrounding $\pi$ orbitals away
from the vacancy are described by a noninteracting tight-binding model,
implicitly considering the renormalization effect of the correlated two-
dimensional massless Dirac electrons Otsuka  _et al._ (2016); Seki  _et al._
(2019). Furthermore, we ignore possible modification of the hopping matrix
elements around the vacancy. Within these assumptions, one can model our
system as the following effective Anderson model:

$\mathcal{H}_{\rm AM}=\mathcal{H}_{\mathbb{I}}+\mathcal{H}_{\mathbb{L}},$ (1)

where $\mathcal{H}_{\mathbb{I}}$ and $\mathcal{H}_{\mathbb{L}}$ describe the
local part around the vacancy and the surrounding part of the system away from
the vacancy, respectively [see Fig. 1(c)].

The local part of the Hamiltonian, $\mathcal{H}_{\mathbb{I}}$, is composed of
the three dangling $\sigma$ orbitals and the three $\pi$ orbitals on the three
carbon sites next to the vacancy:

$\mathcal{H}_{\mathbb{I}}=\mathcal{H}_{t}+\mathcal{H}_{U}+\mathcal{H}_{d},$
(2)

where

$\displaystyle\mathcal{H}_{t}$
$\displaystyle=\varepsilon_{\sigma}\sum_{i\in\mathbb{I}}\sum_{s=\uparrow,\downarrow}c_{i\sigma
s}^{\dagger}c_{i\sigma
s}+\varepsilon_{\pi}\sum_{i\in\mathbb{I}}\sum_{s=\uparrow,\downarrow}c_{i\pi
s}^{\dagger}c_{i\pi s}$ $\displaystyle-
t_{\sigma}\sum_{\scriptsize\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}i,j\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\mathbb{I}}}\sum_{s=\uparrow,\downarrow}c_{i\sigma
s}^{\dagger}c_{j\sigma
s}-t_{\pi}\sum_{\scriptsize\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}i,j\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\mathbb{I}}}\sum_{s=\uparrow,\downarrow}c_{i\pi
s}^{\dagger}c_{j\pi s},$ (3) $\displaystyle\mathcal{H}_{U}$
$\displaystyle=U\sum_{i\in\mathbb{I}}\sum_{\alpha=\sigma,\pi}n_{i\alpha\uparrow}n_{i\alpha\downarrow}+U^{\prime}\sum_{i\in\mathbb{I}}n_{i\sigma}n_{i\pi}$
$\displaystyle-2J\sum_{i\in\mathbb{I}}\left({\bm{S}}_{i\sigma}\cdot{\bm{S}}_{i\pi}+\frac{1}{4}n_{i\sigma}n_{i\pi}\right)$
$\displaystyle+J^{\prime}\sum_{i\in\mathbb{I}}\left(c_{i\sigma\uparrow}^{\dagger}c_{i\sigma\downarrow}^{\dagger}c_{i\pi\downarrow}c_{i\pi\uparrow}+c_{i\pi\uparrow}^{\dagger}c_{i\pi\downarrow}^{\dagger}c_{i\sigma\downarrow}c_{i\sigma\uparrow}\right),$
(4)

and

$\mathcal{H}_{d}=-\frac{1}{2}(U+2U^{\prime}-J)\sum_{i\in\mathbb{I}}\sum_{\alpha=\sigma,\pi}n_{i\alpha}.$
(5)

Here, $\mathbb{I}$ represents the set of site labels $i$ for the three carbon
sites next to the vacancy [see Fig. 1(c)],
$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}i,j\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\mathbb{I}}$
denotes a pair of sites $i$ and $j$ in $\mathbb{I}$, corresponding to the next
nearest neighbors in the honeycomb lattice. $c_{i\alpha s}$ ($c_{i\alpha
s}^{\dagger}$) is the annihilation (creation) operator of an electron at site
$i$ with spin $s$ ($=\uparrow,\downarrow$) and orbital $\alpha$ (=$\sigma$,
$\pi$), $n_{i\alpha s}=c_{i\alpha s}^{\dagger}c_{i\alpha s}$, and
$n_{i\alpha}=n_{i\alpha\uparrow}+n_{i\alpha\downarrow}$.
${\bm{S}}_{i\alpha}=(S_{i\alpha}^{x},S_{i\alpha}^{y},S_{i\alpha}^{z})$ is the
local spin operator for orbital $\alpha$ at site $i$ given by

$S_{i\alpha}^{\nu}=\frac{1}{2}\sum_{s,s^{\prime}=\uparrow,\downarrow}c_{i\alpha
s}^{\dagger}[\mbox{\boldmath{$\sigma$}}^{\nu}]_{ss^{\prime}}c_{i\alpha
s^{\prime}},\quad(\nu=x,y,z)$ (6)

where $\mbox{\boldmath{$\sigma$}}^{x}$, $\mbox{\boldmath{$\sigma$}}^{y}$, and
$\mbox{\boldmath{$\sigma$}}^{z}$ represent the Pauli matrices. We also define
the total spin operator at site $i$,

${\bm{S}}_{i}={\bm{S}}_{i\sigma}+{\bm{S}}_{i\pi}$ (7)

with ${\bm{S}}_{i}=(S_{i}^{x},S_{i}^{y},S_{i}^{z})$, and the total spin
operator over all sites in $\mathbb{I}$,

${\bm{S}}_{\mathbb{I}}=\sum_{i\in\mathbb{I}}{\bm{S}}_{i}$ (8)

with
${\bm{S}}_{\mathbb{I}}=(S_{\mathbb{I}}^{x},S_{\mathbb{I}}^{y},S_{\mathbb{I}}^{z})$.

The first two terms in $\mathcal{H}_{t}$ represent the local energy levels of
the $\sigma$ and $\pi$ orbitals. The last two terms in $\mathcal{H}_{t}$ are
the hopping terms between different sites. Note that there is no hopping
between the two different orbitals because the $\sigma$ orbital and the $\pi$
orbital are orthogonal to each other. For simplicity, we parametrize
$t_{\sigma}=2t_{\pi}$ throughout this paper. $\mathcal{H}_{U}$ represents the
local interaction terms, including the intra-orbital Coulomb interaction $U$,
the inter-orbital Coulomb interaction $U^{\prime}$, the Hund’s coupling $J$,
and the pair hopping $J^{\prime}$. To reduce the number of parameters, here we
assume $U=U^{\prime}+2J$, $J=J^{\prime}$, and $J=0.1U$, although the former is
valid only for atomic orbitals Kanamori (1963). Notice also that
$\mathcal{H}_{d}$ is added to correct the double counting of the interactions.
While $\mathcal{H}_{d}$ does not make any difference when the isolated
$\mathcal{H}_{\mathbb{I}}$ is considered, it becomes important to compensate
the electron density when the surrounding $\pi$ orbitals described by
$\mathcal{H}_{\mathbb{L}}$ is incooporated.

In the Anderson model $\mathcal{H}_{\rm AM}$ in Eq. (1), the effects of the
surrounding $\pi$ orbitals are encapsulated by the noninteracting Hamiltonian,

$\displaystyle\mathcal{H}_{\mathbb{L}}=$ $\displaystyle-t\sum_{\langle
i,j\rangle^{\prime}}\sum_{s=\uparrow,\downarrow}(c_{i\pi s}^{\dagger}c_{j\pi
s}+c_{j\pi s}^{\dagger}c_{i\pi s})$ $\displaystyle-
t_{\pi}\sum_{\mbox{\scriptsize{$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}$}}i,j\mbox{\scriptsize{$\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}$}}^{\prime}}\sum_{s=\uparrow,\downarrow}(c_{i\pi
s}^{\dagger}c_{j\pi s}+c_{j\pi s}^{\dagger}c_{i\pi s}),$ (9)

where $\langle i,j\rangle^{\prime}$ and
$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}i,j\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}^{\prime}$
indicate pairs of nearest-neighbor and next-nearest-neighbor sites $i$ and
$j$, respectively, on the honeycomb lattice except for pairs of sites inside
$\mathbb{I}$, implying that pairs of sites across the boundary between
$\mathbb{I}$ and $\mathbb{L}$ are also included. As indicated in Fig. 1(c), in
the follows, we simply use $\mathbb{L}$ to denote the set of all sites in the
honeycomb lattice but excluding sites $i\in\mathbb{I}$, i.e., all sites in the
surrounding $\pi$ orbital system.

###  II.2 Local multiplets

Setting aside the surrounding $\pi$ orbital system, we shall first examine the
low-lying multiplet structure of the local part of the Hamiltonian
$\mathcal{H}_{\mathbb{I}}$ given in Eq. (2), which is essentially the three-
site two-orbital Hubbard model. The local part of the Hamiltonian
$\mathcal{H}_{\mathbb{I}}$ can be considered as the building block, i.e.,
“impurity site”, of the Anderson model $\mathcal{H}_{\rm AM}$ and the
understanding of the low-lying multiplet structure is essential to understand
the electronic structure of the total system. Since each orbital has nominally
one electron in the effective model of graphene considered here, we impose
that the total number of electrons in $\mathcal{H}_{\mathbb{I}}$ is six.

The model described by $\mathcal{H}_{\mathbb{I}}$ in Eq. (2) possesses the
$\mathcal{C}_{3}$ point group symmetry as well as the spin rotational
symmetry. Therefore, any eigenstate of $\mathcal{H}_{\mathbb{I}}$ can be
simultaneously characterized with the irreducible representation of the point
group and the total spin. As already described in Sec. I, the irreducible
representations of the $\mathcal{C}_{3}$ point group are the symmetric one-
dimensional representation $\mathcal{A}$ and the doubly degenerate
representation $\mathcal{E}$. In this paper, we denote sets of states
belonging to these irreducible representations $\mathcal{A}$ and $\mathcal{E}$
as $\mathbb{A}$ and $\mathbb{E}$, respectively. We also represent the quantum
numbers of ${\bm{S}}_{\mathbb{I}}$ and $S_{\mathbb{I}}^{z}$ as
$S_{\mathbb{I}}$ and $M_{\mathbb{I}}$, respectively. Namely, the eigenvalue of
${\bm{S}}_{\mathbb{I}}\cdot{\bm{S}}_{\mathbb{I}}$ is
$S_{\mathbb{I}}(S_{\mathbb{I}}+1)$. Using these notations, we can refer to the
energy eigenstates with ($\mathbb{C}$, $S_{\mathbb{I}}$)
($\mathbb{C}=\mathbb{A}$, $\mathbb{E}$) as the multiplets characterized by the
irreducible representation $\mathcal{C}$ ($\mathcal{C}=\mathcal{A}$,
$\mathcal{E}$) and the total spin $S_{\mathbb{I}}$.

Figure 2:  (a, b) Single-particle energy level diagrams and the ground-state
electron configurations of $\mathcal{H}_{\mathbb{I}}$ in the noninteracting
limit for (a)
$-t_{\sigma}<\varepsilon_{\pi}-\varepsilon_{\sigma}<t_{\sigma}/2$ and (b)
$t_{\sigma}/2<\varepsilon_{\pi}-\varepsilon_{\sigma}<2t_{\sigma}$ with
$t_{\sigma}=2t_{\pi}>0$. $(a,\alpha)$ [($e,\alpha$)] for orbital
$\alpha\,(=\sigma,\pi)$ characterizes each single-particle energy level with
the irreducible representation $\mathcal{A}$ ($\mathcal{E}$) of the
$\mathcal{C}_{3}$ point group. Green arrows represent up and down electrons.
(c, d) Low-lying multiplet energy diagrams in (c) the weak coupling limit
($U\ll t_{\sigma}$) and (d) the strong coupling limit ($U\gg t_{\sigma}$).
$(\mathbb{C},S_{\mathbb{I}})$ ($\mathbb{C}=\mathbb{A},\mathbb{E}$,
$S_{\mathbb{I}}=0,1,2,3$) characterizes the associated many-body energy
eigenstates, i.e., multiplets, with the irreducible representation
$\mathcal{X}$ ($\mathcal{X}=\mathcal{A}$, $\mathcal{E}$) and the total spin
$S_{\mathbb{I}}$. Here the total number of electrons is assumed to be 6.

In the noninteracting case, the single-particle energy level
$\varepsilon_{\alpha}$ for each orbital $\alpha\,(=\sigma,\pi)$ is split by
the electron hopping $t_{\alpha}$ into a nondegenerate energy level with
$\mathbb{A}$ and doubly degenerate energy levels with $\mathbb{E}$. We denote
these single-particle energy levels as $a\in\mathbb{A}$ and $e\in\mathbb{E}$
with their energy levels $\varepsilon_{\alpha}^{(a)}$ and
$\varepsilon_{\alpha}^{(e)}$, respectively. These energy levels
$\varepsilon_{\alpha}^{(a)}$ and $\varepsilon_{\alpha}^{(e)}$ are easily
evaluated as $\varepsilon_{\sigma}^{(a)}=-2t_{\sigma}+\varepsilon_{\sigma}$
and $\varepsilon_{\sigma}^{(e)}=t_{\sigma}+\varepsilon_{\sigma}$ for the
$\sigma$ orbitals, and $\varepsilon_{\pi}^{(a)}=-2t_{\pi}+\varepsilon_{\pi}$
and $\varepsilon_{\pi}^{(e)}=t_{\pi}+\varepsilon_{\pi}$ for the $\pi$
orbitals. For $t_{\sigma}=2t_{\pi}>0$ and
$-t_{\sigma}<\varepsilon_{\pi}-\varepsilon_{\sigma}<t_{\sigma}/2$, the ground-
state electron configuration with six electrons is open shell with two
electrons occupying the $e$ level of $\pi$ orbitals, as shown in Fig. 2(a).
Notice that in this case the number of electrons in the $\sigma$ ($\pi$)
orbitals, $N_{\sigma}$ ($N_{\pi}$), is $N_{\sigma}=2$ ($N_{\pi}=4$). The first
order perturbation theory with respect to the Coulomb interactions lifts the
six-fold degenerate electron configurations of the open shell $e$ level of
$\pi$ orbitals. Consequently, the ground state becomes ($\mathbb{A},1$), as
shown in Fig. 2(c). When
$t_{\sigma}/2<\varepsilon_{\pi}-\varepsilon_{\sigma}<2t_{\sigma}$ is realized
together with $t_{\sigma}=2t_{\pi}>0$, the $e$ level of $\sigma$ orbitals is
now lower than the $e$ level of $\pi$ orbitals, thus indicating that
$N_{\sigma}=4$ and $N_{\pi}=2$, as shown in Fig. 2(b). However, the resulting
ground state in this case is again ($\mathbb{A},1$) in the weak coupling
limit, as shown in Fig. 2(c).

In the strong coupling limit, the Coulomb interactions force the electron
distribution to be uniform, yielding $N_{\sigma}=N_{\pi}=3$. Moreover, the
Hund’s coupling $J$ makes the local state a spin triplet ($S_{i}=1$) at each
site. Since the kinetic exchange mechanism induces the antiferromagnetic
interactions between neighboring sites, the low-energy effective model for
$\mathcal{H}_{\mathbb{I}}$ in the strong coupling limit is represented by the
three-site spin-1 antiferromagnetic Heisenberg model. The low-lying energy
levels for this effective model are summarized in Fig. 2(d). The ground state
is ($\mathbb{A},0$), which refers to a kind of the valence-bond-type state
Affleck  _et al._ (1987) in the spin-1 antiferromagnetic chain Haldane (1983).
More detailed analysis of the strong coupling limit is provided in Appendix A.

In order to confirm the multiplet energy diagrams shown in Figs. 2(c) and
2(d), we show in Fig. 3 the low-lying energy $\Delta E_{n}$ defined as

$\Delta E_{n}=E_{n}-E_{0},$ (10)

where $E_{n}$ ($n=0,1,2,\cdots$) is the $n$th energy eigenvalues of
$\mathcal{H}_{\mathbb{I}}$ with six electrons and $n=0$ corresponds to the
ground state energy, i.e., $E_{0}\leq E_{1}\leq E_{2}\leq\cdots$, calculated
numerically by fully diagonalizing the Hamiltonian $\mathcal{H}_{\mathbb{I}}$.
Characterizing the energy eigenstates, we find that the low-lying multiplet
structure is precisely reproduced by those shown in Fig. 2(c) in the weak
coupling limit ($U\approx 0$) and in Fig. 2(d) in the strong coupling limit
($U\gg t_{\sigma}$), except for the small energy splitting between the
multiplets with $(\mathbb{A},1)$ and $(\mathbb{E},1)$, which is expected to be
degenerate in the limit of $U/t_{\sigma}\to\infty$. Notice also in Fig. 3 that
the transition from the weak coupling limit to the strong coupling limit,
occurring at $U=3t_{\sigma}$ for the specific parameter used here, is
essentially a level crossing. Indeed, these two states with $(\mathbb{A},1)$
and $(\mathbb{A},0)$ are distinguishable with different eigenvalues of the
total spin $S_{\mathbb{I}}$.

Figure 3:  Low-lying energy $\Delta E_{n}=E_{n}-E_{0}$ of
$\mathcal{H}_{\mathbb{I}}$ as a function of $U$. Here, $E_{n}$ is the $n$th
energy eigenvalues of $\mathcal{H}_{\mathbb{I}}$ with $n=0$ corresponding to
the ground state and $E_{0}\leq E_{1}\leq E_{2}\leq\cdots$.
$(\mathbb{C},S_{\mathbb{I}})$ ($\mathbb{C}=\mathbb{A},\mathbb{E}$,
$S_{\mathbb{I}}=0,1,2,3$) is indicated for some of the energy eigenvalues
$E_{n}$ to characterize the corresponding energy eigenstates. The parameters
used here are $t_{\sigma}=2t_{\pi}>0$,
$\varepsilon_{\sigma}=\varepsilon_{\pi}$, and $J/U=0.1$. The ground state
($n=0$) is $(\mathbb{A},1)$ for $U/t_{\sigma}<3$ and $(\mathbb{A},0)$ for
$U/t_{\sigma}>3$. The total number of electrons is assumed to be 6.

Figure 4(a) shows the intensity plot of
$\langle{\bm{S}}_{\mathbb{I}}\cdot{\bm{S}}_{\mathbb{I}}\rangle\equiv\bar{S}_{\mathbb{I}}(\bar{S}_{\mathbb{I}}+1)$
obtained by numerically diagonalizing $\mathcal{H}_{\mathbb{I}}$ with six
electrons. Here, $\langle\mathcal{O}\rangle$ indicates the expectation value
of an operator $\mathcal{O}$ for the ground state. As expected, the total spin
of the ground state for the small $U/t_{\sigma}$ region is
$\bar{S}_{\mathbb{I}}=1$ and that for the large $U/t_{\sigma}$ region is
$\bar{S}_{\mathbb{I}}=0$. We also find in Fig. 4(b) that, for the large
$U/t_{\sigma}$ region, the total spin at each site, $\bar{S}_{i}$, evaluated
from
$\langle{\bm{S}}_{i}\cdot{\bm{S}}_{i}\rangle\equiv\bar{S}_{i}(\bar{S}_{i}+1)$,
approaches to $\bar{S}_{i}=1$ with increasing $U/t_{\sigma}$, suggesting the
formation of a local spin 1 at each site, where the low-lying states can be
described by the spin-1 antiferromagnetic Heisenberg model. It is also clear
in Fig. 4(a) that the phase boundary between the weak and strong coupling
phases corresponds to a level crossing because the total spin of the ground
state is different.

Figure 4:  Intensity plots of (a)
$\langle{\bm{S}}_{\mathbb{I}}\cdot{\bm{S}}_{\mathbb{I}}\rangle\equiv\bar{S}_{\mathbb{I}}(\bar{S}_{\mathbb{I}}+1)$
and (b)
$\langle{\bm{S}}_{i}\cdot{\bm{S}}_{i}\rangle\equiv\bar{S}_{i}(\bar{S}_{i}+1)$
calculated for the ground state of $\mathcal{H}_{\mathbb{I}}$ with
$t_{\sigma}=2t_{\pi}>0$, $\varepsilon_{\sigma}=\varepsilon_{\pi}$, and
$J/U=0.1$. The values of $\bar{S}_{\mathbb{I}}$ and $\bar{S}_{i}$ are also
indicated in (a) and (b), respectively. The total number of electrons is
assumed to be 6.

There are two remarks in order. First, any multiplet state in both weak and
strong coupling phases cannot be described by a single Slater determinant,
indicating the importance of the correlation effect. Especially, the Hund’s
coupling play an essential role on realizing the valence-bond-type state in
the strong coupling phase. Second, $t_{\pi}$ is the next nearest neighbor
hopping for the $\pi$ orbitals and is estimated as large as $t_{\pi}\sim
t/12$, where $t$ is the nearest neighbor hopping for the $\pi$ orbitals in
graphene Kretinin  _et al._ (2013). On the other hand, the intra-orbital
Coulomb interaction for the $\pi$ orbitals is estimated as large as $\sim
3.5t$ by the constrained random phase approximation Wehling  _et al._ (2011).
Therefore, it is reasonable to expect that the ground state of
$\mathcal{H}_{\mathbb{I}}$ with the parameter set for graphene is in the
strong coupling phase. However, as we will discuss below, the effect of the
surrounding $\pi$ orbitals is dominant to determine the low-lying electronic
structure of the Anderson model $\mathcal{H}_{\rm AM}$.

###  II.3 Hybridization function

It is known that an Anderson model in general is characterized by a local
Hamiltonian for the impurity sites, which corresponds to
$\mathcal{H}_{\mathbb{I}}$ in our case, and a so-called hybridization function
$\Gamma_{ij}(\omega)$ Bulla  _et al._ (1997), mathematically defined as the
Schur complement of the inverse of the hopping matrix (i.e., Green’s function)
for the complementary space of the impurity sites in $\mathbb{I}$, which
corresponds to $\mathbb{L}$ in our case Shirakawa and Yunoki (2014).
Physically, $\Gamma_{ij}(\omega)$ is a quantity characterizing the energy
resolved coupling of the impurity sites to the surrounding bath sites, i.e.,
the surrounding $\pi$ orbital system in $\mathbb{L}$. For the Anderson model
described by the Hamiltonian $\mathcal{H}_{\rm AM}$ given in Eq. (1),
$\Gamma_{ij}(\omega)$ is a three-by-three matrix labeled by the site indices
of the $\pi$ orbitals in $\mathbb{I}$, since the $\sigma$ orbitals in
$\mathbb{I}$ do not hybridize with the surrounding $\pi$ orbital system in
$\mathbb{L}$. Because of the $\mathcal{C}_{3}$ point group symmetry of the
system, $\Gamma_{ij}(\omega)$ can be diagonalized with the diagonal elements
$\tilde{\Gamma}_{a}(\omega)$ and $\tilde{\Gamma}_{e}(\omega)$ corresponding to
the couplings with the $(a,\pi)$ orbital and the $(e,\pi)$ orbitals
($e=e_{1},e_{2}$) in $\mathbb{I}$, respectively. We call the $(a,\pi)$-orbital
system $a$-mode and the $(e,\pi)$-orbital systems $e$-modes.

Figure 5: Hybridization functions $\tilde{\Gamma}_{c}(\omega)$ ($c=a,e$) for
$\mathcal{H}_{\rm AM}$ with $t_{\pi}=t/12$. Green dashed line indicates the
Fermi level at $\omega=0$ for bulk graphene, i.e., undoped graphene. Note that
hybridization functions $\tilde{\Gamma}_{e}(\omega)$ for $e=e_{1}$ and $e_{2}$
are exactly the same and here only one of them is plotted.

Figure 5 shows $\tilde{\Gamma}_{a}(\omega)$ and $\tilde{\Gamma}_{e}(\omega)$
for the Anderson model $\mathcal{H}_{\rm AM}$ given in Eq. (1) with
$t_{\pi}=t/12$ Kretinin  _et al._ (2013). Around the Fermi level for bulk
graphene, i.e., undoped graphene, a diverging behavior appears in
$\tilde{\Gamma}_{e}(\omega)$, which is associated with the zero energy modes
of the $\pi$ orbital system Pereira  _et al._ (2006); Peres  _et al._ (2009),
as described below. This diverging behavior strongly affects the distribution
of the local multiplets coupled to the $e$-modes, as shown later. We should
note here that because of the particle-hole asymmetry due to the presence of
the next-nearest hopping $t_{\pi}$, $\tilde{\Gamma}_{e}(\omega)$ does not
exactly diverge and the maximum location in energy is slightly shifted below
the Fermi level Pereira  _et al._ (2006). In this case, there are two phases
expected for the ground state, according to the previous studies Mitchell and
Fritz (2013). One is the symmetric strong coupling phase where the Kondo
screening state arises. The other is the asymmetric strong coupling phase
where the spin degree of freedom is vanished because of the large potential
difference between the impurity and conduction sites that induces the closed-
shell configuration at the impurity site. The divergent hybridization function
is preferable to the symmetric strong coupling even when the model is
particle-hole asymmetric.

On the other hand, we find in Fig. 5 that $\tilde{\Gamma}_{a}(\omega)$ shows a
pseudogap structure at the Fermi level. An Anderson model with a single
magnetic impurity coupled to the bath sites through a hybridization function
with a pseudogap structure at the Fermi level is called a pseudogap Anderson
model Fritz and Vojta (2013). The pseudogap Anderson model is indeed
considered as an effective model for a hydrogen adatom on graphene Shirakawa
and Yunoki (2014, 2016). Previous studies have revealed that the Kondo
screening is absent in the pseudogap Anderson model Gonzalez-Buxton and
Ingersent (1998); Vojta  _et al._ (2010). Consequently, there appears the free
magnetic moment even at zero temperature. In other words, the hybridization
function with the pseudogap structure is irrelevant to the low-lying
electronic properties of the impurity site. The carrier doping into the
pseudogap Anderson model yields either the Kondo screening of the magnetic
moment at the impurity site or the vanishment of the magnetic moment by adding
or removing electrons at the impurity site Vojta  _et al._ (2010), suggesting
that the free magnetic moment in the adatom graphene would be sensitive
against the carrier doping.

Finally, we note that the diverging part in $\tilde{\Gamma}_{e}(\omega)$ is
sometimes referred to as the zero energy modes Pereira  _et al._ (2006); Peres
_et al._ (2009); Ugeda  _et al._ (2010). The number of the zero energy modes
Pereira  _et al._ (2006); Peres  _et al._ (2009); Ugeda  _et al._ (2010) is
related to the imbalance of the number of sites in the two sublattices of a
bipartite system Lieb (1989). If we consider the system $\mathbb{L}$, the
imbalance of the number of sites in the two sublattices is $2$ ($=3-1$), which
coincides with the degeneracy of $\mathcal{E}$. From this rule, we can readily
understand that $\tilde{\Gamma}_{e}(\omega)$ should show the diverging
behavior because of the presence of the two zero energy modes, while
$\tilde{\Gamma}_{a}(\omega)$ should not exhibit any diverging behavior that
originates from the zero energy modes caused by the sublattice imbalance.

##  III Method

We use the density-matrix renormalization group (DMRG) method White (1992,
1993) to study the ground state properties of the Anderson model described by
the Hamiltonian $\mathcal{H}_{\rm AM}$ in Eq. (1). For this purpose, we employ
the block Lanczos tridiagonalization technique to map the Anderson model
$\mathcal{H}_{\rm AM}$ onto an effective quasi-one-dimensional (Q1D) model
Shirakawa and Yunoki (2014); Allerdt  _et al._ (2015, 2017); Allerdt and
Feiguin (2019), which consists of the three $\sigma$ orbitals and the three
$\pi$ orbitals in $\mathbb{I}$, describing the local part of the Hamiltonian
$\mathcal{H}_{\mathbb{I}}$, and $3(L-1)$ $\pi$ orbitals in $\mathbb{L}$,
generated by the block Lanczos procedure starting from the three $\pi$
orbitals in $\mathbb{I}$, to represent the surrounding $\pi$ orbital system
that couples to the $\pi$ orbitals in $\mathbb{I}$ [see Figs. 6(a) and 6(b)].
Here, $L$ is the maximum number of the block Lanczos iterations. With the
local part of the Hamiltonian $\mathcal{H}_{\mathbb{I}}$ in Eqs. (2)–(5), the
Hamiltonian $\mathcal{H}_{\rm Q1D}$ of the resulting effective Q1D model is
thus given by

$\displaystyle\mathcal{H}_{\rm Q1D}$
$\displaystyle=\mathcal{H}_{\mathbb{I}}+\sum_{s=\uparrow,\downarrow}\mathcal{H}^{s}_{c},$
(11)

where

$\displaystyle\mathcal{H}^{s}_{c}=$
$\displaystyle\sum_{l=1}^{L}\sum_{m=1}^{3}\sum_{m^{\prime}=1}^{3}[\mbox{\boldmath{$\varepsilon$}}_{l}]_{m,m^{\prime}}\tilde{c}_{lms}^{\dagger}\tilde{c}_{lm^{\prime}s}$
$\displaystyle+\sum_{l=1}^{L-1}\sum_{m=1}^{3}\sum_{m^{\prime}=1}^{3}[\mbox{\boldmath{$\tau$}}_{l}]_{m,m^{\prime}}\tilde{c}_{lms}^{\dagger}\tilde{c}_{l+1m^{\prime}s}+{\rm
h.c.}$ (12)

and $\tilde{c}_{lms}$ ($m=1,2,3$, $l=1,2,\cdots,L$) is an electron
annihilation operator defined by the linear combination of $c_{i\pi s}$ in
$\mathcal{H}_{\rm AM}$ through

$\tilde{c}_{lms}=\sum_{i}U^{l}_{i,m}c_{i\pi s}.$ (13)

Note that $\tilde{c}_{lms}$ ($m=1,2,3$) for $l=1$ is identical to $c_{i\pi s}$
for $i\in\mathbb{I}$ and therefore the form of $\mathcal{H}_{\mathbb{I}}$,
including the interaction part $\mathcal{H}_{U}$, remains unchanged under this
transformation.

Figure 6:  Schematic depictions of (a) the Anderson model $\mathcal{H}_{\rm
AM}$ given in Eq. (1) for a single vacancy in graphene, (b) the effective Q1D
three-leg Anderson model $\mathcal{H}_{\rm Q1D}$ obtained by the block Lanczos
transformation with the $\pi$ orbitals in $\mathbb{I}$ as the initial block
Lanczos bases, and (c) the symmetrized Q1D three-leg Anderson model
$\mathcal{H}_{\rm Q1D}^{\prime}$ obtained by the unitary transformation in
Eqs. (21)–(23). $m$ ($=1,2,3$) and $c$ ($=a,e_{1},e_{2}$) stand for labels of
legs in the Q1D ladder representations, and $l$ stands for a label of rungs
corresponding to the block Lanczos basis generated by the block Lanczos
iteration. Blue spheres in (a)–(c) denote the $\sigma$ orbitals in
$\mathbb{I}$ represented by $(c_{i\sigma s},c_{i\sigma s}^{\dagger})$. Orange
and red spheres in (a), (b), and (c) denote the $\pi$ orbitals in both
$\mathbb{I}$ and $\mathbb{L}$ represented by $(c_{i\pi s},c_{i\pi
s}^{\dagger})$, $(\tilde{c}_{lms},\tilde{c}_{lms}^{\dagger})$, and
$(\tilde{f}_{l\pi s},\tilde{f}_{l\pi s}^{\dagger})$, respectively. Note that
$l=1$ corresponds to the initial bases for the block Lanczos iteration and
thus $(\tilde{c}_{lms},\tilde{c}_{lms}^{\dagger})$ for $l=1$ in (b) are
exactly the same as $(c_{i\pi s},c_{i\pi s}^{\dagger})$ for $i\in\mathbb{I}$.
The region indicated by the blue dashed rectangle corresponds to the local
Hamiltonian $\mathcal{H}_{\mathbb{I}}$, which remains unchanged under the
transformation in (b). However, the form of $\mathcal{H}_{\mathbb{I}}$ is
changed in (c). The left (right) part of the sites in (b) and (c) is assigned
to represent orbitals with up spin (down spin).

Let us now introduce the matrix
$[{\bm{U}}]_{i,3l+m}\equiv[{\bm{U}}_{l}]_{i,m}\equiv U^{l}_{i,m}$ with
$m=1,2,3$, $l=1,2,\cdots,L$, and $i\in\mathbb{I}+\mathbb{L}$ in Eq. (13). To
obtain the form of the Hamiltonian $\mathcal{H}_{c}^{s}$ in Eq. (12), we
should notice that the block Lanczos transformation is applied to the
surrounding part of the Hamiltonian $\mathcal{H}_{\mathbb{L}}$ defined in Eq.
(9). The matrix ${\bm{U}}$ represents the block Lanczos basis that transforms
the Hamiltonian matrix of $\mathcal{H}_{\mathbb{L}}$ for the $\pi$ orbitals
with spin $s$, i.e.,

$[{\bm{H}}_{\mathbb{L}}^{s}]_{ij}=\left\\{\begin{array}[]{ll}-t&i,j\in\langle
i,j\rangle\\\
-t_{\pi}&i,j\in\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}i,j\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}\\\
0&\text{otherwise}\\\ \end{array}\right.$ (14)

into a tri-block-diagonal matrix form, i.e.,

${\bm{U}}^{\dagger}{\bm{H}}_{\mathbb{L}}^{s}{\bm{U}}=\left(\begin{array}[]{cccccc}\mbox{\boldmath{$\varepsilon$}}_{1}&\mbox{\boldmath{$\tau$}}_{1}&{\bm{0}}&{\bm{0}}&\cdots&{\bm{0}}\\\
\mbox{\boldmath{$\tau$}}_{1}^{\dagger}&\mbox{\boldmath{$\varepsilon$}}_{2}&\mbox{\boldmath{$\tau$}}_{2}&{\bm{0}}&\cdots&{\bm{0}}\\\
{\bm{0}}&\mbox{\boldmath{$\tau$}}_{2}^{\dagger}&\mbox{\boldmath{$\varepsilon$}}_{3}&\mbox{\boldmath{$\tau$}}_{3}&{}\hfil&\vdots\\\
{\bm{0}}&{\bm{0}}&\mbox{\boldmath{$\tau$}}_{3}^{\dagger}&\ddots&\ddots&{}\hfil\\\
\vdots&\vdots&{}\hfil&\ddots&\mbox{\boldmath{$\varepsilon$}}_{L-1}&\mbox{\boldmath{$\tau$}}_{L-1}\\\
{\bm{0}}&{\bm{0}}&\cdots&{}\hfil&\mbox{\boldmath{$\tau$}}_{L-1}^{\dagger}&\mbox{\boldmath{$\varepsilon$}}_{L}\\\
\end{array}\right),$ (15)

where $\mbox{\boldmath{$\varepsilon$}}_{l}$ and $\mbox{\boldmath{$\tau$}}_{l}$
are the $3\times 3$ matrices appearing in Eq. (12). We can obtain these
matrices recursively through the block Lanczos iteration

${\bm{U}}_{l+1}\mbox{\boldmath{$\tau$}}_{l+1}^{\dagger}={\bm{H}}_{\mathbb{L}}^{s}{\bm{U}}_{l}-{\bm{U}}_{l}\mbox{\boldmath{$\varepsilon$}}_{l}-{\bm{U}}_{l-1}\mbox{\boldmath{$\tau$}}_{l},$
(16)

where $\mbox{\boldmath{$\varepsilon$}}_{1}=0$, ${\bm{U}}_{0}=0$,
${\bm{U}}_{1}=(\bm{e}_{i_{1}},\bm{e}_{i_{2}},\bm{e}_{i_{3}})$ with
$\bm{e}_{i}$ being a column unit vector, i.e., $(\bm{e}_{i})_{k}=\delta_{i,k}$
and $\\{i_{1},i_{2},i_{3}\\}\equiv\mathbb{I}$, and the left hand side is
obtained by the QR decomposition of the right hand side. Note that the
particular form of $\bm{U}_{1}$ is due to the initial block Lanczos bases
chosen intentionally to be the three $\pi$ orbitals in $\mathbb{I}$. In
addition, as shown in Fig. 6(b), we can divide the system into two blocks
using the spin degrees of freedom to reduce the local Hilbert space Shirakawa
and Yunoki (2014). More details of the derivation and the technical
information for Eqs. (11)–(16) can be found in Ref. Shirakawa and Yunoki
(2014).

Because of the $\mathcal{C}_{3}$ rotational symmetry, we can find that the
matrix elements of $\mbox{\boldmath{$\varepsilon$}}_{l}$ and
$\mbox{\boldmath{$\tau$}}_{l}$ satisfy

$\displaystyle[\mbox{\boldmath{$\varepsilon$}}_{l}]_{11}=[\mbox{\boldmath{$\varepsilon$}}_{l}]_{22}=[\mbox{\boldmath{$\varepsilon$}}_{l}]_{33},$
(17)
$\displaystyle[\mbox{\boldmath{$\varepsilon$}}_{l}]_{12}=[\mbox{\boldmath{$\varepsilon$}}_{l}]_{23}=[\mbox{\boldmath{$\varepsilon$}}_{l}]_{31}=[\mbox{\boldmath{$\varepsilon$}}_{l}]_{21}=[\mbox{\boldmath{$\varepsilon$}}_{l}]_{32}=[\mbox{\boldmath{$\varepsilon$}}_{l}]_{13},$
(18)
$\displaystyle[\mbox{\boldmath{$\tau$}}_{l}]_{11}=[\mbox{\boldmath{$\tau$}}_{l}]_{22}=[\mbox{\boldmath{$\tau$}}_{l}]_{33},$
(19)

and

$\displaystyle[\mbox{\boldmath{$\tau$}}_{l}]_{12}=[\mbox{\boldmath{$\tau$}}_{l}]_{23}=[\mbox{\boldmath{$\tau$}}_{l}]_{31}=[\mbox{\boldmath{$\tau$}}_{l}]_{21}=[\mbox{\boldmath{$\tau$}}_{l}]_{32}=[\mbox{\boldmath{$\tau$}}_{l}]_{13}.$
(20)

Therefore, the unitary transformation of the electron operators

$\displaystyle\tilde{f}_{las}$
$\displaystyle=(\tilde{c}_{l1s}+\tilde{c}_{l2s}+\tilde{c}_{l3s})/\sqrt{3},$
(21) $\displaystyle\tilde{f}_{le_{1}s}$
$\displaystyle=(2\tilde{c}_{l1s}-\tilde{c}_{l2s}-\tilde{c}_{l3s})/\sqrt{6},$
(22)

and

$\displaystyle\tilde{f}_{le_{2}s}$
$\displaystyle=(\tilde{c}_{l2s}-\tilde{c}_{l3s})/\sqrt{2}$ (23)

for $l=1,2,\cdots,L$ transforms $\mathcal{H}_{\rm Q1D}$ into a decoupled form
(apart from the local part in $\mathbb{I}$) of the Q1D model

$\displaystyle\mathcal{H}_{\rm Q1D}^{\prime}$
$\displaystyle=\mathcal{H}_{\mathbb{I}}^{\prime}+\sum_{s=\uparrow,\downarrow}\mathcal{H}^{s}_{f},$
(24)

where

$\displaystyle\mathcal{H}^{s}_{f}=$
$\displaystyle\sum_{l=1}^{L}\sum_{c=a,e_{1},e_{2}}\tilde{\varepsilon}_{lc}\tilde{f}_{lcs}^{\dagger}\tilde{f}_{lcs}$
$\displaystyle+\sum_{l=1}^{L-1}\sum_{c=a,e_{1},e_{2}}\tilde{\tau}_{lc}\tilde{f}_{lcs}^{\dagger}\tilde{f}_{l+1cs}+{\rm
h.c.},$ (25)

as shown schematically in Fig. 6(c). Note that because of the unitary
transformation, the form of the local part of the Hamiltonian
$\mathcal{H}_{\mathbb{I}}$, including the interaction part, is changed to
$\mathcal{H}_{\mathbb{I}}^{\prime}$. We should also note that $c=a$
($c=e_{1},e_{2}$) indeed corresponds to the orbitals with the $\mathcal{A}$
($\mathcal{E}$) irreducible representation of the $\mathcal{C}_{3}$ point
group.

The systems sizes studied here are up to $L=47$, i.e., $6(L+1)=288$ sites in
the Q1D model including the spin degrees of freedom, which corresponds to
approximately $7000$ sites in the honeycomb lattice arround the vacancy
Shirakawa and Yunoki (2014). We keep the number $\chi$ of the density matrix
eigenstates up to $\chi=120(L+1)$ to obtain the reasonable convergence. For
example, we have checked that the results do not change much even if we
increase the number of the density matrix eigenstates kept up to
$\chi=144(L+1)$, and a typical error of the ground state energy is estimated
as large as $10^{-5}t$.

Finally, we must pay special attention to the treatment of the total number of
electrons in this Q1D representation. In this study, the number of electrons,
$N_{e}$, is determined so as to reproduce the chemical potential $\mu$ for the
bulk system. For a given chemical potential $\mu$, the ground potential
$\Omega(N,\mu)=E_{0}(N)-\mu N$ becomes minimum at $N_{e}$ if $N_{e}$ satisfies

$\frac{E_{0}(N_{e})-E_{0}(N_{e}-2)}{2}<\mu<\frac{E_{0}(N_{e}+2)-E_{0}(N_{e})}{2},$
(26)

where $E_{0}(N)$ is the ground state energy of the Q1D model with $N$
electrons. The condition (26) may also be written as

$N_{e}=\underset{N}{\mathrm{min}}\Omega(N,\mu)$ (27)

for a given $\mu$. To obtain a target chemical potential $\mu$ for a given set
of parameters ($\varepsilon_{\sigma}$, $\varepsilon_{\pi}$, $t_{\sigma}$,
$t_{\pi}$, $U$, $J$, $L$), we search $N_{e}$ that satisfies the condition in
Eq. (27) by calculating $E_{0}(N)$ for several values of $N$.

##  IV Results

In this section, we provide the numerical results obtained by the block
Lanczos DMRG method for the ground state of the Anderson model
$\mathcal{H}_{\rm AM}$ in Eq. (1). The results for undoped and doped cases are
shown in Sec. IV.1 and Sec. IV.2, respectively. The spin correlation functions
between the impurity and conduction sites are also examined in Sec. IV.3.
Here, we set the hopping parameters $t=6t_{\sigma}=12t_{\pi}$ Kretinin  _et
al._ (2013), although the realistic value for $t_{\sigma}$ would be larger
because the dangling $\sigma$ orbitals should be expanded around the vacancy.
We also set $\varepsilon_{\sigma}=\varepsilon_{\pi}=0$.

###  IV.1 Undoped case

As described at the end of the previous section, we first have to determine
the number of electrons, $N_{e}$, in the Q1D model $\mathcal{H}_{\rm Q1D}$ so
as to reproduce the chemical potential $\mu=2t_{\pi}$ of the undoped bulk
system. To this end, we calculate the ground state energy $E_{0}(N_{e})$ for a
given $L$ with varying the number $N_{e}$ of electrons, and find that the
condition $N_{e}:3(L+1)=7:6$ satisfies Eq. (27) with $\mu=2t_{\pi}$ for the
parameter region studied here. Note that $3(L+1)$ is the total number of sites
in the Q1D model $\mathcal{H}_{\rm Q1D}$, excluding the spin degrees of
freedom.

Figure 7 shows the $U$ dependence of the expectation values of local
operators, including the local density $n_{i\alpha}$ of orbital $\alpha$
($\alpha=\sigma,\pi$) at site $i\in\mathbb{I}$, the local spin squared
${\bm{S}}_{i}\cdot{\bm{S}}_{i}$ at site $i\in\mathbb{I}$, and the total spin
squared ${\bm{S}}_{\mathbb{I}}\cdot{\bm{S}}_{\mathbb{I}}$ over sites
$i\in\mathbb{I}$, for the ground state of the Anderson model $\mathcal{H}_{\rm
AM}$. We find that there are two distinguishable phases: $\langle
n_{i\sigma}\rangle=4/3$ in a weak coupling region ($U\leq t$) and $\langle
n_{i\sigma}\rangle=1$ in a strong coupling region ($U>t$). We refer to these
two phases as the weak coupling phase and the strong coupling phase,
respectively.

Figure 7:  Expectation values $\langle\mathcal{O}\rangle$ of operators
$\mathcal{O}=n_{i\sigma}$, $n_{i\pi}$, ${\bm{S}}_{i}\cdot{\bm{S}}_{i}$, and
${\bm{S}}_{\mathbb{I}}\cdot{\bm{S}}_{\mathbb{I}}$, where $i\in\mathbb{I}$,
calculated for the ground state of the Anderson model $\mathcal{H}_{\rm AM}$.
The parameters used here are $L=35$, $J=0.1U$, $t=6t_{\sigma}=12t_{\pi}$, and
$\mu=2t_{\pi}$, corresponding to the undoped case.

The abrupt change of the expectation values at $U\approx t$ found in Fig. 7
suggests that the transition from the weak to the strong coupling phase is of
the first order due to a level crossing. This is because of the symmetry
sector of the number of $\sigma$ electrons,
$N_{\sigma}=\sum_{i\in\mathbb{I}}n_{i\sigma}$. In our model, the symmetry
sectors with even $N_{\sigma}$ and odd $N_{\sigma}$ are orthogonal to each
other. Notice that the pair hopping term in $\mathcal{H}_{\mathbb{I}}$ changes
$N_{\sigma}$ by $\pm 2$ and thus does not mix the even and odd sectors. This
also implies that the strong coupling phase cannot be accessible
perturbatively from the noninteracting limit.

In the weak coupling phase, we find that
$\langle{\bm{S}}_{\mathbb{I}}\cdot{\bm{S}}_{\mathbb{I}}\rangle\sim 3$,
corresponding to $\bar{S}_{\mathbb{I}}\sim 1.3$, which is different from
$S_{\mathbb{I}}=1$ for the ground-state local multiplet $(\mathbb{A},1)$ in
the weak coupling phase found in Fig. 2(c) and Fig. 4(a). The deviation is
attributed to the fact that the total number of electrons,
$N_{\mathbb{I}}=\sum_{i\in\mathbb{I}}\sum_{\alpha=\sigma,\pi}n_{i\alpha}$, is
larger than six. Indeed, we find here that $\langle
N_{\mathbb{I}}\rangle=3(\langle n_{i\sigma}\rangle+\langle
n_{i\pi}\rangle)\sim 7.4$ for $U=0.5t$.

On the other hand, in the strong coupling phase, we find that
$\langle{\bm{S}}_{\mathbb{I}}\cdot{\bm{S}}_{\mathbb{I}}\rangle\sim 2$,
corresponding to $\bar{S}_{\mathbb{I}}\sim 1$, which is again different from
$\bar{S}_{\mathbb{I}}=0$ for the ground-state local multiplet $(\mathbb{A},0)$
in the strong coupling limit found in Fig. 2(d) and Fig. 4(a). One should
notice here that the total number $\langle N_{\mathbb{I}}\rangle$ of electrons
is now as large as 6. However, the expectation value of the local spin squared
$\langle{\bm{S}}_{i}\cdot{\bm{S}}_{i}\rangle$ increases with $U$ and becomes
as larger as 1.5 (i.e., $\bar{S}_{i}\sim 0.82$) at $U=4t$, suggesting the
tendency to form the local site spin 1 with further increasing $U$, as
expected in the strong coupling limit shown in Fig. 4(b). Therefore, these
results reveal that the dominant local multiplet structure found here in the
ground state of $\mathcal{H}_{\rm AM}$ is different from the ground-state
local multiplet $(\mathbb{A},0)$ shown in Fig. 2(d) and Fig. 4(a) due to the
effect of the surrounding $\pi$ electrons in the conduction band, which
deserves more analysis.

Figure 8 shows the projected weight $P(\mathbb{C},S_{\mathbb{I}})$ of the
ground state of the Anderson model $\mathcal{H}_{\rm AM}$ in the strong
coupling phase onto the local multiplet states $(\mathbb{C},S_{\mathbb{I}})$
with $\mathbb{C}=\mathbb{A},\mathbb{E}$ and $S_{\mathbb{I}}=0,1,2$ of
$\mathcal{H}_{\mathbb{I}}$ in the strong coupling limit [see Fig. 2(d)]. Note
that the projected weight for the local multiplet state $(\mathbb{A},3)$ is
not shown in Fig. 8 (and also in Fig. 11) because we find that it is quite
small (less than $10^{-4}$). The explicit form of the local multiplet states
and the definition of the projected weight $P(\mathbb{C},S_{\mathbb{I}})$ are
given in Eqs. (36)–(42) and Eq. (49), respectively, in Appendix A. It is
clearly found in Fig. 8 that the projected weight for $(\mathbb{E},1)$ is
dominant over that for $(\mathbb{A},0)$ at large $L$, suggesting that the
coupling to the surrounding $\pi$ electrons in the conduction band through the
hybridization function $\tilde{\Gamma}_{e}(\omega)$ enhances the local
multiplet component with $(\mathbb{E},1)$.

Figure 8:  $L$ dependence of the projected weight
$P(\mathbb{C},S_{\mathbb{I}})$ of the ground state of the Anderson model
$\mathcal{H}_{\rm AM}$ onto the local multiplet states
$(\mathbb{C},S_{\mathbb{I}})$ of $\mathcal{H}_{\mathbb{I}}$ in the strong
coupling limit. Here, the ground state of the Anderson model $\mathcal{H}_{\rm
AM}$ is obtained by mapping it to the Q1D model $\mathcal{H}_{\rm Q1D}$ with
$L$ at $U=3t$ in the strong coupling phase. The other parameters used are
$J=0.1U$, $t=6t_{\sigma}=12t_{\pi}$, and $\mu=2t_{\pi}$, corresponding to the
undoped case.

To examine whether or not the local spin $\bar{S}_{\mathbb{I}}$ around the
vacancy found in Fig. 7 is screened by the surrounding $\pi$ electrons, we now
calculate the expectation value of $S^{z}_{\mathbb{I}}$ for the ground state
when an external magnetic field $h_{L}$ is applied on $\mathbb{I}$ in the
Anderson model $\mathcal{H}_{\rm AM}$:

$\mathcal{H}_{\rm AM}\to\mathcal{H}_{\rm AM}-h_{L}S^{z}_{\mathbb{I}}.$ (28)

Here, the existence or absence of the local free magnetic moment in the
thermodynamic limit can be analyzed by scaling the magnetic field $h_{L}=t/L$
with increasing $L$ when the ground state is calculated by mapping
$\mathcal{H}_{\rm AM}$ into the Q1D model $\mathcal{H}_{\rm Q1D}$ with $L$.
Figure 9 shows the $L$ dependence of $\langle S_{\mathbb{I}}^{z}\rangle$. In
the weak coupling phase ($U=0.5t$), we find that $\langle
S_{\mathbb{I}}^{z}\rangle$ approaches to zero with decreasing $h_{L}$ (i.e.,
increasing the system size $L$), indicating that the local spin
$\bar{S}_{\mathbb{I}}\sim 1.3$ around the vacancy is completely screened by
the surrounding $\pi$ electrons. In contrast, in the strong coupling phase
($U=3t$), Fig. 9 clearly shows that $\langle S_{\mathbb{I}}^{z}\rangle$ in the
limit of $L\to\infty$ approaches to a finite value, as large as $0.55$,
indicating that the local spin $\bar{S}_{\mathbb{I}}\sim 1$ around the vacancy
is partially screened by the surrounding $\pi$ electrons.

Figure 9:  Finite size scaling of $\langle S_{\mathbb{I}}^{z}\rangle$ by
applying an external magnetic field $h_{L}=t/L$ for the weak coupling phase at
$U=0.5t$ and the strong coupling phase at $U=3t$. The other parameters are the
same as in Fig. 8.

###  IV.2 Doped case

Let us now examine the ground state properties for the doped case. For this
purpose, we set that the number $N_{e}$ of electrons is $N_{e}=3(L+1)$ in the
Q1D model $\mathcal{H}_{\rm Q1D}$ with $L$. This corresponds to a case of hole
doping. The results of the local quantities are summarized in Fig. 10. As in
the undoped case, we find that there are two distinguishable phases, the weak
coupling phase and the strong coupling phase, separated at $U\approx 1.8t$.

Figure 10:  Same as Fig. 7 except for $L=39$ and $N_{e}=3(L+1)$, corresponding
to the doped case.

Figures 11 and 12 show the $L$ dependence of the projected weight
$P(\mathbb{C},S_{\mathbb{I}})$ and $\langle S_{\mathbb{I}}^{z}\rangle$ induced
by the external magnetic field $h_{L}=t/L$, respectively, calculated at $U=3t$
in the strong coupling phase. These results correspond to Figs. 8 and 9 for
the undoped case. Comparing these results, we find that the overall features
in the strong coupling phase are qualitatively the same as those for the
undoped case. Namely, the emergence of the local free magnetic moment is
robust under the carrier doping and the local spin $\bar{S}_{\mathbb{I}}\sim
1$ around the vacancy is only partially screened by the surrounding $\pi$
electrons. Indeed, the difference between the undoped and doped cases is
rather qualitative. For example, the unscreened spin $\langle
S_{\mathbb{I}}^{z}\rangle$ by the surrounding $\pi$ electrons for the doped
case is as large as $\langle S_{\mathbb{I}}^{z}\rangle\sim 0.4$ (see Fig. 12),
which is slightly smaller than that for the undoped case shown in Fig. 9.
Interestingly, this tendency has also been observed in the experiments Nair
_et al._ (2013).

Figure 11:  Same as Fig. 8 except for $N_{e}=3(L+1)$, corresponding to the
doped case.

On the other hand, the noticable difference is found in the weak coupling
phase. As shown in Fig. 10, $\langle{n}_{i\sigma}\rangle\sim 2/3$ and
$\langle{n}_{i\pi}\rangle\sim 1$ for $i\in\mathbb{I}$ in the weak coupling
phase, and thus the total number $N_{\mathbb{I}}$ of electrons in $\mathbb{I}$
is $\langle N_{\mathbb{I}}\rangle\sim 5$, which is less than that for the
undoped case (see Fig. 7) and is also less than 6 assumed in Sec. II.2. This
is probably the reason why the expectation value of total spin squared
${\bm{S}}_{\mathbb{I}}\cdot{\bm{S}}_{\mathbb{I}}$ is smaller than $2$, but
$\langle{\bm{S}}_{\mathbb{I}}\cdot{\bm{S}}_{\mathbb{I}}\rangle\sim 1.2$,
corresponding to $\bar{S}_{\mathbb{I}}\sim 0.7$, in the weak coupling phase.
However, as in the undoped case, the local magnetic moment in the weak
coupling phase vanishes due to the coupling with the surrounding $\pi$
electron system, as shown in Fig. 12. Therefore, we can conclude that the weak
coupling phase is always nonmagnetic for the undoped and doped cases.

Figure 12:  Same as Fig. 9 except for $U=t$ in the weak coupling phase and
$N_{e}=3(L+1)$, corresponding to the doped case.

###  IV.3 Asymptotic behavior of spin correlation function

Let us now examine the spin correlation function
$\langle{\bm{S}}_{\alpha}\cdot\tilde{\bm{S}}_{lc}\rangle$ between the impurity
sites and the conduction sites in the Q1D model $\mathcal{H}_{\rm Q1D}$. Here,
${\bm{S}}_{\alpha}$ ($\alpha=\sigma,\pi$) is the total spin operator of the
$\alpha$ orbitals around the vacancy, i.e.,

${\bm{S}}_{\sigma}=\sum_{i\in\mathbb{I}}{\bm{S}}_{i\sigma}$ (29)

and

${\bm{S}}_{\pi}=\sum_{i\in\mathbb{I}}{\bm{S}}_{i\pi},$ (30)

and
$\tilde{\bm{S}}_{lc}=(\tilde{S}^{x}_{lc},\tilde{S}^{y}_{lc},\tilde{S}^{z}_{lc})$
for $c=a,e_{1},e_{2}$ is the local spin operator of the conduction sites in
the Q1D model $\mathcal{H}_{\rm Q1D}$, i.e.,

$\tilde{S}_{lc}^{\mu}=\frac{1}{2}\sum_{s,s^{\prime}=\uparrow,\downarrow}\tilde{f}_{lcs}^{\dagger}\left[\mbox{\boldmath{$\sigma$}}^{\mu}\right]_{ss^{\prime}}\tilde{f}_{lcs^{\prime}}$
(31)

for $l\geq 2$ [see Fig. 6(c)]. Note that
$\sum_{c=a,e_{1},e_{2}}\tilde{\bm{S}}_{1c}={\bm{S}}_{\pi}$ if we use Eq. (31)
for $l=1$ because the orbitals at $l=1$ correspond to the $\pi$ orbitals in
$\mathbb{I}$ around the vacancy, as shown in Fig. 6. Since the results for the
$c=e_{1}$ and $e_{2}$ modes are identical, here we show only one of these
results and refer to them simply as the results for the $c=e$ mode. We also
note that the spin correlation function
$\langle{\bm{S}}_{\alpha}\cdot\tilde{\bm{S}}_{lc}\rangle$ is plotted as a
function of the distance $d$ between the impurity site and the $l$th
conduction site given by $l-1$ in the Q1D model.

Figure 13 shows the spin correlation function
$\langle{\bm{S}}_{\sigma}\cdot\tilde{\bm{S}}_{lc}\rangle$ calculated for the
ground states of the doped and undoped cases in the strong coupling phase. We
can find in Fig. 13 that
$\langle{\bm{S}}_{\sigma}\cdot\tilde{\bm{S}}_{lc}\rangle$ decays as or even
slower than $d^{-1}$ for both $c=a$ and $e$ modes. Note also that such
behavior does not depend qualitatively on the carrier doping, and thus the
spin structure of the $\sigma$ orbitals is robust against the carrier doping.
However, the overall values of the spin correlation function for the $e$ mode
is approximately $10$ times larger than those for the $a$ mode, suggesting
that the local spin ${\bm{S}}_{\sigma}$ of the $\sigma$ orbitals around the
vacancy is mostly coupled to the surrounding $\pi$ orbitals with the $c=e$
mode.

Figure 13:  Log-log plots of the absolute values of the spin correlation
function $\langle{\bm{S}}_{\sigma}\cdot\tilde{\bm{S}}_{lc}\rangle$ with the
$c=a$ and $e$ modes for the doped case (left panels) and the undoped case
(right panels) with various values of $L$ in the Q1D model $\mathcal{H}_{\rm
Q1D}$. The other parameters used are $U=3t$, $J=0.1U$, and
$t=6t_{\sigma}=12t_{\pi}$. Green dashed line (blue dotted line) represents a
function proportional to $(l-1)^{-1}$ [$(l-1)^{-2}$], for comparison.

As an effective model of the Anderson model $\mathcal{H}_{\rm AM}$, we
consider in Appendix B a two-orbital Anderson model
$\mathcal{H}^{({\text{II}})}$ described in Eq. (51), referred to as “model II”
in Appendix B. As schematically shown in Fig. 15(b), the impurity sites in
this model are composed of one of the lattice sites in the honeycomb lattice
that is replaced with a Hubbard site (denoted as impurity site B) and an
additional Hubbard site (denoted as impurity site A) attached on top of the
impurity site B. These impurity sites contains the inter-orbital Coulomb
interactions with no hopping and are coupled to the conduction sites through
the impurity site B with a diverging hybridization function. Therefore,
corresponding the impurity sites A and B in $\mathcal{H}^{({\text{II}})}$ to
the $\sigma$ and $\pi$ orbitals in $\mathbb{I}$ for $\mathcal{H}_{\rm AM}$,
respectively, we can consider this model $\mathcal{H}^{({\text{II}})}$ as a
simple version of the Anderson model $\mathcal{H}_{\rm AM}$.

As shown in Figs. 17(a) and 17(b) for the doped and undoped cases,
respectively, we find that the spin correlation function between the impurity
site A and the conduction sites in the Q1D representation of
$\mathcal{H}^{({\text{II}})}$ [see Fig. 15(d)] decays approximately as
$d^{-1}$, which resemble the results shown in Figs. 13(a) and 13(b). In the
strong coupling phase of this two-orbital Anderson model
$\mathcal{H}^{({\text{II}})}$, the Hund’s coupling favors the formation of the
local spin-1 state at the impurity sites, while the effect of the surrounding
electrons in the conduction sites partially screen this spin-1 state, thus
leading to the emergence of the spin-1/2 free magnetic moment. The partial
screening in this model is robust even for the doped case because of the large
hybridization. These features are indeed similar to the underlying picture
obtained for the Anderson model $\mathcal{H}_{\rm AM}$ in the previous
sections. Namely, the local spin $\bar{S}_{\mathbb{I}}\sim 1$ around the
vacancy is partially screened by the surrounding $\pi$ electrons in the
conduction sites and the residual free magnetic moment is robust against the
doping.

Figure 14 shows the spin correlation function
$\langle{\bm{S}}_{\pi}\cdot\tilde{\bm{S}}_{lc}\rangle$ for the ground states
of the doped and undoped cases of the Q1D model $\mathcal{H}_{\rm Q1D}$ in the
strong coupling phase. We find in Figs. 14(c) and 14(d) that
$\langle{\bm{S}}_{\pi}\cdot\tilde{\bm{S}}_{le}\rangle$ decays approximately as
$d^{-2}$ for both undoped and doped cases. The spin correlation function
$\langle{\bm{S}}_{\pi}\cdot\tilde{\bm{S}}_{la}\rangle$ for the doped case in
Fig. 14(a) also exhibits the $d^{-2}$ decay. However,
$\langle{\bm{S}}_{\pi}\cdot\tilde{\bm{S}}_{la}\rangle$ for the undoped case in
Fig. 14(b) is different and is instead somewhat enhanced as compared with the
other cases: it deviates from the $d^{-2}$ decay but rather approaches to the
$d^{-1}$ decay.

Figure 14:  Same as Fig. 13 except for the spin correlation function
$\langle{\bm{S}}_{\pi}\cdot\tilde{\bm{S}}_{lc}\rangle$ with the $c=a$ and $e$
modes.

As described in Appendix B, the $d^{-2}$ decay of the spin correlation
function implies the absence of the free magnetic moment. We find in Figs.
17(c) and 17(d) that the spin correlation function between the impurity site B
and the conduction sites decays as $d^{-2}$ for both doped and undoped cases
of the two-orbital Anderson model $\mathcal{H}^{({\text{II}})}$. In this
model, the screening mainly occurs at the impurity site B among the two
impurity sites via the diverging hybridization function. This is indeed
similar to the features found here for the Anderson model $\mathcal{H}_{\rm
AM}$ because the local spin $\bar{S}_{\mathbb{I}}\sim 1$ around the vacancy is
partially screened by the surrounding $\pi$ electrons in the conduction band
via the diverging hybridization function $\tilde{\Gamma}_{e}(\omega)$, leading
to the $d^{-2}$ decay of the spin correlation function
$\langle{\bm{S}}_{\pi}\cdot\tilde{\bm{S}}_{le}\rangle$.

In contrast, the spin correlation function
$\langle{\bm{S}}_{\pi}\cdot\tilde{\bm{S}}_{la}\rangle$ in the largest system
for the undoped case shown in Fig. 14(b) decays as $d^{-1.43}$, certainly
slower than $d^{-2}$, which is a hallmark of the presence of the free magnetic
moment (see the discussion in Appendix B). This can be attributed to the fact
that the hybridization function $\tilde{\Gamma}_{a}(\omega)$ exhibits the
pseudogap structure at the Fermi level for the undoped case, as shown in Fig.
5. It is known that the Kondo screening does not occur in the undoped
pseudogap Anderson model because the density of states vanishes at the Fermi
energy Bulla  _et al._ (1997); Fritz and Vojta (2013); Shirakawa and Yunoki
(2016); Gonzalez-Buxton and Ingersent (1998); Vojta  _et al._ (2010). In this
case, there are two possible phases that arises when the carriers are
introduced by varying the chemical potential away from the pseudogap position.
One is the Kondo screening phase simply because the density of state becomes
finite when the chemical potential is shifted from the pseudogap position
Vojta  _et al._ (2010). The other is the asymmetric strong coupling phase, in
which the electronic state at the impurity site becomes either empty or doubly
occupied Vojta  _et al._ (2010). Both in the Kondo screening and asymmetric
strong coupling phases, the magnetic moment vanishes. In either case,
therefore, the free magnetic moment is sensitive against the carrier doping
(see Appendix B for the numerical study of a simple model for the pseudogap
Anderson system). This might explain why the residual free magnetic moment
around the vacancy is slightly reduced when the carriers are introduced in the
strong coupling phase of the Anderson model $\mathcal{H}_{\rm AM}$ (see Figs.
9 and 12).

##  V Summary and Discussion

We have studied the ground state properties of the effective Anderson model
for a single vacancy in graphene. The Anderson model considered here is
composed of the three dangling $\sigma$ orbitals as well as the three $\pi$
orbitals of carbon atoms around the vacancy, treated as the impurity sites,
which are coupled to the noninteracting conduction electrons in the
surrounding $\pi$ orbitals of carbon atoms forming the honeycomb lattice with
the single vacancy. Employing the block Lanczos DMRG method, we have found
that the ground state in the reasonable parameter region for graphene,
belonging to the strong coupling phase, shows the emergence of the free
magnetic moment with its spin as large as $\sim 1/2$ around the vacancy. Our
systematic analysis of the projected weight onto the local multiplet states
has shown that the local spin ${S}_{\mathbb{I}}=1$ multiplet with the doubly
degenerate $\mathcal{E}$ representation of the $\mathcal{C}_{3}$ point group
is dominant in the ground state, and thus the residual free magnetic moment is
attributed to the partial screening of this ${S}_{\mathbb{I}}=1$ multiplet by
the sorrounding $\pi$ electrons in the conduction band. We have also found
that the emergent local free magnetic moment is robust against the hole
doping. The calculations of the spin correlation function have shown that the
screening involves mostly the local $\pi$ orbitals around the vacancy and thus
the residual free magnetic moment is composed mainly of the local $\sigma$
orbitals. Because these $\sigma$ orbitals are decoupled with no happing terms
to the surrounding $\pi$ orbital system, the emergent local free magnetic
moment is not sensitive to any change of the surrounding $\pi$ orbital system
such as the carrier doping. However, our calculations of the spin correlation
function have also shown the precursor indicating another contribution to the
emergent local free magnetic moment in the undoped case. This contribution
arises due to the pseudogap structure of the hybridization function
$\tilde{\Gamma}_{a}(\omega)$ in the $\pi$ orbital system and thus sensitive to
the carrier doping. The presence of these two contributions to the emergent
local free magnetic moment fits nicely the experimental observations Nair  _et
al._ (2013); Zhang  _et al._ (2016) and their scenario proposed in Ref. Nair
_et al._ (2013).

It should be noted that the dominant contribution to the local multiplet
structure around the vacancy, having the local spin 1 and the $\mathcal{E}$
irreducible representation, is compatible with the occurrence of the in-plane
Jahn-Teller distortion found in the previous ab-initio DFT calculations Paz
_et al._ (2013); Padmanabhan and Nanda (2016), in a sense that the Jahn-Teller
distortion occurs through the coupling with the in-plane $e$ vibration mode
Casartelli  _et al._ (2013), similar to the so-called $\mathbb{E}\otimes e$
problem Bersuker (2001). However, our results imply that the Jahn-Teller
distortion is not necessarily required to form the local spin triplets and the
residual free magnetic moment as large as spin $S=1/2$.

In contrast, the out-of-plane distortion generates a finite hybridization
between the neighboring $\sigma$ and $\pi$ orbitals. As a consequence, a
nonmagnetic ground state can be realized because the residual free magnetic
moment composed mainly of the $\sigma$ orbitals around the vacancy would be
screened by the surrounding $\pi$ electrons via this finite hybridization
Kanao  _et al._ (2012); Shirakawa and Yunoki (2014). This conjecture is also
in accordance with the previous ab-initio DFT results Paz  _et al._ (2013);
Padmanabhan and Nanda (2016), showing that the out-of-plane distortion gives
rise to a nonmagnetic ground state. This may bring an interesting
functionality by controlling the out-of-plane strain field suggested by the
ab-initio DFT calculation Dharma-wardana and Zgierski (2008).

The scanning tunneling spectroscopy measurements of a single vacancy in
graphene have revealed two peaks in the tunneling spectrum as a function of
bias voltage that are originated from the spin-polarized states around the
vacancy Zhang  _et al._ (2016). They also found that these two peaks are quite
robust against the carrier doping into graphene Nair  _et al._ (2013); Zhang
_et al._ (2016). This is in sharp contrast to the magnetism in graphene
induced by hydrogen absorption González-Herrero  _et al._ (2016), where the
magnetic moment is sensitive to the carrier doping. This is also supported by
the theoretical analysis on the carrier doping in the pseudogap Kondo model
Vojta  _et al._ (2010). Our results in this paper as well as the previous
calculations of effective models for a single adatom in graphene Shirakawa and
Yunoki (2014, 2016) suggest that the magnetism observed in these two systems
is caused by different mechanisms. Here, we have focused on static physical
quantities because the calculation of dynamical quantities is rather
demanding. However, it is highly valuable to examine the dynamical quantities
such as the local density of states and the response to an external field
including a gauge field, which can shed light on the relation between the
local multiplet structures and their experimental observation. These studies
are left for the future.

## Acknowledgements

The authors are grateful to Beom Hyn Kim for fruitful discussions. A part of
the numerical simulations has been performed using the HOKUSAI supercomputer
at RIKEN (Project ID: Q21532) and also supercomputer Fugaku installed in RIKEN
R-CCS. This work is supported by Grant-in-Aid for Scientific Research (B) (No.
JP18H01183) and Grant-in-Aid for Scientific Research (A) (No. JP21H04446) from
MEXT, Japan. This work is also supported in part by the COE research grant in
computational science from Hyogo Prefecture and Kobe City through Foundation
for Computational Science.

##  Appendix A Local multiplet states in the strong coupling limit and
projected weight

In this appendix, we shall construct explicitly the eigenstates of the local
part $\mathcal{H}_{\mathbb{I}}$ of the Anderson model $\mathcal{H}_{\rm AM}$
in the strong coupling limit, assuming that the number of electrons is six. We
also provide the definition of the projected weight
$P(\mathbb{C},S_{\mathbb{I}})$ with $\mathbb{C}=\mathbb{A},\mathbb{E}$.

###  A.1 Local multiplet states

As described in Sec. II.2, the low-lying states in the strong coupling limit
form the local spin one at each carbon site that is occupied by two electrons.
Therefore, the local Hilbert space at site $i$ ($i=1,2,3$) $\in\mathbb{I}$
[see Fig. 1(c)] can be spanned by the following bases:

$\displaystyle|1\rangle_{i}=c_{i\sigma\uparrow}^{\dagger}c_{i\pi\uparrow}^{\dagger}|0\rangle,$
$\displaystyle|0\rangle_{i}=\frac{1}{\sqrt{2}}(c_{i\sigma\uparrow}^{\dagger}c_{i\pi\downarrow}^{\dagger}+c_{i\sigma\downarrow}^{\dagger}c_{i\pi\uparrow}^{\dagger})|0\rangle,$
(32)
$\displaystyle|-1\rangle_{i}=c_{i\sigma\downarrow}^{\dagger}c_{i\pi\downarrow}^{\dagger}|0\rangle,$

where $|S_{i}^{z}\rangle_{i}$ with $S_{i}^{z}=-1,0,1$ is the simultaneous
eigenstate of the total spin operator ${\bm{S}}_{i}$ defined in Eq. (7) and
the $z$-component of ${\bm{S}}_{i}$ at site $i$ with its eigenvalues 1 and
$S_{i}^{z}$, respectively, and $c_{i\alpha s}^{\dagger}$ is the creation
operator of an electron at site $i$ with spin $s$ ($=\uparrow,\downarrow$) and
orbital $\alpha$ (=$\sigma$, $\pi$).

The effective low-energy Hamiltonian is then given by the three-site spin-one
antiferromagnetic Heisenberg model,

$\mathcal{H}_{\rm S}=J_{\rm
ex}\sum_{\scriptsize\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}i,j\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\mathbb{I}}}{\bm{S}}_{i}\cdot{\bm{S}}_{j},$
(33)

where $J_{\rm ex}$ is an effective exchange interaction,
$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}i,j\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}_{\mathbb{I}}$
denotes a pair of neighboring sites $i$ and $j$ around the vacancy, and
${\bm{S}}_{i}$ is the spin-1 operator at site $i\in\mathbb{I}$ with its
eigenstates being given by Eq. (32). For the three-site case, this Hamiltonian
is particularly simplified by introducing the total spin operator
${\bm{S}}_{\mathbb{I}}=\sum_{i\in\mathbb{I}}{\bm{S}}_{i}=(S_{\mathbb{I}}^{x},S_{\mathbb{I}}^{y},S_{\mathbb{I}}^{z})$
as

$\mathcal{H}_{\rm S}=\frac{J_{\rm
ex}}{2}\left({\bm{S}}_{\mathbb{I}}\cdot{\bm{S}}_{\mathbb{I}}-6\right).$ (34)

Therefore, we can readily find that the eigenvalues of $\mathcal{H}_{\rm S}$
are determined solely by the quantum number of ${\bm{S}}_{\mathbb{I}}$,
denoted as $S_{\mathbb{I}}$. The eigenvalues $E(S_{\mathbb{I}})$ of
$\mathcal{H}_{\rm S}$ are thus given as

$\begin{split}&E(S_{\mathbb{I}}=0)=-3J_{\rm ex},\\\
&E(S_{\mathbb{I}}=1)=-2J_{\rm ex},\\\ &E(S_{\mathbb{I}}=2)=0,\\\
&E(S_{\mathbb{I}}=3)=3J_{\rm ex}.\\\ \end{split}$ (35)

Because of the $120^{\circ}$ rotational symmetry around the center of the
cluster, the eigenstates of $\mathcal{H}_{\rm S}$ are also characterized by
the irreducible representations of the $\mathcal{C}_{3}$ point group.
Therefore, we can denote the energy eigenstate as
$|C,S_{\mathbb{I}},M_{\mathbb{I}}\rangle$, where $C=A$ ($C=E_{k}$ with
$k=1,2$) indicates the symmetric (doubly degenerate) irreducible
representation $\mathcal{A}$ ($\mathcal{E}$) of the $\mathcal{C}_{3}$ point
group, and $M_{\mathbb{I}}$ denotes the quantum number of the $z$-component
$S^{z}_{\mathbb{I}}$ of the total spin operator ${\bm{S}}_{i}$. The lowest
eigenstate is $(\mathbb{A},0)$ [for notation, see Sec. II.2 and Fig. 2(d)] and
is given by

$\displaystyle|A,0,0\rangle$
$\displaystyle=\frac{1}{\sqrt{6}}\left[|0\rangle_{1}|1\rangle_{2}|-1\rangle_{3}-|0\rangle_{1}|-1\rangle_{2}|1\rangle_{3}\right]$
$\displaystyle+\frac{1}{\sqrt{6}}\left[|-1\rangle_{1}|0\rangle_{2}|1\rangle_{3}-|1\rangle_{1}|0\rangle_{2}|-1\rangle_{3}\right]$
$\displaystyle+\frac{1}{\sqrt{6}}\left[|1\rangle_{1}|-1\rangle_{2}|0\rangle_{3}-|-1\rangle_{1}|1\rangle_{2}|0\rangle_{3}\right],$
(36)

where $|S_{i}^{z}\rangle_{i}$ with $S_{i}^{z}=-1,0,1$ and $i=1,2,3$ is given
in Eq. (32).

The second lowest eigenstates are those with $S_{\mathbb{I}}=1$, which are
futher classified by the irreducible representation $C$, i.e.,
$(\mathbb{A},1)$ and $(\mathbb{E},1)$. The eigenstate for $(\mathbb{A},1)$
with $M_{\mathbb{I}}^{z}=0$ is given by

$\displaystyle|A,1,0\rangle$
$\displaystyle=\frac{1}{\sqrt{15}}\left[|0\rangle_{1}|1\rangle_{2}|-1\rangle_{3}+|0\rangle_{1}|-1\rangle_{2}|1\rangle_{3}\right]$
$\displaystyle+\frac{1}{\sqrt{15}}\left[|-1\rangle_{1}|0\rangle_{2}|1\rangle_{3}+|1\rangle_{1}|0\rangle_{2}|-1\rangle_{3}\right]$
$\displaystyle+\frac{1}{\sqrt{15}}\left[|1\rangle_{1}|-1\rangle_{2}|0\rangle_{3}+|-1\rangle_{1}|1\rangle_{2}|0\rangle_{3}\right]$
$\displaystyle-\frac{\sqrt{3}}{\sqrt{5}}|0\rangle_{1}|0\rangle_{2}|0\rangle_{3}.$
(37)

The eigenstates with different $M_{\mathbb{I}}$ ($=\pm 1$) are obtained by
multiplying $S_{\mathbb{I}}^{\pm}=S_{\mathbb{I}}^{x}\pm{\rm
i}S_{\mathbb{I}}^{y}$ operators. Similarly, the other two eigenstates for
$(\mathbb{E},1)$ with $M_{\mathbb{I}}=0$ are given by

$\displaystyle|E_{1},1,0\rangle$
$\displaystyle=\frac{2}{2\sqrt{3}}\left[|0\rangle_{1}|1\rangle_{2}|-1\rangle_{3}+|0\rangle_{1}|-1\rangle_{2}|1\rangle_{3}\right]$
$\displaystyle-\frac{1}{2\sqrt{3}}\left[|-1\rangle_{1}|0\rangle_{2}|1\rangle_{3}+|1\rangle_{1}|0\rangle_{2}|-1\rangle_{3}\right]$
$\displaystyle-\frac{1}{2\sqrt{3}}\left[|1\rangle_{1}|-1\rangle_{2}|0\rangle_{3}+|-1\rangle_{1}|1\rangle_{2}|0\rangle_{3}\right]$
(38)

and

$\displaystyle|E_{2},1,0\rangle$
$\displaystyle=\frac{1}{2}\left[|-1\rangle_{1}|0\rangle_{2}|1\rangle_{3}+|1\rangle_{1}|0\rangle_{2}|-1\rangle_{3}\right]$
$\displaystyle-\frac{1}{2}\left[|1\rangle_{1}|-1\rangle_{2}|0\rangle_{3}+|-1\rangle_{1}|1\rangle_{2}|0\rangle_{3}\right].$
(39)

The third eigenstates are those for $(\mathbb{E},2)$. The eigenstates with
$M_{\mathbb{I}}=0$ are given by

$\displaystyle|E_{1},2,0\rangle$
$\displaystyle=\frac{2}{2\sqrt{3}}\left[|0\rangle_{1}|1\rangle_{2}|-1\rangle_{3}-|0\rangle_{1}|-1\rangle_{2}|1\rangle_{3}\right]$
$\displaystyle-\frac{1}{2\sqrt{3}}\left[|-1\rangle_{1}|0\rangle_{2}|1\rangle_{3}-|1\rangle_{1}|0\rangle_{2}|-1\rangle_{3}\right]$
$\displaystyle-\frac{1}{2\sqrt{3}}\left[|1\rangle_{1}|-1\rangle_{2}|0\rangle_{3}-|-1\rangle_{1}|1\rangle_{2}|0\rangle_{3}\right]$
(40)

and

$\displaystyle|E_{2},2,0\rangle$
$\displaystyle=\frac{1}{2}\left[|-1\rangle_{1}|0\rangle_{2}|1\rangle_{3}-|1\rangle_{1}|0\rangle_{2}|-1\rangle_{3}\right]$
$\displaystyle-\frac{1}{2}\left[|1\rangle_{1}|-1\rangle_{2}|0\rangle_{3}-|-1\rangle_{1}|1\rangle_{2}|0\rangle_{3}\right].$
(41)

Finally, the fourth eigenstates are those for $(\mathbb{A},3)$. The
eigenstates with $M_{\mathbb{I}}=0$ is given by

$\displaystyle|A,3,0\rangle$
$\displaystyle=\frac{1}{\sqrt{10}}\left[|0\rangle_{1}|1\rangle_{2}|-1\rangle_{3}+|0\rangle_{1}|-1\rangle_{2}|1\rangle_{3}\right]$
$\displaystyle+\frac{1}{\sqrt{10}}\left[|-1\rangle_{1}|0\rangle_{2}|1\rangle_{3}+|1\rangle_{1}|0\rangle_{2}|-1\rangle_{3}\right]$
$\displaystyle+\frac{1}{\sqrt{10}}\left[|1\rangle_{1}|-1\rangle_{2}|0\rangle_{3}+|-1\rangle_{1}|1\rangle_{2}|0\rangle_{3}\right]$
$\displaystyle+\frac{\sqrt{2}}{\sqrt{10}}|0\rangle_{1}|0\rangle_{2}|0\rangle_{3}.$
(42)

###  A.2 Projected weight

Let $\mathcal{Q}^{\dagger}_{CS_{\mathbb{I}}M_{\mathbb{I}}}$ be an operator
that generates the eigenstate
$|C,S_{\mathbb{I}},M_{\mathbb{I}}\rangle_{\mathbb{I}}$ from the vacuum state
$|0\rangle_{\mathbb{I}}$ in the Hilbert space of the region $\mathbb{I}$,
i.e.,

$\displaystyle|C,S_{\mathbb{I}},M_{\mathbb{I}}\rangle_{\mathbb{I}}=\mathcal{Q}^{\dagger}_{CS_{\mathbb{I}}M_{\mathbb{I}}}|0\rangle_{\mathbb{I}},$
(43)

where $|C,S_{\mathbb{I}},M_{\mathbb{I}}\rangle_{\mathbb{I}}$ ($C=A$, $E_{1}$,
$E_{2}$) is given in Eqs. (36)-(42). Replacing all
$|S_{1}^{z}\rangle_{1}|S_{2}^{z}\rangle_{2}|S_{3}^{z}\rangle_{3}$ in
$|C,S_{\mathbb{I}},M_{\mathbb{I}}\rangle_{\mathbb{I}}$ with the electron
creation operators as in Eqs. (32),
$\mathcal{Q}^{\dagger}_{CS_{\mathbb{I}}M_{\mathbb{I}}}$ are now represented as
the linear combination of products of six electron creation operators. We then
introduce operators $\mathcal{N}_{0}$ and
$\mathcal{P}_{CS_{\mathbb{I}}M_{\mathbb{I}}}^{\dagger}$ defined respectively
by

$\mathcal{N}_{0}=\prod_{i\in\mathbb{I}}\prod_{\alpha=\sigma,\pi}\prod_{s=\uparrow,\downarrow}c_{i\alpha
s}c_{i\alpha s}^{\dagger}$ (44)

and

$\mathcal{P}^{\dagger}_{CS_{\mathbb{I}}M_{\mathbb{I}}}=\mathcal{Q}^{\dagger}_{CS_{\mathbb{I}}M_{\mathbb{I}}}\mathcal{N}_{0}.$
(45)

We can readily find that $\mathcal{N}_{0}$ is a projection operator onto the
vacuum state in the region $\mathbb{I}$. Note also that
$\mathcal{P}^{\dagger}_{CS_{\mathbb{I}}M_{\mathbb{I}}}$ can be represented as
the linear combination of products of eighteen electron creation and
annihilation operators.

Let $|\Phi_{0}\rangle$ be the ground state of the whole system of the Anderson
model $\mathcal{H}_{\rm AM}$. In general, the ground state $|\Phi_{0}\rangle$
can be written as

$|\Phi_{0}\rangle=\sum_{n,l}\phi_{nl}|\psi_{n}\rangle_{\mathbb{I}}|\tilde{\psi}_{l}\rangle_{\mathbb{L}},$
(46)

where $|\psi_{n}\rangle_{\mathbb{I}}$ and
$|\tilde{\psi}_{l}\rangle_{\mathbb{L}}$ are orthonormalized basis states in
the regions $\mathbb{I}$ and $\mathbb{L}$, respectively, and we assume that
the ground state $|\Phi_{0}\rangle$ is normalized. The reduced density matrix
operator $\rho_{\mathbb{I}}$ in the region $\mathbb{I}$ is then given by

$\rho_{\mathbb{I}}={\rm Tr}_{\mathbb{L}}|\Phi_{0}\rangle\langle\Phi_{0}|,$
(47)

where ${\rm Tr}_{\mathbb{L}}\cdots$ denotes the trace over all the basis
states in the region $\mathbb{L}$. It is now easy to show that
${}_{\mathbb{I}}\langle\psi_{n}|\rho_{\mathbb{I}}|\psi_{n}\rangle_{\mathbb{I}}=\sum_{l}|\phi_{nl}|^{2}$.
The projected weight $P(\mathbb{C},S_{\mathbb{I}})$ discussed in Secs. IV.1
and IV.2 are defined by

$\begin{split}P(\mathbb{A},S_{\mathbb{I}})&=\sum_{M_{\mathbb{I}}=-S_{\mathbb{I}}}^{S_{\mathbb{I}}}\langle
A,S_{\mathbb{I}},M_{\mathbb{I}}|\rho_{\mathbb{I}}|A,S_{\mathbb{I}},M_{\mathbb{I}}\rangle,\\\
P(\mathbb{E},S_{\mathbb{I}})&=\sum_{k=1}^{2}\sum_{M_{\mathbb{I}}=-S_{\mathbb{I}}}^{S_{\mathbb{I}}}\langle
E_{k},S_{\mathbb{I}},M_{\mathbb{I}}|\rho_{\mathbb{I}}|E_{k},S_{\mathbb{I}},M_{\mathbb{I}}\rangle.\end{split}$
(48)

To reduce the numerical complexity of the DMRG calculations, we evaluate the
following equivalent quantities

$\begin{split}P(\mathbb{A},S_{\mathbb{I}})&=\sum_{M_{\mathbb{I}}=-S_{\mathbb{I}}}^{S_{\mathbb{I}}}\langle\Phi_{0}|\mathcal{P}_{AS_{\mathbb{I}}M_{\mathbb{I}}}^{\dagger}\mathcal{P}_{AS_{\mathbb{I}}M_{\mathbb{I}}}|\Phi_{0}\rangle,\\\
P(\mathbb{E},S_{\mathbb{I}})&=\sum_{k=1}^{2}\sum_{M_{\mathbb{I}}=-S_{\mathbb{I}}}^{S_{\mathbb{I}}}\langle\Phi_{0}|\mathcal{P}_{E_{k}S_{\mathbb{I}}M_{\mathbb{I}}}^{\dagger}\mathcal{P}_{E_{k}S_{\mathbb{I}}M_{\mathbb{I}}}|\Phi_{0}\rangle\end{split}$
(49)

by using the multitarget technique Hallberg (1995), for which the ground state
$|\Phi_{0}\rangle$ and
$\mathcal{P}_{C{S}_{\mathbb{I}}M_{\mathbb{I}}}|\Phi_{0}\rangle$ are included
as the target states.

##  Appendix B Asymptotic behavior of spin correlation function for reference
models

In order to better understand the results for the Anderson model
$\mathcal{H}_{\rm AM}$ in Sec. IV.3, here in this appendix we consider two
simpler Anderson models (model I and model II) as reference systems and
provide the results of the spin correlation function between the impurity and
conduction sites for these models.

As shown schematically in Fig. 15(a), model I is a single-impurity Anderson
model with a single impurity site attached to one of the conduction sites in
the honeycomb lattice. The Hamiltonian $\mathcal{H}^{(\text{I})}$ of model I
is given by

$\displaystyle\mathcal{H}^{(\text{I})}$ $\displaystyle=Un_{{\rm
H}\uparrow}n_{{\rm H}\downarrow}-\frac{U}{2}(n_{{\rm H}\uparrow}+n_{{\rm
H}\downarrow}),$ $\displaystyle-V\sum_{s=\uparrow,\downarrow}(c_{{\rm
H}s}^{\dagger}c_{i_{0}\pi s}+c_{i_{0}\pi s}^{\dagger}c_{{\rm H}s})$
$\displaystyle-t\sum_{\langle i,j\rangle}\sum_{s=\uparrow,\downarrow}(c_{i\pi
s}^{\dagger}c_{j\pi s}+c_{j\pi s}^{\dagger}c_{i\pi s})$ $\displaystyle-
t_{\pi}\sum_{\mbox{\scriptsize{$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}$}}i,j\mbox{\scriptsize{$\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}$}}}\sum_{s=\uparrow,\downarrow}(c_{i\pi
s}^{\dagger}c_{j\pi s}+c_{j\pi s}^{\dagger}c_{i\pi s}),$ (50)

where $c_{{\rm H}s}$ ($c_{{\rm H}s}^{\dagger}$) is the annihilation (creation)
operator of an electron at the impurity site H with spin
$s\,(=\uparrow,\downarrow)$ and $n_{{\rm H}s}=c_{{\rm H}s}^{\dagger}c_{{\rm
H}s}$. $c_{i\pi s}$ ($c_{i\pi s}^{\dagger}$) is the annihilation (creation)
operator of an electron at conduction site $i$ with spin
$s\,(=\uparrow,\downarrow)$. $i_{0}$ labels the conduction site that is
connected to the impurity site H through the hopping $V$. As in the Anderson
model $\mathcal{H}_{\rm AM}$, the nearest and next nearest neighboring
hoppings $t$ and $t_{\pi}$, respectively, are considered in the last two terms
in Eq. (50).

Figure 15:  Schematic depictions of (a) model I (a single-impurity Anderson
model) and (b) model II (a two-impurity Anderson model). The impurity sites
are indicated with “H” in (a) and “A” and “B” in (b). Schematic depictions of
the Q1D representation of (c) model I and (d) model II. Blue and orange
spheres denote the impurity and conduction sites, respectively. Black and bold
magenta lines indicate lattice bonds connecting two neighboring sites via the
hopping and two-body interactions, respectively. Local density of states at
$l=2$ without the impurity site in (c) shows a pseudogap structure, while
local density of states at $l=3$ without the impurity sites in (d) exhibits a
diverging structure, at Fermi level for the undoped case. The impurity site at
$l=1$ in (c) is exactly the same as the impurity site $H$ in (a), and the
impurity sites at $l=1$ and 2 in (d) are exactly the same as the impurity
sites A and B, respectively, in (b).

As shown schematically in Fig. 15(b), model II is a two-impurity (or two-
orbital) Anderson model, where the impurity sites are composed of one of the
lattice sites in the honeycomb lattice (referred to as the impurity site B)
and an additional impurity site (referred to as the impurity site A) attached
on top of the impurity site B through the inter-orbital interactions without
any hopping. Therefore, only the impurity site B is connected to the
conduction sites. Comparing the Anderson model $\mathcal{H}_{\rm AM}$ in Eq.
(1), the impurity site A (B) mimics the $\sigma$ ($\pi$) orbitals in the
impurity sites (i.e., in the region $\mathbb{I}$) for $\mathcal{H}_{\rm AM}$.
The Hamiltonian $\mathcal{H}^{(\text{II})}$ of model II is given by

$\displaystyle\mathcal{H}^{(\text{II})}$ $\displaystyle=U\sum_{\alpha={\rm
A},{\rm
B}}n_{i_{0}\alpha\uparrow}n_{i_{0}\alpha\downarrow}+U^{\prime}n_{i_{0}{\rm
A}}n_{i_{0}{\rm B}}$ $\displaystyle-2J\left({\bm{S}}_{i_{0}{\rm
A}}\cdot{\bm{S}}_{i_{0}{\rm B}}+\frac{1}{4}n_{i_{0}{\rm A}}n_{i_{0}{\rm
B}}\right)$ $\displaystyle+J^{\prime}\left(c_{i_{0}{\rm
A}\uparrow}^{\dagger}c_{i_{0}{\rm A}\downarrow}^{\dagger}c_{i_{0}{\rm
B}\downarrow}c_{i_{0}{\rm B}\uparrow}+c_{i_{0}{\rm
B}\uparrow}^{\dagger}c_{i_{0}{\rm B}\downarrow}^{\dagger}c_{i_{0}{\rm
A}\downarrow}c_{i_{0}{\rm A}\uparrow}\right)$
$\displaystyle-\frac{1}{2}(U+2U^{\prime}-J)\sum_{\alpha={\rm A},{\rm
B}}n_{i_{0}\alpha}$ $\displaystyle-t\sum_{\langle
i,j\rangle}\sum_{s=\uparrow,\downarrow}(c_{i\pi s}^{\dagger}c_{j\pi s}+c_{j\pi
s}^{\dagger}c_{i\pi s})$ $\displaystyle-
t_{\pi}\sum_{\mbox{\scriptsize{$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}$}}i,j\mbox{\scriptsize{$\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}$}}}\sum_{s=\uparrow,\downarrow}(c_{i\pi
s}^{\dagger}c_{j\pi s}+c_{j\pi s}^{\dagger}c_{i\pi s}),$ (51)

where $c_{i_{0}\alpha s}$ ($c_{i_{0}\alpha s}^{\dagger}$) is the annihilation
(creation) operator of an electron at the impurity site $\alpha\,(={\rm
A},{\rm B})$ with spin $s\,(=\uparrow,\downarrow)$, $n_{i_{0}\alpha
s}=c^{\dagger}_{i_{0}\alpha s}c_{i_{0}\alpha s}$, and
$n_{i_{0}\alpha}=n_{i_{0}\alpha\uparrow}+n_{i_{0}\alpha\downarrow}$ with
$i_{0}$ labeling the location of the impurity site B in the honeycomb lattice
as well as the impurity site A. $c_{i\pi s}$ ($c_{i\pi s}^{\dagger}$) is the
annihilation (creation) operator of an electron at conduction site $i$ with
spin $s\,(=\uparrow,\downarrow)$ and $c_{i_{0}\pi s}\equiv c_{i_{0}{\rm B}s}$.
${\bm{S}}_{i_{0}\alpha}=(S_{i_{0}\alpha}^{x},S_{i_{0}\alpha}^{y},S_{i_{0}\alpha}^{z})$
is the spin operator at the impurity site $\alpha$ given by

$S_{i_{0}\alpha}^{\nu}=\frac{1}{2}\sum_{s,s^{\prime}=\uparrow,\downarrow}c_{i_{0}\alpha
s}^{\dagger}[\mbox{\boldmath{$\sigma$}}^{\nu}]_{ss^{\prime}}c_{i_{0}\alpha
s^{\prime}}.\quad(\nu=x,y,z)$ (52)

The sums indicated by $\langle i,j\rangle$ and
$\mathopen{\hbox{\set@color${\langle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\langle}$}}i,j\mathclose{\hbox{\set@color${\rangle}$}\kern-1.94444pt\leavevmode\hbox{\set@color${\rangle}$}}$
run over all pairs of nearest and next nearest pairs of sites $i$ and $j$,
respectively, in the honeycomb lattice, including the impurity site B. Notice
that the fourth term in Eq. (51) is introduced to correct the double counting
of the interactions, as in the case of $\mathcal{H}_{\rm AM}$.

Using the Lanczos transformation for the single-particle hopping terms in the
honeycomb lattice, model I described by the Hamiltonian
$\mathcal{H}^{(\text{I})}$ is mapped onto the following Q1D model:

$\displaystyle\tilde{\mathcal{H}}^{(\text{I})}_{\rm Q1D}$
$\displaystyle=U\tilde{n}_{1\uparrow}\tilde{n}_{1\downarrow}-\frac{U}{2}(\tilde{n}_{1\uparrow}+\tilde{n}_{1\downarrow})$
$\displaystyle+V\sum_{s=\uparrow,\downarrow}(\tilde{c}_{1s}^{\dagger}\tilde{c}_{2s}+\tilde{c}_{2s}^{\dagger}\tilde{c}_{1s})$
$\displaystyle+\sum_{l=2}^{L}\sum_{s=\uparrow,\downarrow}\varepsilon_{l}\tilde{c}_{ls}^{\dagger}\tilde{c}_{ls}$
$\displaystyle+\sum_{l=2}^{L-1}\sum_{s=\uparrow,\downarrow}t_{l}(\tilde{c}_{ls}^{\dagger}\tilde{c}_{l+1s}+\tilde{c}_{l+1s}^{\dagger}\tilde{c}_{ls}),$
(53)

where $\tilde{c}_{ls}$ ($\tilde{c}_{ls}^{\dagger}$) is the electron creation
(annihilation) operator generated by the $(l-1)$th Lanczos iteration $(l\geq
2)$ with $\tilde{c}_{1s}\equiv c_{{\rm H}s}$ and $\tilde{c}_{2s}\equiv
c_{i_{0}\pi s}$, $\tilde{n}_{ls}=\tilde{c}_{ls}^{\dagger}\tilde{c}_{ls}$, and
the generalized Lanczos coefficients $\varepsilon_{l}$ and $t_{l}$ are
determined through the Lanczos procedure Shirakawa and Yunoki (2014). The
schematic depiction of model I in the Q1D representation described by the
Hamiltonian $\tilde{\mathcal{H}}^{(\text{I})}_{\rm Q1D}$ is shown in Fig.
15(c).

Similarly, model II described by the Hamiltonian $\mathcal{H}^{(\text{II})}$
is mapped onto the following Q1D model:

$\displaystyle\tilde{\mathcal{H}}^{(\text{II})}_{\rm Q1D}$
$\displaystyle=U\sum_{l=1,2}\tilde{n}_{l\uparrow}\tilde{n}_{l\downarrow}+U^{\prime}\tilde{n}_{1}\tilde{n}_{2}$
$\displaystyle-2J\left(\tilde{\bm{S}}_{1}\cdot\tilde{\bm{S}}_{2}+\frac{1}{4}\tilde{n}_{1}\tilde{n}_{2}\right)$
$\displaystyle+J^{\prime}\left(\tilde{c}_{1\uparrow}^{\dagger}\tilde{c}_{1\downarrow}^{\dagger}\tilde{c}_{2\downarrow}\tilde{c}_{2\uparrow}+\tilde{c}_{2\uparrow}^{\dagger}\tilde{c}_{2\downarrow}^{\dagger}\tilde{c}_{1\downarrow}\tilde{c}_{1\uparrow}\right)$
$\displaystyle-\frac{1}{2}(U+2U^{\prime}-J)\sum_{l=1,2}\tilde{n}_{l}$
$\displaystyle+\sum_{l=2}^{L}\sum_{s=\uparrow,\downarrow}\varepsilon_{l}\tilde{c}_{ls}^{\dagger}\tilde{c}_{ls}$
$\displaystyle+\sum_{l=2}^{L-1}\sum_{s=\uparrow,\downarrow}t_{l}(\tilde{c}_{ls}^{\dagger}\tilde{c}_{l+1s}+\tilde{c}_{l+1s}^{\dagger}\tilde{c}_{ls}),$
(54)

where $\tilde{n}_{l}=\tilde{n}_{l\uparrow}+\tilde{n}_{l\downarrow}$,
$\tilde{\bm{S}}_{l}$ is the spin operator given in Eq. (52) but for
$\tilde{c}_{ls}$ instead of $c_{i_{0}\alpha s}$, $\tilde{c}_{1s}\equiv
c_{i_{0}{\rm A}s}$, and $\tilde{c}_{2s}\equiv c_{i_{0}{\rm B}s}$. The
schematic depiction of model II in the Q1D representation described by the
Hamiltonian $\tilde{\mathcal{H}}^{(\text{II})}_{\rm Q1D}$ is shown in Fig.
15(d). Notice that the generated Lanczos coefficients $\varepsilon_{l}$ and
$t_{l}$ in Eq. (54) are exactly the same as those in Eq. (53). Therefore, the
difference between models I and II given in Eqs. (53) and (54) is a type of
couplings between sites $l=1$ and 2 represented by operators $\tilde{c}_{1s}$
and $\tilde{c}_{2s}$, respectively: these two sites are connected through the
hopping in model I while they are connected through the interactions in model
II, as also indicated in Figs. 15(c) and 15(d), besides the site at $l=2$
being one of the impurity sites in model II.

Model I is called a pseudogap Anderson model, in which the impurity site is
connected to the conduction band through a hybridization function with a
pseudogap structure at Fermi level, corresponding to the undoped case in model
I when the chemical potential $\mu$ is at $2t_{\pi}\,(\equiv\mu_{0})$. The
ground state of model I in the undoped case is in the local moment phase where
a free magnetic moment exists Shirakawa and Yunoki (2014). In the doped case
with $\mu\neq\mu_{0}$, the ground state is in either the Kondo screening phase
(i.e., the symmetric strong coupling phase) or the asymmetric strong coupling
phase, where no free magnetic moment appears Vojta  _et al._ (2010).

Figure 16:  Log-log plots of the absolute values of the spin correlation
function $\langle\tilde{{\bm{S}}}_{1}\cdot\tilde{\bm{S}}_{l}\rangle$ between
the impurity and conduction sites for (a) the doped case and (b) the undoped
case with various values of $L$ in the Q1D model $\tilde{\mathcal{H}}_{\rm
Q1D}^{\rm(I)}$. The number $N_{e}$ of electrons is set to be $N_{e}=L$ in (a),
corresponding to the doped case, while $N_{e}$ is adjusted so as to minimize
$E_{0}(N_{e})-\mu_{0}N_{e}$ with the chemical potential $\mu_{0}=0.25t_{\pi}$
in (b), corresponding to the undoped case, where $E_{0}(N_{e})$ is the ground
state energy of $\tilde{\mathcal{H}}_{\rm Q1D}^{\rm(I)}$ with $N_{e}$
electrons in the one-dimensional cluster of size $L$. The other parameters
used are $U=3t$, $V=t$, and $t_{\pi}=t/12$. Red, orange, green, cyan, and blue
symbols correspond to the results for $L=12$, $24$, $48$, $96$, and $192$,
respectively. Solid and dashed lines represent functions proportional to
$(l-1)^{-1}$ and $(l-1)^{-2}$, respectively, for comparison.

Figure 16 shows the spin correlation function
$\langle\tilde{\bm{S}}_{1}\cdot\tilde{\bm{S}}_{l}\rangle$ between the impurity
site ($l=1$) and the conduction sites ($l\geq 2$) in the Q1D representation of
the single-impurity Anderson model $\tilde{\mathcal{H}}^{(\text{I})}_{\rm
Q1D}$. We find in Fig. 16(a) that the spin correlation function
$\langle\tilde{\bm{S}}_{1}\cdot\tilde{\bm{S}}_{l}\rangle$ for the doped case
decays faster than $d^{-1}$, where $d=l-1$ is the ditance between the impurity
site and the conduction site labeled by $l$ ($l\geq 2$) in the Q1D model
$\tilde{\mathcal{H}}^{(\text{I})}_{\rm Q1D}$. In contrast,
$\langle\tilde{\bm{S}}_{1}\cdot\tilde{\bm{S}}_{l}\rangle$ for the undoped case
in Fig. 16(b) shows somewhat nonmonotonic behavior as a function of the
distance $d$. If we focus on the distances near the impurity site in the
largest system (denoted by blue squares),
$\langle\tilde{\bm{S}}_{1}\cdot\tilde{\bm{S}}_{l}\rangle$ decays approximately
as $d^{-1}$. If we focus on the spin correlation function at the maximum
distance of a given system size $L$, i.e.,
$\langle\tilde{\bm{S}}_{1}\cdot\tilde{\bm{S}}_{L}\rangle$, the envelope seems
to decay slower that $d^{-1}$. Although, it is not easy to determine the
exponent precisely in Fig. 16(b), we can safely claim the much slower decay of
$\langle\tilde{\bm{S}}_{1}\cdot\tilde{\bm{S}}_{l}\rangle$ for the undoped case
than for the doped case.

Figure 17:  Log-log plots of the absolute values of the spin correlation
functions (a,b) $\langle\tilde{{\bm{S}}}_{1}\cdot\tilde{\bm{S}}_{l}\rangle$
and (c,d) $\langle\tilde{{\bm{S}}}_{2}\cdot\tilde{\bm{S}}_{l}\rangle$ between
the impurity and conduction sites for the doped case (left panels) and the
undoped case (right panels) with various values of $L$ in the Q1D model
$\tilde{\mathcal{H}}_{\rm Q1D}^{\rm(II)}$. The number $N_{e}$ of electrons is
set to be $N_{e}=L$ in (a,c), corresponding to the doped case, while $N_{e}$
is adjusted so as to minimize $E_{0}(N_{e})-\mu_{0}N_{e}$ with the chemical
potential $\mu_{0}=0.25t_{\pi}$ in (b,d), corresponding to the undoped case,
where $E_{0}(N_{e})$ is the ground state energy of $\tilde{\mathcal{H}}_{\rm
Q1D}^{\rm(II)}$ with $N_{e}$ electrons in the one-dimensional cluster of size
$L$. The other parameters used are $U=3t$, $J=0.1U$, and $t_{\pi}=t/12$. Red,
orange, green, cyan, and blue symbols correspond to the results for $L=12$,
$24$, $48$, $96$, and $192$, respectively. Solid and dashed lines represent
functions proportional to $(l-2)^{-1}$ and $(l-2)^{-2}$, respectively, for
comparison.

In model II, the local spin triplet is favored at the impurity sites at $l=1$
and $2$ because of the Hund’s coupling in the strong coupling phase (e.g.,
with $U=3t$, $J=0.1U$, and $t_{\pi}=t/12$) for the undoped case, and this
spin-1 impurity is coupled to the conduction band via the hybridization
function that diverges at Fermi level Shirakawa and Yunoki (2014). The Fermi
level is shifted away from this diverging point of the hybridization function
in the doped system. However, in both undoped and doped cases, the impurity
spin is similarly screened partially by the conduction electrons and there
appears residual unscreened free magnetic moment with its spin as large as 1/2
Mattis (1967); Nozières, Ph. and Blandin, A. (1980).

Figures 17(a) and 17(b) show the spin correlation function
$\langle\tilde{\bm{S}}_{1}\cdot\tilde{\bm{S}}_{l}\rangle$ between the impurity
site ($l=1$) and the conduction sites ($l\geq 3)$. We find that in both
undoped and doped cases
$\langle\tilde{\bm{S}}_{1}\cdot\tilde{\bm{S}}_{l}\rangle$ decays rather close
to $d^{-1}$ with the distance $d=l-2$ and shows a similar behavior to that for
the undoped case in model I [see Fig. 16(b)], where the free magnetic moment
appears. In contrast, as shown in Figs. 17(c) and 17(d), the spin correlation
function $\langle\tilde{\bm{S}}_{2}\cdot\tilde{\bm{S}}_{l}\rangle$ from the
second impurity site ($l=2$), which is directly connected to the conduction
sites, decays approximately as $d^{-2}$, specially focusing on the region
close to the impurity site for the larger systems. This $d^{-2}$ decay is
expected for the Fermi liquid and thus it is consistent with the partial
screening of the impurity spin at $l=2$.

When we compare these results with those for the Anderson model
$\mathcal{H}_{\rm AM}$ in Sec. IV.3, we have to remind that the accessible
system size is very limited for the Anderson model $\mathcal{H}_{\rm AM}$. For
example, the nearly $d^{-2}$ decay of the spin correlation functions
$\langle\tilde{\bm{S}}_{1}\cdot\tilde{\bm{S}}_{l}\rangle$ in Fig. 16(a) and
$\langle\tilde{\bm{S}}_{2}\cdot\tilde{\bm{S}}_{l}\rangle$ in Figs. 17(c) and
17(d) can still be observed even in the moderate systems as large as $L=48$,
which is the largest accessible system size for the Anderson model
$\mathcal{H}_{\rm AM}$ in Sec. IV.3. Indeed, we can observe the spin
correlation function consistent with the $d^{-2}$ behavior in Figs. 14(a),
14(c), and 14(d). On the other hand, the spin correlation function
$\langle\tilde{\bm{S}}_{1}\cdot\tilde{\bm{S}}_{l}\rangle$ for
$\mathcal{H}^{({\text{II}})}_{\rm Q1D}$ with $L=48$, denoted by green symbols
in Figs. 17(a) and 17(b), appears to be constant rather than decaying
algebraically with some power, which is similar to the behavior of the spin
correlation function $\langle{\bm{S}}_{\sigma}\cdot\tilde{\bm{S}}_{lc}\rangle$
($c=a,e$) in Fig. 13.

The difference between the doped and undoped cases in these systems is most
difficult to distinguish, since the exponent for the undoped case deviates
from not only $d^{-1}$ but also $d^{-2}$ even in the largest system of model
I, as shown in Fig. 16(b). Rather, the exponent seems to vary between $d^{-1}$
and $d^{-2}$ for different system sizes. Therefore, the clear determination of
the decay exponent is difficult even in the simplest model such as model I,
based on the results available at this moment. However, we can assert that the
nontrivial decay exponent slower than $d^{-2}$ is the hallmark of the presence
of the free magnetic moment Mitchell  _et al._ (2011).

## References

  * Yazyev and Helm (2007) Oleg V. Yazyev and Lothar Helm, “Defect-induced magnetism in graphene,” Phys. Rev. B 75, 125408 (2007).
  * Dharma-wardana and Zgierski (2008) M.W.C. Dharma-wardana and Marek Z. Zgierski, “Magnetism and structure at vacant lattice sites in graphene,” Physica E: Low-dimensional Systems and Nanostructures 41, 80–83 (2008).
  * Dai  _et al._ (2011) X. Q. Dai, J. H. Zhao, M. H. Xie, Y. N. Tang, Y. H. Li, and B. Zhao, “First-principle study of magnetism induced by vacancies in graphene,” The European Physical Journal B 80, 343–349 (2011).
  * Paz  _et al._ (2013) Wendel S. Paz, Wanderlã L. Scopel, and Jair C.C. Freitas, “On the connection between structural distortion and magnetism in graphene with a single vacancy,” Solid State Communications 175-176, 71–75 (2013), special Issue: Graphene V: Recent Advances in Studies of Graphene and Graphene analogues.
  * Padmanabhan and Nanda (2016) Haricharan Padmanabhan and B. R. K. Nanda, “Intertwined lattice deformation and magnetism in monovacancy graphene,” Phys. Rev. B 93, 165403 (2016).
  * Valencia and Caldas (2017) A. M. Valencia and M. J. Caldas, “Single vacancy defect in graphene: Insights into its magnetic properties from theoretical modeling,” Phys. Rev. B 96, 125431 (2017).
  * El-Barbary  _et al._ (2003) A. A. El-Barbary, R. H. Telling, C. P. Ewels, M. I. Heggie, and P. R. Briddon, “Structure and energetics of the vacancy in graphite,” Phys. Rev. B 68, 144107 (2003).
  * Lehtinen  _et al._ (2004) P. O. Lehtinen, A. S. Foster, Yuchen Ma, A. V. Krasheninnikov, and R. M. Nieminen, “Irradiation-induced magnetism in graphite: A density functional study,” Phys. Rev. Lett. 93, 187202 (2004).
  * Ney  _et al._ (2011) Andreas Ney, Pagona Papakonstantinou, Ajay Kumar, Nai-Gui Shang, and Nianhua Peng, “Irradiation enhanced paramagnetism on graphene nanoflakes,” Applied Physics Letters 99, 102504 (2011).
  * Nair  _et al._ (2012) R. R. Nair, M. Sepioni, I-Ling Tsai, O. Lehtinen, J. Keinonen, A. V. Krasheninnikov, T. Thomson, A. K. Geim, and I. V. Grigorieva, “Spin-half paramagnetism in graphene induced by point defects,” Nature Physics 8, 199–202 (2012).
  * Nair  _et al._ (2013) 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, “Dual origin of defect magnetism in graphene and its reversible switching by molecular doping,” Nature Communications 4, 2010 (2013).
  * Sepioni  _et al._ (2010) M. Sepioni, R. R. Nair, S. Rablen, J. Narayanan, F. Tuna, R. Winpenny, A. K. Geim, and I. V. Grigorieva, “Limits on intrinsic magnetism in graphene,” Phys. Rev. Lett. 105, 207205 (2010).
  * Zhang  _et al._ (2016) Yu Zhang, Si-Yu Li, Huaqing Huang, Wen-Tian Li, Jia-Bin Qiao, Wen-Xiao Wang, Long-Jing Yin, Ke-Ke Bai, Wenhui Duan, and Lin He, “Scanning tunneling microscopy of the $\pi$ magnetism of a single carbon vacancy in graphene,” Phys. Rev. Lett. 117, 166801 (2016).
  * Casartelli  _et al._ (2013) M. Casartelli, S. Casolo, G. F. Tantardini, and R. Martinazzo, “Spin coupling around a carbon atom vacancy in graphene,” Phys. Rev. B 88, 195424 (2013).
  * Jr  _et al._ (2019) Max Pinheiro Jr, Daniely V. V. Cardoso, Adélia J. A. Aquino, Francisco B. C. Machado, and Hans Lischka, “The characterization of electronic defect states of single and double carbon vacancies in graphene sheets using molecular density functional theory,” Molecular Physics 117, 1519–1531 (2019).
  * González-Herrero  _et al._ (2016) Héctor González-Herrero, José M. Gómez-Rodríguez, Pierre Mallet, Mohamed Moaied, Juan José Palacios, Carlos Salgado, Miguel M. Ugeda, Jean-Yves Veuillen, Félix Yndurain, and Iván Brihuega, “Atomic-scale control of graphene magnetism by using hydrogen atoms,” Science 352, 437–441 (2016).
  * Kanao  _et al._ (2012) Taro Kanao, Hiroyasu Matsuura, and Masao Ogata, “Theory of defect-induced kondo effect in graphene: Numerical renormalization group study,” Journal of the Physical Society of Japan 81, 063709 (2012).
  * Cazalilla  _et al._ (2012) M. A. Cazalilla, A. Iucci, F. Guinea, and A. H. Castro Neto, “Local moment formation and kondo effect in defective graphene,” (2012), arXiv:1207.3135 [cond-mat.str-el] .
  * Mitchell and Fritz (2013) Andrew K. Mitchell and Lars Fritz, “Kondo effect with diverging hybridization: Possible realization in graphene with vacancies,” Phys. Rev. B 88, 075104 (2013).
  * Haldane (1983) F. D. M. Haldane, “Nonlinear field theory of large-spin heisenberg antiferromagnets: Semiclassically quantized solitons of the one-dimensional easy-axis néel state,” Phys. Rev. Lett. 50, 1153–1156 (1983).
  * Affleck  _et al._ (1987) Ian Affleck, Tom Kennedy, Elliott H. Lieb, and Hal Tasaki, “Rigorous results on valence-bond ground states in antiferromagnets,” Phys. Rev. Lett. 59, 799–802 (1987).
  * White (1992) Steven R. White, “Density matrix formulation for quantum renormalization groups,” Phys. Rev. Lett. 69, 2863–2866 (1992).
  * White (1993) Steven R. White, “Density-matrix algorithms for quantum renormalization groups,” Phys. Rev. B 48, 10345–10356 (1993).
  * Shirakawa and Yunoki (2014) Tomonori Shirakawa and Seiji Yunoki, “Block lanczos density-matrix renormalization group method for general anderson impurity models: Application to magnetic impurity problems in graphene,” Phys. Rev. B 90, 195109 (2014).
  * Allerdt  _et al._ (2015) Andrew Allerdt, C. A. Büsser, G. B. Martins, and A. E. Feiguin, “Kondo versus indirect exchange: Role of lattice and actual range of rkky interactions in real materials,” Phys. Rev. B 91, 085101 (2015).
  * Allerdt  _et al._ (2017) A. Allerdt, A. E. Feiguin, and S. Das Sarma, “Competition between kondo effect and rkky physics in graphene magnetism,” Phys. Rev. B 95, 104402 (2017).
  * Allerdt and Feiguin (2019) Andrew Allerdt and Adrian E. Feiguin, “A numerically exact approach to quantum impurity problems in realistic lattice geometries,” Frontiers in Physics 7 (2019), 10.3389/fphy.2019.00067.
  * Novoselov  _et al._ (2005) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, “Two-dimensional gas of massless dirac fermions in graphene,” Nature 438, 197–200 (2005).
  * Lieb (1989) Elliott H. Lieb, “Two theorems on the hubbard model,” Phys. Rev. Lett. 62, 1201–1204 (1989).
  * Pereira  _et al._ (2006) Vitor M. Pereira, F. Guinea, J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, “Disorder induced localized states in graphene,” Phys. Rev. Lett. 96, 036801 (2006).
  * Peres  _et al._ (2009) N. M. R. Peres, Shan-Wen Tsai, J. E. Santos, and R. M. Ribeiro, “Scanning tunneling microscopy currents on locally disordered graphene,” Phys. Rev. B 79, 155442 (2009).
  * Ugeda  _et al._ (2010) M. M. Ugeda, I. Brihuega, F. Guinea, and J. M. Gómez-Rodríguez, “Missing atom as a source of carbon magnetism,” Phys. Rev. Lett. 104, 096804 (2010).
  * Otsuka  _et al._ (2016) Yuichi Otsuka, Seiji Yunoki, and Sandro Sorella, “Universal quantum criticality in the metal-insulator transition of two-dimensional interacting dirac electrons,” Phys. Rev. X 6, 011029 (2016).
  * Seki  _et al._ (2019) Kazuhiro Seki, Yuichi Otsuka, Seiji Yunoki, and Sandro Sorella, “Fermi-liquid ground state of interacting dirac fermions in two dimensions,” Phys. Rev. B 99, 125145 (2019).
  * Kanamori (1963) Junjiro Kanamori, “Electron Correlation and Ferromagnetism of Transition Metals,” Progress of Theoretical Physics 30, 275–289 (1963).
  * Kretinin  _et al._ (2013) A. Kretinin, G. L. Yu, R. Jalil, Y. Cao, F. Withers, A. Mishchenko, M. I. Katsnelson, K. S. Novoselov, A. K. Geim, and F. Guinea, “Quantum capacitance measurements of electron-hole asymmetry and next-nearest-neighbor hopping in graphene,” Phys. Rev. B 88, 165427 (2013).
  * Wehling  _et al._ (2011) T. O. Wehling, E. Şaşıoğlu, C. Friedrich, A. I. Lichtenstein, M. I. Katsnelson, and S. Blügel, “Strength of effective coulomb interactions in graphene and graphite,” Phys. Rev. Lett. 106, 236805 (2011).
  * Bulla  _et al._ (1997) R Bulla, Th Pruschke, and A C Hewson, “Anderson impurity in pseudo-gap fermi systems,” Journal of Physics: Condensed Matter 9, 10463–10474 (1997).
  * Fritz and Vojta (2013) Lars Fritz and Matthias Vojta, “The physics of kondo impurities in graphene,” Reports on Progress in Physics 76, 032501 (2013).
  * Shirakawa and Yunoki (2016) Tomonori Shirakawa and Seiji Yunoki, “Density matrix renormalization group study in energy space for a single-impurity anderson model and an impurity quantum phase transition,” Phys. Rev. B 93, 205124 (2016).
  * Gonzalez-Buxton and Ingersent (1998) Carlos Gonzalez-Buxton and Kevin Ingersent, “Renormalization-group study of anderson and kondo impurities in gapless fermi systems,” Phys. Rev. B 57, 14254–14293 (1998).
  * Vojta  _et al._ (2010) M. Vojta, L. Fritz, and R. Bulla, “Gate-controlled kondo screening in graphene: Quantum criticality and electron-hole asymmetry,” EPL (Europhysics Letters) 90, 27006 (2010).
  * Bersuker (2001) Isaac B. Bersuker, “Modern aspects of the jahn-teller effect theory and applications to molecular problems,” Chemical Reviews 101, 1067 (2001).
  * Hallberg (1995) Karen A. Hallberg, “Density-matrix algorithm for the calculation of dynamical properties of low-dimensional systems,” Phys. Rev. B 52, R9827–R9830 (1995).
  * Mattis (1967) D. C. Mattis, “Symmetry of ground state in a dilute magnetic metal alloy,” Phys. Rev. Lett. 19, 1478–1481 (1967).
  * Nozières, Ph. and Blandin, A. (1980) Nozières, Ph. and Blandin, A., “Kondo effect in real metals,” J. Phys. France 41, 193–211 (1980).
  * Mitchell  _et al._ (2011) Andrew K. Mitchell, Michael Becker, and Ralf Bulla, “Real-space renormalization group flow in quantum impurity systems: Local moment formation and the kondo screening cloud,” Phys. Rev. B 84, 115120 (2011).