PHYSICAL REVIEW B 92, 155140 (2015)

Emergence of a Kondo singlet state with Kondo temperature well beyond 1000 K
in a proton-embedded electron gas
Yasutami Takada,1 Ryo Maezono,2 and Kanako Yoshizawa1
1

Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581, Japan
2
School of Information Science, JAIST, 1-1 Asahidai, Nomi, Ishikawa 923-1292, Japan(Received 31 October 2014; revised manuscript received 9 July 2015; published 23 October 2015)
Hydrogen in metals has attracted much attention for a long time from both basic scientific and technologicalpoints of view. Its electronic state has been investigated in terms of a proton embedded in the electron gas mostly
by the local density approximation (LDA) to the density functional theory. At high electronic densities, it iswell described by a bare proton H+ screened by metallic electrons (charge resonance), while at low densities
two electrons are localized at the proton site to form a closed-shell negative ion H− protected from surrounding
metallic electrons by the Pauli exclusion principle. However, no details are known about the transition from H+
to H− in the intermediate-density region. Here, by accurately determining the ground-state electron distributionn(r) by the use of LDA and diffusion Monte Carlo simulations with the total electron number up to 170,
we obtain a complete picture of the transition, in particular, a sharp transition from short-range H+ screening
charge resonance to long-range Kondo-type spin-singlet resonance, the emergence of which is confirmed by thepresence of an anomalous Friedel oscillation characteristic to the Kondo singlet state with the Kondo temperature
TK well beyond 1000 K. This study not only reveals interesting competition between charge and spin resonances,enriching the century-old paradigm of metallic screening to a point charge, but also discovers a high-TK systemlong sought in relation to the development of exotic superconductivity in the quantum critical regime.
DOI: 10.1103/PhysRevB.92.155140

PACS number(s): 71.10.Ca, 71.15.Mb, 71.45.Gm, 75.20.Hr

I. INTRODUCTION

Physics in the heavy-fermion superconductors has been
understood by the concept of quantum criticality in a systemof regularly arrayed dense Kondo impurities (Kondo lattice)[1–3] and the spin-fluctuation mechanism is believed to be
responsible for superconductivity, as inferred from the strong
correlation between the superconducting transition temperature Tc and the Kondo temperature TK [4–7]. More specifically,
Tc is of the order of 0.1TK , as shown in Fig. 1 plotted basedon the information available in the literature, from which we
can conceive an idea that high-Tc superconductivity will beobtained if we can discover a Kondo system with very high
TK . In fact, the recently discovered plutonium compounds suchas PuCoGa5 with Tc = 18.5 K and TK ≈ 260 K [8,9] may be
regarded as the successful realization of this idea. Thus, weshould make further pursuit of this idea by searching for a
new class of Kondo systems with TK higher than 1000 K.
Theoretically, this search can be done by the first-principles
quantitative determination of TK for the composite system of
an impurity atom embedded in a metal.
As first suggested by Debye and Hückel [10], an atomic
nucleus charge +Ze in a metal is screened by accumulationof metallic electrons which is regarded as a charge resonanceand well described by the linear-response theory. This conceptof metallic screening is prevailing for a century, but because
Z is not infinitesimally small, we need to consider nonlinear
effects in the screening, including the contribution from spin
fluctuations. The spin contribution will be enhanced, if Z isan odd integer, such as Z = 1 (case of a proton), in which aspin-polarized bound state might appear at the impurity atom.
With the above basic scientific issues in mind, we haveconcerned with the problem of hydrogen impurity in metals
which attracts long attention from a technological point of view[11], such as hydrogen storage in solids, sensor applications,and catalysis. Its electronic state has been investigated in
1098-0121/2015/92(15)/155140(11)terms of a proton immersed into an interacting many-electron
system plus a compensating background (electron gas: EG)
since 1970s because this is an ideal system to study important
topics related to an impurity in metals, such as the embedding
energy, a key quantity in the effective-medium theory [12–14].
It is also studied from a motivation to improve on the local
density approximation (LDA) to the density functional theory(DFT) by using the electron distribution n(r) obtained byquantum Monte Carlo methods [15–17]. Nonlinear metallic
screening is another extensively examined topic in this system[18–20], but no serious attention has been paid so far to thespin resonance effect.
The homogeneous EG with the average density n0 is
specified by a single parameter rs , defined by rs = (3/4π n0 )1/3
in units of the Bohr radius aB . (We use atomic units hereafter.)
Its characteristic energy is the Fermi energy εF , given by
kF2 /2 (= 1.84rs−2 ) with kF (= 1.92rs−1 ) the Fermi momentum[Fig. 2(a)]. Hydrogen, on the other hand, has two typicalenergies [Fig. 2(b)], the 1s level ε1s (= 0.5) and the electron
affinity εA (= 0.0278). Ratio of εF to ε1s or εA determines thecharacter of the ground state in the proton-embedded EG; forhigh n0 corresponding to εF  ε1s (or rs  2), the 1s levelis buried in the continuum of EG and thus no electrons arebound to a bare proton H+ , leading to a charge resonance (CR)state in which H+ is dielectrically screened by accumulation ofitinerant electrons near the Fermi level μ [Fig. 2(c)]. For low n0
with εF  εA (or rs  10), H+ captures two antiparallel-spin
electrons to form H− . This closed-shell negative ion resides in
EG with repelling other electrons owing to the Pauli exclusion
principle [Fig. 2(d)], but if εF increases and reaches as high as
εA , the Fermi pressure from EG to the ion becomes so largethat the localized electrons in H− begin to spill out into EG.
Then, a crucial question is whether this state at εF ≈ εA is thesame as that in Fig. 2(c) or not.
Intuitively, for εA  εF  ε1s , we can imagine a spinpolarized state made of a single electron with either up or

155140-1©2015 American Physical Society

YASUTAMI TAKADA, RYO MAEZONO, AND KANAKO YOSHIZAWA

50.0

Tc [K]The Hamiltonian for IAM, HA , is written as [22]

 †
†
†
†
HA =
εk ckσ ckσ + Ed
cdσ cdσ + U cd↑ cd↑ cd↓ cd↓

Tc = 0.1TK
PuCoGa5

10.0
5.0

U6 Fe

kσ

PuRhGa5+

CeCoIn5
UPd2 Al3
CeRhIn5
UBe13
URu 2 Si 2
CeCu 2Si2
CeCu 2Ge2
UPt3
CePd2 Si2 CeIrIn CeRh2 Si2

1.0
0.5

CeNi2 Ge2

5

10

T K [K]50 100500

FIG. 1. (Color online) Superconducting transition temperature
Tc versus Kondo temperature TK (a characteristic energy scale forspin fluctuations) in heavy-fermion superconductors.

down spin captured by H+ , but in view of the concept ofspin screening to form a Kondo singlet [1,21] in the impurity
Anderson model (IAM) [22], we anticipate the emergence of
not a spin-polarized but a Kondo-type spin-singlet resonance(SSR) state SSR [Fig. 2(e)]. Because there is no clear
distinction between conduction and localized electrons, thisSSR state is composed of only itinerant electrons near μ
without a local spin moment, similar to CR, but an importantdifference exists in the screening length; for CR, it is theThomas-Fermi length λTF ≈ kF−1 , but for SSR, the Kondoscreening length ξK is much longer than kF−1 , leading to ananomalous Friedel oscillation [23–26]. Then, the main aim of
this paper is to confirm this conjecture about the emergence of
SSR in the proton-embedded electron gas with determining TK
from first principles, but this confirmation is not an easy task
due to the existence of various difficulties, as we shall explainin the following in some detail.

1s (d)

F

A

↑

↑

+

↑

ε1s or1s(e)

ε1s +εA

ε1s
1s

↑

F

A

1s

H

H+

H Affinity:
εA = 0.0278

kwhere 0 is the Slater determinant made of conductionelectron orbitals and the parameters a0 and ak are to be
determined variationally. As Eq. (2) clearly shows, SSR is
not described by a single Slater determinant but is a correlated
many-body state. It must also be noted that continuum
conduction states around the Fermi level μ are indispensablefor the construction of this SSR .
In our work, we are not allowed to employ this simple
model HA . Instead, we have to start with the first-principles
Hamiltonian H , described in atomic units as
H =−

 ∇2
i

i

2

+

+


1
1
vext (r i ) + CNZ , (3)+2 i=j |r i − r j |i

in first quantization. Here, we have considered a neutral atomof atomic number Z at the origin of coordinates immersed
into the jellium sphere of radius R and average density n0 .
The number of electrons contained in the jellium sphere is
4π R 3 n0 /3 = (R/rs )3 , so that the total electron number N isequal to Z + (R/rs )3 , satisfying global neutrality, from which
we obtain R = (N − Z)1/3 rs . In Eq. (3), the external potential
working on an electron vext (r) is composed of the potentialfrom the nucleus and that from the positive background,
written as
N − Z 3R 2 − r 2Z−θ (R − |r|)
|r|
2
R3
N −Z
−
θ (|r| − R),
|r|

H

H
↑

(1)

vext (r) = −↑F(b) Hydrogen (H)H 1s level:
ε 1s= 0.5↑(c)Fermi energy: εF =1.84/rs 2
εk
μ
εF
k
kF

σ
†
†
Vd k (cdσ ckσ + ckσ cdσ ),in second quantization with use of the annihilation operator
ckσ for a conduction electron with wave vector k, spin σ , andone-body band energy εk , while cdσ is an operator to destroy
a localized d electron with spin σ at the energy level Ed . Thed electrons not only interact to each other at the localized sitewith the strength U , but also hybridize with the conduction
electrons with the strength Vd k . Based on HA , SSR is given
as [27,28]
††SSR = a0 0 +ak (cd↑ ck↑ + cd↓ ck↓ )0 ,(2)↑(a) Electron gas (EG)


kσ

5

0.1
1

PHYSICAL REVIEW B 92, 155140 (2015)

FIG. 2. (Color online) Proton-embedded EG with characteristicenergies in (a) and (b). Three possible ground states are schematicallyillustrated in (c)–(e) corresponding, respectively, to a bare proton H+
screened by metallic electrons with the screening length ≈kF−1 (a CR
state), a closed-shell ion H− confined in EG, and an SSR state in which
if the 1s (virtual) level of H is temporarily occupied by a single up(down-) spin electron, down- (up-) spin clouds of itinerant electrons
are formed around H for spin screening with the screening lengthξK  kF−1 . These temporary states are superposed with interchanging
the roles of spins to make a resonance state SSR .

(4)

with θ (x) the Heaviside function and CNZ represents
the Coulomb self-energy stemming from both nucleusbackground and intrabackground interactions, given by
CNZ =3 Z(N − Z)2/3
3 (N − Z)5/3
+
.
2
rs
5
rs

(5)

In solving Eq. (3), we impose the fixed boundary condition
to make the wave function vanish at |r i | = R. Irrespective ofwhether we include the constant term CNZ in H or not, there
is no problem of divergence in this finite-N system, but CNZ isneeded in order to achieve the mutual global cancellation inenergies between the electron-background attractive potentials

155140-2

EMERGENCE OF A KONDO SINGLET STATE WITH KONDO . . .

and the repulsive potentials of both electron-electron and intrabackground interactions for the bulk (N → ∞) system [29].
By comparing H in Eq. (3) with HA , we readily see theirdifferences; first, there is no predetermined localized statein H and thus a (virtual) localized state, if any, must be
determined in the first place in constructing SSR with relatedparameters such as Ed , U , and Vd k , if necessary. Second andmore importantly, the long-range Coulomb interaction, whichis completely neglected in HA , works among metallic electronsin H , making solution of the problem quite difficult. In fact,even in the homogeneous EG without the embedded atom, theproblem is sufficiently complicated in both variational [30–32]
and Green’s-function [33,34] approaches. Incidentally, thislong-range Coulomb interaction is the source to bring about
CR and therefore it is indispensable for discussing competition
between CR and SSR. Since this discussion constitutes anotherimportant aim of this paper, we can never neglect this longrange Coulomb interaction, making all theoretical and computational techniques developed so far for HA useless to H .
Usually, the first-principles Hamiltonian is solved by either
diffusion Monte Carlo (DMC) simulations or DFT-based methods. The former is an excellent method to obtain fairly accurate
results for the ground state, but it can never directly treat
SSR , because DMC simulations can be done only for finite-N
systems in which all levels are discrete, while in constructing
SSR , we need continuum conduction states which are allowed
only in the bulk (N → ∞) system. The latter methods caneasily treat the bulk system, but the ground-state physicalquantities are calculated in terms of a single Slater determinantmade of Kohn-Sham (KS) orbitals introduced in DFT, so that it
is not clear at all as to how much the obtained quantities reflect
the highly correlated many-body nature of SSR and how
accurate they are, especially because in actual calculations wealways have to resort to some approximation to the exchangecorrelation energy functional E xc [n(r)] such as LDA.
Faced with those difficulties, we have decided to focus
on n(r) rather than the wave function SSR itself, mainly
because DFT can, in principle, provide exact n(r) andthe corresponding ground-state energy E0 by projecting thereal interacting many-body system to a fictitious auxiliary
noninteracting system in which n(r) can be calculated with useof a single Slater determinant made of KS orbitals, even if weknow nothing about SSR in the real system. This nontrivial
assertion, one of the central theorems in DFT, is rigorously
proved as long as the ground state is nondegenerate [35], as is
the case for SSR. Of course, information obtained only throughn(r) and E0 is limited and useless for discussing transport [36]and excited-state properties, but we claim that it is still plenty
enough for our purpose of distinguishing between CR andSSR states and determining TK in SSR. In the context of DFT,
given exact n(r), it is also an interesting issue to clarify how
KS orbitals, which are defined in the fictitious system and thus
have basically no physical relevance, behave so as to correctly
provide n(r) in a strongly correlated state such as SSR. Thus,
this clarification constitutes an additional aim of this paper.
In implementing calculations of n(r), we adopt the following strategies: (i) In actual LDA calculations, we employ
the local spin-density approximation (LSDA) [37] by choosing a spin-dependent exchange-correlation energy functional
E xc [nσ (r)] written in terms of the spin-resolved electron

PHYSICAL REVIEW B 92, 155140 (2015)

distribution nσ (r), so that we can check a possibility of thespin-polarized (i.e., non-spin-singlet) ground state by detecting
the difference between n↑ (r) and n↓ (r). (ii) We assess theresults in LSDA at finite-N systems in comparison withthose in fixed-node DMC [38] with taking N up to 170, amuch larger size than those in previous variational Monte
Carlo (VMC) calculations [16,17]. (iii) We obtain n(r) inthe bulk system by the calculation in LSDA and also by
extrapolation of DMC data to N → ∞. The extrapolated
results for n(r) in DMC are independent of N and thus we canassume that they will be free from any restrictions incurredfrom the fixed-node approximation, the only approximation
adopted in DMC simulations, because the fixed-node positionsin DMC are prescribed by R [and consequently by N for
given rs due to R = (N − Z)1/3 rs ] in the fixed boundary
condition, indicating that independence of N also suggests
independence of the postulated node positions. (iv) We check
whether the obtained n(r) at N → ∞ exhibits the behavior
characteristic to the Kondo SSR state or not. More specifically,
we look for modulation of the Friedel-oscillation period, avery important inherent property of the anomalous Friedel
oscillation, to which we call serious attention for the first timein quantitatively determining ξK and consequently TK fromfirst principles.
In accordance with those strategies, we have investigatedn(r) in finite-N systems in both LSDA and DMC to find spinunpolarized ground states with a strong size effect and a seriesof magic numbers (10, 60, 170, . . .) of N at which convergent
results are easily obtained. These features can be explainedin terms of the emergence of SSR with its long-range natureof ξK . Its emergence is also signaled in LSDA at N → ∞for rs  2 by the appearance of a strange shallow bound KS
orbital below the conduction band bottom with an unusuallylong binding radius. In spite of the large size effect, for eachN , n(r) in LSDA agrees very well with that in DMC.
By summarizing the results thus calculated, we haveobtained a ground-state diagram in (rs , N −1 ) space, shownin Fig. 3, exhibiting sharp sequential transitions among CR,
SSR, and closed-shell H− ion confinement states. Contrary to
the previous explanation [12], we claim that the very shallow
bound KS orbital found in LSDA at intermediate densities is
not a physical H− but appears just to describe the long-rangechange of n(r) over ξK in SSR in the form of an envelope
density. We also find that hydrogen is most stably embeddedin EG in the SSR region, especially, optimally firmly at rs ≈ 4with the Kondo temperature TK ≈ 2100 K, indicating our
success in discovering a long-sought high-TK system.
In Sec. II, we explain the calculation methods in both LSDA
and DMC. In Sec. III, we show the calculated results and in
Sec. IV we discuss on the obtained results, together with theirimplications and future directions. Finally, in Sec. V, we givea summary of this paper.

II. CALCULATION METHODS
A. LSDA in the finite-N system

Let us consider the neutral system of a single nucleus withatomic number Z embedded in the N -electron jellium sphereof radius R. Its Hamiltonian H is given in Eq. (3). In LSDA

155140-3

YASUTAMI TAKADA, RYO MAEZONO, AND KANAKO YOSHIZAWA

1/N

B. VMC

Isolated H - ion

1/2

With use of the lowest-N KS orbitals thus obtained, we candefine the Slater determinant 0 (r 1 , . . . ,r N ), with which thetrial many-body ground-state wave function (r 1 , . . . ,r N ) for
the VMC calculation can be constructed in the Slater-Jastrow
type as [40]

1/10

1/60
1/170
State c:
H + induced
chargeresonance

1/∞

PHYSICAL REVIEW B 92, 155140 (2015)

State e:Spin-singlet
resonance
(SSR)

2

0

4

6

rs

8

10

State d:Closed-shellH - ionconfinement12

(r 1 , . . . ,r N ) = exp[J (r 1 , . . . ,r N )]0 (r 1 , . . . ,r N ), (12)14FIG. 3. (Color online) Ground-state diagram in the protonembedded electron gas in (rs , N −1 ) space with N the total electron
number, indicating sharp but size-dependent sequential transitions
among CR, SSR, and closed-shell ion H− . Rigorously speaking,
SSR is defined only at N → ∞, but the states for finite N directly
connected to SSR at N → ∞ are also called SSR. In the bulk
system, the screening length ξK is predicted to diverge at the CR-SSR
boundary or rs ≈ 1.97 with the change of rs , signaling the sharptransition. For an isolated H− ion at N = 2, we have used the exact
data for n(r) [39] to determine rs = 5.04 by averaging the local
rs (r){= [3/4π n(r)]1/3 } over the weight of n(r) itself.

to DFT, the KS equation is written as−∇ 2 /2 + vσKS (r) φiσ (r) = εiσ φiσ (r),where vσxc (r; [nσ ]) is derived from E xc [nσ ] through thefunctional derivative as(8)With use of the lowest-Nσ KS orbitals, nσ (r) is given bynσ (r) =Nσ


|φiσ (r)|2 ,

i

+

(9)

(10)

at |r| = R = (N − Z)1/3 rs . By using those converged quantities, we can calculate E0 (N,Z) the ground-state energy
including the constant term CNZ byE0 (N,Z) =εiσ +d r vext (r) − vσKS (r) nσ (r)
+

1
2



σ

dr dr

u3 (r i ,r j ,r i − r j ).

(13)

The actual choice of the forms for u1 (r i ), u2 (r i − r j ), andu3 (r i ,r j ,r i − r j ) as well as their optimization is done byadopting the CHAMP-code package [41] as it is. Then, theexpectation value A of an operator A is given by(14)
A VMC = |A| / | .

By putting A = n̂(r) = i δ(r − r i ) in Eq. (14), we obtainn(r) in VMC.
C. DMC

Starting with the variationally optimized wave function 
thus determined, we can further improve on the ground-state
wave function by considering the diffusion equation for (τ )in the imaginary time τ as∂(τ )(15)= (H − ER )(τ ),∂τwhere ER is the reference energy to be adjusted to E0 in thecourse of DMC simulations by removing the τ dependence
from the asymptotic form of (τ ) at τ → ∞. Note that theformal solution to Eq. (15) is written ase−(En −ER )τ |n n | ,(16)
(τ ) =−nn↓ (r). The spin density nσ (r)and n(r) is the sum of n↑ (r) and 
and consequently Nσ with N = σ Nσ should be determinedby the self-consistent solution of Eqs. (6)–(9), together withthe fixed boundary condition

iσ

i>j

i>j

i=1

φiσ (r) = 0,(6)where εiσ and φiσ are the energy level and the normalized
wave function for KS orbital i and spin σ , respectively, andvσKS (r) is the KS potential, determined byn(r )vσKS (r) = vext (r) + d r+ vσxc (r; [nσ ]), (7)|r − r |vσxc (r; [nσ ]) = δE xc [nσ ]/δnσ (r).

where the Jastrow function J (r 1 , . . . ,r N ) contains the terms
to describe electron-nucleus correlation u1 (r i ), two-electron
correlation u2 (r i − r j ), and three-body nucleus-two-electron
correlation u3 (r i ,r j ,r i − r j ) as
u1 (r i ) +u2 (r i − r j )J (r 1 , . . . ,r N ) =n(r)n(r )+ E xc [nσ ] + CNZ .
|r − r |(11)

where {n } is the normalized mutually orthogonal complete
set of eigenfunctions for H with the corresponding setof eigenenergies {En }. Then, as long as 0 | = 0, theasymptotic τ -independent wave function  is reduced tothe true ground-state wave function 0 , apart from thenormalization factor.
It is appropriate to add a comment on the condition of0 | = 0 here; by invoking the Anderson’s orthogonality
theorem [42], one may argue that 0 | vanishes in SSR,
but this is not correct for the reasons below; (i) both 0 and 
include the effect of the impurity atom, while the Anderson’s
theorem concerns with the relation between the wave
functions with and without the impurity. (ii) Due to thepresence of the Jastrow factor J ,  is not simply given by
the single Slater determinant 0 on which the Anderson’s
theorem is proved. (iii) DMC simulations are done for finite

155140-4

EMERGENCE OF A KONDO SINGLET STATE WITH KONDO . . .

N, while the Anderson’s theorem becomes valid only at
N → ∞.
In order to avoid the notorious fermion sign problem, we
employ the fixed-node approximation in DMC simulations.
This approximation may bring about undesirable errors in ,
but we try to minimize them by seeking for N -independent
results by exploiting the fact that the node positions dependon N in the fixed boundary condition, leading to the hope thatunphysical node-position-dependent effects will be removedby extracting the N -independent results. In performing actual
fixed-node DMC simulations at a fixed N, we adopt CHAMP
again to obtain the stably converged asymptotic wave function
(r 1 , . . . ,r N ). Then, n(r) in DMC is estimated by a secondorder approximation to the exact expectation value [40], whichamounts ton(r)DMC = 2 |n̂(r)| / | − n(r)VMC ,(17)

PHYSICAL REVIEW B 92, 155140 (2015)

the convergence of integrals:
δE =

 kF

δ (k) k 2
(2l + 1)
dk lσ
π 2
0
lσ
iσ ∈BS1Z
δn(r)δn(r )−d r δn(r) +
dr dr
r
2
|r − r |



n0
xc
xc
d r vσ (r; [nσ ])nσ (r) − vσ (r; [n0 /2])
−
2
σ


εiσ +

+ E xc [nσ ] − E xc [n0 /2] − EaZ ,(20)

where BS stands for the set of possible bound states. The
difference in the contribution of KS energies from continuum
states is treated by the consideration of change in the densityof states, δlσ (k)/π .

where n(r)VMC = n̂(r) VMC .

III. RESULTS FOR PROTON IMMERSION
A. LSDA for the bulk system

D. LSDA in the bulk system

Contrary to VMC and DMC, LSDA allows us to directly
treat the bulk (N → ∞) system, in which almost all statesin KS orbitals are continuum ones, for which we may write
εi = k2 /2 with momentum k and φiσ (r) = Rklσ (r)Ylm (r/r)
with the spherical harmonics Ylm (r/r) and the radial wave
function Rklσ (r) satisfying the following boundary condition
at r (≡|r|) → ∞:Rklσ (r) → cos[δlσ (k)]jl (kr) − sin[δlσ (k)]nl (kr),(18)

apart from a normalization factor, where jl (kr) and nl (kr) are
the spherical Bessel functions and δlσ (k) is the phase shift ofangular momentum l to be determined under the condition ofδlσ (∞) = 0, ensuring that an electron behaves as a free particle
at k (≡|k|) → ∞. There is a possibility of finding bound states
below the bottom of the conduction band (εiσ < 0) among
KS orbitals for which φiσ (r) → 0 at r → ∞. The Levinson
theorem [43] dictates that the total numberof the bound
states in KSorbitals NBS is given by lσ (2l + 1)δlσ (0)/π ,while Z = lσ (2l + 1)δlσ (kF )/π by the Friedel sum rule
[44].
E. Embedding energy

The embedding energy δE is defined as the difference ofthe ground-state energies between the atom-embedded EG andthe system of separated homogeneous EG and neutral atom.
Usually this concept is relevant only to the bulk EG. Thus,
denoting the ground-state energy of the isolated neutral atom
as EaZ , we can obtain δE throughδE = lim [E0 (N,Z) − E0 (N − Z,0)] − EaZ .
N→∞

(19)

Since δE is of the order O(1) in comparison with E0 of theorder O(N), due care must be exerted in order to accurately
evaluate δE at N → ∞. For this purpose, we rewrite Eq. (19)with using δn(r) [= n(r) − n0 ] and δlσ (k), the derivative ofδlσ (k) with respect to k, into the following form with ensuring

In line with the previous result [45], the ground statein the proton-embedded (Z = 1) bulk EG is perfectly spin
unpolarized, i.e., n↑ (r) = n↓ (r) at every r, at least for rs < 15
and our results on n(r) and δE in LSDA are in good agreementwith those in previous calculations [12–15,18,19,46]. In Fig. 4,
the obtained n(r) normalized by n(0) is plotted as a functionof 2kF r/π for rs = 1, 4, and 14, together with the s-wave
phase shift δs (k) which is spin independent and gives by far
the largest contribution among all δlσ (k)’s.
For rs < 1.97, we obtain δlσ (0) = 0 and thus NBS is zero,leading to the typical n(r) in CR with H+ screened by metallic
electrons in a short range. For rs  1.97, on the other hand,δs (0) = π and thus NBS = 2, seemingly implying the sudden
appearance of H − ion confined in EG at rs = 1.97 [12]. This
must be true, if nL (r) the localized-electron distribution isabout the same as n(r) for |r| smaller than the H − ion range,
as is the case for rs  12.5, in which δs (k) is distinctive andalmost a universal function of k/kF , as seen in Fig. 4(b). Notethat the deep dip in n(r) just outside the ion region, as seen for
rs = 14 in Fig. 4(a), is a typical electron profile describing therepulsive action of localized closed-shell electrons to exclude
metallic electrons from the ion region by the Pauli exclusion
principle.
For 1.97 < rs  12.5, however, the “localized” electrons
behave much differently; at rs = 4, for example, they extend
long up to |r| ≈ 82 and concomitantly nL (0) is much smaller
than n(0), indicating deep and massive penetration of itinerant
electrons into the proton site, but such penetration would
never be allowed due to the Pauli exclusion principle, if
the closed-shell H − ion were firmly constructed. Thus, byremembering that KS orbitals in DFT are not necessarily tiedwith real physical entities but just introduced for mathematicalconvenience to correctly reproduce n(r), we can assume that
this “localized” orbital does not represent a real localized state
but just describes the long-range change of n(r) in the formof an envelope density over ξK in SSR. Notice that δs (kF ) inthis density region lies between 0.7 × (π/2) and 1.3 × (π/2),
which means that δs (kF ) is close to π/2, a value expected in
the Kondo resonance state in IAM [22].

155140-5

YASUTAMI TAKADA, RYO MAEZONO, AND KANAKO YOSHIZAWA

10

-2

10

-3

(a) LSDA for the bulk system
rs =1

0.003

0.001
02
0.007

rs =14

0.006

n(0)=0.320
n0 =0.0000870
localized KS state

0.9

2
2k F r / π
rs =13
rs =14
rs =150.8

3

δs (k) / π

0.4
0.3

rs =1.8
rs =1.5

0.2

rs =1

8

10

r [aB]

1214

16(b) rs = 4

N =1700.003: DMC (N=170): VMC (N=170): LSDA(N=170): LSDA(N ∞)0.002

rs =10
rs =6

0.001
0

2

4

6

8

10

12

14

r [aB]

16

18

20

22

FIG. 5. (Color online) Examples of the calculated electron distribution n(r) obtained by DMC, VMC, and LSDA for (a) N = 60
and (b) 170 at rs = 4. The result in LSDA at N → ∞ is also shown.

implying that N = 60 is large enough to obtain the convergentn(r) near the proton site.

0.1(b) LSDA for the bulk system

0.0

0.2

0.4

6

0.004

4

rs =4
rs =3
rs =2.07

0.5

4

0.005

0.7
0.6

: DMC (N=60): VMC (N=60): LSDA(N=60): LSDA(N ∞)0.002n(0)=0.364
n 0 =0.00373“localized” KS state

1.0

N =60

0.004

rs =4

1

(a) rs = 4

0.005

n(0)=1.233
n 0 = 3 3 = 0.239
4πrs

10 -4
10 -500.006

n(r) [a B-3 ]

10

-1

0.007

n(r) [aB-3 ]n(r) /n (0)10 0PHYSICAL REVIEW B 92, 155140 (2015)

k / kF

0.6

0.8

C. Cusp theorem and the on-top density

1.0

According to the cusp theorem [47], n(r) near a nucleus of
atomic number Z behaves rigorously in the manner as

FIG. 4. (Color online) Electron distribution n(r) normalized byn(0) in (a) and spin-independent s-wave phase shift δs (k) in (b)obtained in LSDA for the proton-embedded bulk EG. In (a), thedotted curves show the contribution from the KS bound state which
is absent at rs = 1 must be a real physical state representing H− at
rs = 14, and appears only as mathematical convenience for describing
the long-range decrease in amplitude of the Friedel oscillation in SSR
at rs = 4. The label “localized” with a quotation mark indicates thissituation.
B. Comparison between LSDA and DMC

Before going into a more detailed discussion on SSR, let
us assess the accuracy of LSDA in comparison with DMC,
specifically at intermediate densities. At rs = 4, for example,
in Fig. 5, we see a good agreement between LSDA and DMC
for n(r) at any r, including the sphere boundary, at each N ,although the results in VMC do not match so well, assuring theimportance to perform DMC for taking the expectation valuesin accordance with Eq. (17).
As for N dependence or the size effect, we find that N =170 is not large enough to attain convergence in n(r) for |r| 
3.5 in both LSDA and DMC. For smaller |r|, however, noappreciable difference is seen between N = 60 and 170 in
DMC (and among N = 60, 170, and ∞ in LSDA) for rs  4,n(r) −−→ ncusp (r) ≡ n(0) exp(−2Zr),
r ≈0

(21)

where ncusp (r) exhibits strictly a linear change with r insemilog plots [the dashed lines in Figs. 6(a) and 6(b)]. By
exploiting this linear behavior, we can rather easily andaccurately determine the on-top density n(0) in DMC bylooking at the data in the region of 0.2  r  0.5, in spite of
the scattered nature of data points for n(r) at r  0.2 due to therapid increase in energy scale determined by the r −1 -Coulomb
potential.
For rs  11, n(0) thus obtained is about the same as theexact one for an isolated H− ion [39], as shown in Fig. 6(b)for rs = 11–14, indicating that the state is very close to the
H − ion confinement state. A more detailed observation on
the DMC data reveals a systematic change of the behavior atr  0.2 with the increase of rs ; for rs  12.5, the majorityof data points deviate upward, while opposite is the case forrs  12.5. The upward deviation indicates that the metallic
electrons rather easily penetrate into the core of the H− ion,
which is not allowed, once the closed-shell structure is solidly
constructed. On the other hand, the downward deviation isconsistent with the formation of the closed-shell ion, leading
to the conclusion that the transition to the H − ion confinement
state occurs at rs ≈ 12.5, the same rs as that in LSDA.

155140-6

EMERGENCE OF A KONDO SINGLET STATE WITH KONDO . . .

n(r) [a B-3 ]

0.40
0.30

: DMC (N= 60): DMC (N=170): LSDA(N= 60): LSDA(N=170): LSDA(N ∞)0.20
0.15 (a) r = 4
s

N = 60 & 170

0.10: Isolated H : cusp line for H : cusp line for DMC

0.05
0

0.2

n(r) [a B-3 ]

0.40
0.30

0.4

0.6

0.8

r [aB]

1.0

0.40(b) rs = 11-14 0.30
N = 60

0.40
0.30

rs=11

0.40
0.30

rs=12

0.40
0.30

0.20
0.15
0.10
0.05
0.04

rs=13

0.20
0.15

0.05
0.04

rs=14

0.10

0.05
0.04

: DMC (N=60): LSDA(N=60): Isolated H : cusp line for H -

0.05
0.04
0

0.2

0.4

0.6

r [aB]

0.8

0.05
0.04
1.0

FIG. 6. (Color online) Short-range part of n(r) in DMC andLSDA at (a) rs = 4 and (b) rs = 11−14. For comparison, the exact
result for an isolated negative hydrogen ion H− [39] is also plotted.

PHYSICAL REVIEW B 92, 155140 (2015)

N → ∞ no jump is seen in LSDA even at rs = 1.97 at which
the transition is known to occur through the abrupt change inδs (0).
In an extensive search for favorable N at which we
can easily obtain the convergent ground state, we findthat this jump becomes much enhanced at the magic
numbers (10, 60, 170, . . .) of N . Let us consider thereason for this fact; in LSDA, all KS levels are discrete at finite N and the stacking sequences of the occupied levels are (1s,2p,2s), (1s,2p,2s,3d,4f,5g,3p,3s), and(1s,2p,2s,3d,4f,3p,3s,5g,4d,6h,7i,5f,8k,4p,4s) for N =10, 60, and 170, respectively, reflecting the competition
between the −r −1 and r 2 potentials in Eq. (4) and the outermosts orbital, situating very near the Fermi level μ at each magic
number N , plays a key role in stabilizing the ground state. We
recognize that this outermost s orbital mimics the SSR state
at finite N , at least for R < ξK , in view of the fact that theactual SSR state in Kondo physics is situated at μ, extending
very long over the range ξK with an s-wave character. Then, thejump in n(0) is related to the transition from the empty s-orbital
state (corresponding to CR) to the s-orbital occupied one(corresponding to SSR) with changing N in 8 → 10, 58 → 60,
and 168 → 170 at rs = 3.85, 2.65, and 2.87, respectively, in
LSDA. In DMC, convergent results are also easily obtained forthe same series of N and the transition occurs, for example, at
rs ≈ 2.65 for N = 58 → 60, just as in LSDA. The next N inthis series is 340 associated with the 5s-orbital empty-occupied
transition, but DMC at this N is currently out of our reach.
At the second transition into the H− confinement state, achange in the stacking sequence occurs in LSDA. Thosetransition points, along with the second transition point
obtained in DMC by the data shown in Fig. 6(b), provide
the ground-state diagram in (rs , N −1 ) space in Fig. 3.
D. Embedding energy

In Fig. 7, we plot n(0) by changing rs and N and find thatthe results do not depend on N for rs  3, but they do forrs  3, i.e., in the CR-SSR transition region. The transition issignaled by a jump in n(0) in both LSDA and DMC, although
its magnitude decreases with increasing N and eventually atn(0) [aB-3]

1.2: DMC (N=58 or 60): DMC (N=170)

1.0
0.8
0.6

Transition at N ∞ at rs =1.97 in LSDA
Transition (N=58 60) at rs =2.65 in LSDA
Transition (N=168 170) at rs =2.87 in LSDA
Isolated H-

2nd Transition
at rs ~
~12.5

0.4

In calculating δE by use of Eq. (19) at some finite N ,we need E0 (N − 1,0) in addition to E0 (N,1), but for aneven number of N , the ground state of the (N − 1)-electron
system is necessarily a spin-polarized one. Nevertheless, forsmoothly connecting to the spin-unpolarized ground state
at N → ∞, it is better to calculate in the spin-unpolarizedsituation. In LSDA, by a spin-symmetrized sum of KS orbitals
by fractionally occupying the orbitals at the Fermi level, we
can obtain this needed spin-unpolarized ground-state energy
E0 (N − 1,0) with which we can calculate δE for finite-N
systems.
In DMC, however, we cannot adopt this procedure. Thus,
we first calculate the spin-unpolarized ground-state energy
E0 (N − 2,0) in DMC, starting with the Slater determinant
0 (r 1 , . . . ,r N−2 ) without the outermost s orbital. Then, we
estimate E0 (N − 1,0) by
E0 (N − 1,0) =1 2 3 4 5 6 7r 8 9 10 11 12 13 14s

FIG. 7. (Color online) On-top density n(0) plotted as a functionof rs for N = 60 (and 170 only for rs = 4) in DMC and for N = 60,
170, and ∞ in LSDA, together with the exact results for an isolated
H− ion (N = 2). The transition points are indicated by arrows. For
the data in DMC, errors are within the size of square symbols.

N −1
αN−1
,
E0 (N − 2,0) 1 +N −2(N − 1)1/3(22)

with the coefficient αN−1 determined through E0 (N − 1,0)and E0 (N − 2,0) in LSDA. We have deduced this approximation scheme by considering that the leading term in
E0 (N − 1,0) is, in general, in proportion to (N − 1) due to155140-7

YASUTAMI TAKADA, RYO MAEZONO, AND KANAKO YOSHIZAWA
0.0042: LSDA (N ∞): DMC (N=58 or 60): DMC (N=170)

0.6

0.0038

0.4 Transition at N ∞ at rs =1.97 in LSDA

0.0036

Transition (N=58 60) at rs =2.65 in LSDA
Transition (N=168 170) at rs =2.87 in LSDA
2nd Transition
Isolated H at rs ~~12.5

0.00320.00300.0028s

0.00262

FIG. 8. (Color online) Embedding energy δE calculated in thesame situations as those for n(0) in Fig. 7. Errors in DMC are within
the size of square symbols.

extensiveness of the total energy as well as the subleading
term in proportion to (N − 1)2/3 due to the surface-energy
contribution. In CR in which the outermost s orbital is empty
even for the system with Z = 1 and N electrons, instead of
E0 (N − 2,0), we calculate E0 (N,0) in DMC, with which weestimate E0 (N − 1,0) by a similar strategy.
In Fig. 8, δE is given as a function of rs , exhibiting the N
dependence similar to that in n(0) in Fig. 7, including jumps
at CR-SSR transitions. Our results for δE agree reasonably
well with previous ones [12–19], although due attention to
the size dependence was not paid previously. These results of
δE demonstrate that hydrogen is most stably embedded in EG
in the form of SSR. Thus, the concept of SSR is deemed toplay a key role in hydrogen storage in metals and hydrogen is
expected to reside at a site with rs ≈ 4 in an inhomogeneous
metal.

4

6

8

10

12

14

r [aB]

16

18

20

22

FIG. 9. (Color online) Long-range part of n(r) in both DMC andLSDA obtained by the elimination of the sphere-boundary effect.

An example of n(r) obtained through the above procedure
at rs = 4 is given in Fig. 9, in which we find good convergencefor r  7 in changing N by the comparison between N = 60
and 170, illustrating that this procedure doubles the sizeconvergent range of r as compared with that in Fig. 5.
Admittedly, discrepancy is seen in the oscillation amplitude
between DMC and LSDA, but overall good agreement andsize convergence are obtained in the Friedel-oscillation phase,indicating that accurate enough information is now available
on the node positions in δn(r) up to r ≈ 13. Incidentally, theoutermost s orbital contributes much to the oscillation behavior
in Fig. 9, assuring its importance in the SSR-density region.
In order to unambiguously confirm the emergence of SSR,
let us examine the Friedel oscillation in n(r) in the light of its
general behavior, known as [23–26]
1 
cos 2kF r − 3π/2 + 2δs(0) F (r/ξK )4π 2 r 3− cos(2kF r − 3π/2) .
(26)n(r) −−−→n0 +E. Anomalous Friedel oscillationrkF−1

Basically, the Friedel oscillation due to the presence of a
proton at the origin is a concept defined in the bulk system,
but if we try to discuss it with use of the data in DMC,
we need to eliminate the sphere-boundary effect from n(r)obtained in finite-N systems. For this purpose, we adopt
the following procedure [17]; we first calculate the charge
distribution nN,Z (r) in the finite-N system corresponding to
E0 (N,Z) and then we estimate n(r) through the cancellationof the sphere-boundary effect by subtracting nN−1,0 (r) from it
as(23)

In LSDA, we can employ Eq. (23) by using nN−1,0 (r) which
is obtained in the spin-unpolarized situation simultaneouslywith E0 (N − 1,0), but in DMC, only nN−2,0 (r) associatedwith E0 (N − 2,0) is available. Thus, we estimate nN−1,0 (r)from this nN−2,0 (r) asβ(r) DMCnDMC
nN−2,0 (λr),
N−1,0 (r) = λ: DMC (N= 60): DMC (N=170): LSDA(N= 60): LSDA(N=170): LSDA (N ∞): n0 = 0.00373

0.00340.0-0.11 2 3 4 5 6 7 8 9 10 11 12 13 14rn(r) ≈ nN,1 (r) − nN−1,0 (r) + n0 .

rs = 4 kF = 0.480

0.0040n(r) [a B-3 ]δE [ha]0.8

0.2

PHYSICAL REVIEW B 92, 155140 (2015)(24)

where λ is a parameter to correct the difference in thesphere radius, defined as λ = [(N − 2)/(N − 1)]1/3 , and ther-dependent exponent β(r) is determined with use of the data
for the electron distribution in LSDA as
LSDAβ(r) = ln nLSDA
ln λ.
(25)
N−1,0 (r) nN−2,0 (λr)

Here, only the s-wave contribution, which indeed dominates
others in the present case, is considered and δs(0) is the s-wave
phase shift at the Fermi level produced by the potential
scattering without the Kondo-resonance effect. In Eq. (26),F (r/ξK ) ≡ 1 in CR, but in SSR it gradually decreases from
1 for r  ξK to −1 for r  ξK due to physics of asymptotic
freedom [24]. Actually in SSR, F (x) is known to be a universal
scaling function, as explicitly given in Fig. 10(a).
By appropriately choosing the branch of tan−1 x, we canrewrite Eq. (26) inton(r) −−−→ n0 −rkF−1A(r)cos[2kF r + δ(r)],4π 2 r 3

with the amplitude A(r) and the phase δ(r), given byA(r) = 1 − 2F (r/ξK ) cos 2δs(0) + F (r/ξK )2 ,δ(r) = tan−1

1 − F (r/ξK ) cos 2δs(0)F (r/ξK ) sin 2δs(0)

.

(27)

(28)
(29)

If δs(0) is in the range (0,π/2), δ(r) increases gradually fromδs(0) for small r to δs(0) + π/2 for large r in accordance

155140-8

EMERGENCE OF A KONDO SINGLET STATE WITH KONDO . . .

with the change of F (r/ξK ). On the other hand, if δs(0) isin the range (π/2,π ), δ(r) decreases gradually from δs(0)for small r to δs(0) − π/2 for large r. With taking care ofsuch a gradual change in δ(r), we can determine ri the ith
zero of δn(r) [≡ n(r) − n0 ] by 2kF ri + δi = iπ + π/2 with
δi ≡ δ(ri ) and i = 1,2,3, . . . . Then, the half period of theFriedel oscillation i in SSR is given by i = ri+1 − ri =[1 − (δi+1 − δi )/π ](0) with (0) (≡π/2kF ) the half period in
CR. Because δi+1 is slightly different from δi , i is modified
from (0) by the amount of (δi+1 − δi )/π . This anomalous
Friedel oscillation or the oscillation-period modulation effectis an important consequence of the presence of SSR, but thismodulation effect has not been well recognized even in the
Kondo-physics community due to the fact that the effect istotally absent for δs(0) being equal to a multiple of π/2 [see thedenominator in Eq. (29)], which happened to be assumed in
the previous model calculations [24,25].

In Fig. 10(b), we plot δn(r)r 3 /Z as a function of 2kF r/π at
rs = 4 to check whether the modulation effect exists or not in
our first-principles calculations. In DMC, the first four zerosof δn(r) in units of (0) , 2kF ri /π , are given by 0.994 ± 0.005,
1.970 ± 0.008, 2.949 ± 0.013, and 3.927 ± 0.019 with theerrors estimated by the distribution of data points aroundδn(r) = 0 at each ri . Then, (δi+1 − δi )/π the modulation
factors for i = 1, 2, and 3 are, respectively, obtained as
0.024 ± 0.005, 0.021 ± 0.009, and 0.022 ± 0.016 in which
the errors are estimated by the inclusion of covariance between
ri and ri+1 . Those results, at least definitely those for i = 1
and 2, assure the existence of the shortening of the Friedeloscillation period, confirming the emergence of SSR.
Once the data for {ri } are known, we can independently
calculate the modulation factors through Eq. (29) as a functionof ξK and δs(0) with using F (x) in Fig. 10(a). The results areshown in Fig. 10(c), from which we find that our DMC data forthe modulation factors (actually the shortening factors in thiscase) for i = 1–3 agree very well with those obtained at ξK ≈35/kF = 73 and δs(0) = 0.86 × (π/2) [see the horizontal andvertical double-dotted-dashed lines in Fig. 10(c)], indicatingthat by quantitatively analyzing first few modulation factors,
we can determine both ξK and δs(0) uniquely, even if only thedata for the system size much shorter than ξK are available.
This determination is made possible due to the fact that
F (r/ξK ) changes most rapidly at r ≈ ξK /10. By using ξK thus
determined and the Fermi velocity vF , the Kondo temperature
TK is estimated by TK = vF /ξK ≈ 0.0066 Hartree = 2100 K
with about 10% errors, revealing that this is an astonishingly
high-TK system. With such a high TK , the system will not
exhibit the prominent Kondo effects such as the ln T anomaly
in the experiment for T around the room temperature or below;
it just behaves as a usual Fermi liquid [48].
The shortening effect in the Friedel-oscillation period is
also found in LSDA, as seen in Fig 11, in which δn(r)r 3 isplotted for both Z = 1 and −1 in the bulk system at rs = 4. For
large enough r outside the SSR binding radius ξK , the Friedeloscillation phase for Z = 1 coincides with that for Z = −1;
the result with Z = −1 is plotted to represent the behavior for
the fictitious very tightly bound H − ion confinement state inwhich the anomalous Friedel oscillation is absent. In LSDA,

0.08

δn(r)r 3, n L(r) r3

FIG. 10. (Color online) (a) Universal scaling function characterizing SSR [24,25]. (b) Friedel oscillation in SSR in the protonembedded EG at rs = 4, as seen by the plot of δn(r)r 3 /Z in both DMC
and LSDA to detect the half period of the oscillation in comparison
with that in CR obtained in the system with Z = 0.01, a fictitious tiny
charge to ensure that the ground state can be accurately obtained in
LSDA as a CR state even at rs = 4. (c) First three modulation factors
in the half period of the anomalous Friedel oscillation, given throughEq. (29) with use of F (x) in (a), plotted as a function of the s-wave
potential scattering phase shift δs(0) for various values of kF ξK .

PHYSICAL REVIEW B 92, 155140 (2015)nL (r)r 3 :“localized” Kohn-Sham state@Z=1

0.04
0

-0.04
-0.08
-0.12
0

Z=1

LSDA (N5∞) @ rs =4

δn(r)=n(r)−n0 k F =0.480

Z= -1
10

2kF r/π

15

20

25

FIG. 11. (Color online) Plot of δn(r)r 3 in LSDA in the bulk
system with both Z = 1 and −1 at rs = 4. The contribution fromthe “localized” KS state obtained at Z = 1, δnL (r)r 3 , is also given
and seen as playing a role of an envelope density.

155140-9

YASUTAMI TAKADA, RYO MAEZONO, AND KANAKO YOSHIZAWA

the modulation of the Friedel oscillation is brought about by thecontribution from the “localized” density nL (r). Incidentally,
the oscillation-period shortening factors in LSDA are muchsmaller than those in DMC; (δ2 − δ1 )/π = 0.0133 and allothers are less than 0.004, indicating δs(0) ≈ 0.014 × (π/2)
in LSDA. [We note that δs (kF ) in Fig. 4(b) is equal toδs(0) + π/2 ≈ 0.507π .] We can easily understand the reason forthis large difference in δs(0) between DMC and LSDA; because
this phase shift is directly connected with the wave function, its
accurate value will not be obtained by LSDA in which the wave
function in the fictitious noninteracting system is qualitatively
different from the true correlated SSR wave function. Due toδs(0) ≈ 0 in LSDA, we cannot employ the diagram in Fig. 10(c)to determine ξK very accurately. Therefore, we estimate ξKfrom the extent of nL (r), which is 82, giving TK to be about
1900 K. In relation to nL (r), εBS the binding energy for this“localized” KS state in LSDA is given as 0.0115 Hartree =3600 K, which is about twice as large as TK . Thus, although
εBS has no direct physical meaning, this quantity seems to bea good measure for the magnitude of TK . In our calculationsin LSDA, the values for εBS are 97, 2200, 3500, 2200, and 720
K for rs = 2.07, 3, 6, 8, and 10, respectively, indicating thatin the majority of the SSR region, namely, for 3  rs  8, we
may expect TK to be well beyond 1000 K.
IV. DISCUSSION

Five comments on this work are in order:(i) In SSR, we find a good semiquantitative agreement
between LSDA and DMC, but this is by no means fortuitous
because this can be understood by the long-range nature ofξK which makes the density variation associated with the SSR
state slow, validating the use of LSDA for the calculationof n(r).
(ii) From our present perspective, we may regard ourprevious study in LSDA on the spin-polarized ground states
for second-period atoms in periodic table [45] as a successful
extension to multichannel Kondo systems [1], in which the
Hund’s-rule coupling plays a crucial role in producing thespin-polarized ground states.
(iii) In a short term, the immediate next target of research isa hydrogen molecule H2 immersed in EG [12] as a function of
the interatomic distance to pursue a new concept in chemicalbonding [49,50] in metals in the light of SSR. In fact, Bonev
and Ashcroft [51] have already found an interesting bistability
between paired and unpaired states for rs > 3.2. In a longer
term, we can expect fruitful research on new aspects in the
Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction [1] anddilute magnetic semiconductors [52].
(iv) If the number of immersed protons is increased up to a
macroscopic level and those protons are arranged into a latticein a metal with the SSR-density region, we might be able toobtain a Kondo lattice with high TK . Furthermore, if the lattice
constant might be so arranged as to tune the RKKY interaction
to be about the same as TK for realizing the quantum-

PHYSICAL REVIEW B 92, 155140 (2015)

critical situation in the Doniach phase diagram [53], then as
indicated in Fig. 1 in the optimum situation, we speculate thatsuperconductivity occurs at a rather high temperature in the
spin-fluctuation mechanism in a metallic hydrogen alloy. This
speculation may motivate people engaged in hydrogen-based
physics, chemistry, and technology to synthesize hydrogenalloys exhibiting exotic superconductivity at ambient pressure.
In this regard, the transition-metal-hydride system such as
TiH2 and ZrH2 [54–56], usually used for secondary batteries,
might be a promising candidate, although it seems thatthe metallic electron densities in the metal hydrides so far
synthesized are too high for our purpose.
(v) The basic reason for about 100 times increase of TK in
the hydrogen system compared with those in the f -electron
heavy-fermion systems is the overall increase of energy scales
as seen by the large difference in the magnitudes of 1s and 4f
energy levels. For the same reason of increased energy scales,
superconductivity with Tc over 100 K has been discussed insolid hydrogen at pressures of about 500 GPa in the conventional phonon mechanism [57–60]. Very recently, stimulatedby the experimental result [61], similar discussions are madeon H2 S [62,63] as well as H3 S [64] at pressures of about
200 GPa. Note that the metallic densities in those systems arefound to be in the CR region, namely, rs < 2 (or typicallyrs ≈ 1.4). Therefore, our present proposal has nothing to dowith those conventional theories for superconductivity. We
also emphasize that high metallic densities realized at suchvery high pressures in solid H2 or sulfur hydrides are notneeded in our proposal for finding exotic superconductivity.
V. SUMMARY

In conclusion, by employing both DMC and LSDA, wehave revealed the emergence of SSR in the proton-embedded
electron gas by the confirmation of an anomalous Friedel
oscillation characteristic to the Kondo-type spin-singlet statewith quantitatively determining TK and emphasized its stability
in embedding hydrogen into the electron gas. Our work
necessitates to modify the paradigm of CR in metallic
screening to a point charge of Z in the sense that, dependingon the metallic electron density and Z, SSR takes the placeof CR. This work also provides a first concrete example toshow how the KS orbitals behave to represent the exact n(r) instrongly correlated electron systems such as those describing
Kondo physics, even though they are not always physically
relevant.
ACKNOWLEDGMENTS

Y.T. thanks A. Savin and E. K. U. Gross for valuable
discussions. This work is partially supported by Innovative
Area “Materials Design through Computics: Complex Correlation and Non-Equilibrium Dynamics” (Grant No. 22104011)
from Ministry of Education, Culture, Sports, Science, andTechnology, Japan.

[1] A. C. Hewson, in The Kondo Problem to Heavy Fermions,Cambridge Studies in Magnetism, edited by D. Edwards and155140-10

D. Melville (Cambridge University Press, Cambridge, England,1993).

EMERGENCE OF A KONDO SINGLET STATE WITH KONDO . . .
[2] G. R. Stewart, Rev. Mod. Phys. 73, 797 (2001); 78, 743 (2006).
[3] H. v. Löhneysen, A. Rosch, M. Vojta, and P. Wölfle, Rev. Mod.
Phys. 79, 1015 (2007).
[4] P. Gegenwart, Q. Si, and F. Steglich, Nat. Phys. 4, 186 (2008).
[5] M. B. Maple, E. D. Bauer, V. S. Zapf, and J. Wosnitza,
in Superconductivity, edited by K. H. Bennemann and J. B.
Ketterson (Springer, Berlin, 2008), Chap. 13.
[6] O. Stockert, S. Kirchner, F. Steglich, and Q. Si, J. Phys. Soc.
Jpn. 81, 011001 (2012).
[7] F. Steglich, O. Stockert, S. Wirth, C. Geibel, H. Q. Yuan, S.
Kirchner, and Q. Si, J. Phys.: Conf. Ser. 449, 012028 (2013).
[8] E. D. Bauer, J. D. Thompson, J. L. Sarrao, L. A. Morales, F.
Wastin, J. Rebizant, J. C. Griveau, P. Javorsky, P. Boulet, E.
Colineau, G. H. Lander, and G. R. Stewart, Phys. Rev. Lett. 93,
147005 (2004).
[9] N. J. Curro, T. Caldwell, E. D. Bauer, L. A. Morales, M. J.
Graf, Y. Bang, A. V. Balatsky, J. D. Thompson, and J. L. Sarrao,
Nature (London) 434, 622 (2005).
[10] P. Debye and E. Hückel, Physikalische Zeitschrift 24, 185(1923).
[11] A. Pundt and R. Kirchheim, Annu. Rev. Mater. Res. 36, 555(2006) and the references therein.
[12] J. K. Nørskov, Phys. Rev. B 20, 446 (1979).
[13] M. J. Stott and E. Zaremba, Phys. Rev. B 22, 1564 (1980).
[14] M. J. Puska, R. M. Nieminen, and M. Manninen, Phys. Rev. B
24, 3037 (1981).
[15] J.-H. Song, Ph.D. thesis, Oregon State University, 2004,
http://hdl.handle.net/1957/29170[16] G. Sugiyama, L. Terray, and B. J. Alder, J. Stat. Phys. 52, 1221
(1988).
[17] A. I. Duff and J. F. Annett, Phys. Rev. B 76, 115113 (2007).
[18] C. O. Almbladh, U. von Barth, Z. D. Popovic, and M. J. Stott,Phys. Rev. B 14, 2250 (1976).
[19] E. Zaremba, J. H. Rose, L. M. Sander, and H. B. Shore, J. Phys.
F: Metal Phys. 7, 1763 (1977).
[20] P. Jena, A. K. Gupta, and K. S. Singwi, Phys. Rev. B 18, 2723
(1978).
[21] K. Yosida, Phys. Rev. 147, 223 (1966).
[22] P. W. Anderson, Phys. Rev. 124, 41 (1961).
[23] F. Mezei and G. Grüner, Phys. Rev. Lett. 29, 1465 (1972).
[24] I. Affleck, L. Borda, and H. Saleur, Phys. Rev. B 77, 180404(R)(2008).
[25] G. Bergmann, Phys. Rev. B 78, 195124 (2008).
[26] Y. Tao and G. Bergmann, Eur. Phys. J. B 85, 42 (2012).
[27] C. M. Varma and Y. Yafet, Phys. Rev. B 13, 2950 (1976).
[28] A. M. Oleś and K. A. Chao, Phys. Status Solidi B 98, 271 (1980).
[29] F. Sottile and P. Ballone, Phys. Rev. B 64, 045105 (2001).
[30] Y. Takada, Phys. Rev. A 28, 2417 (1983).
[31] Y. Takada, Phys. Rev. B 35, 6923 (1987).
[32] Y. Takada, Phys. Rev. B 43, 5979 (1991).

PHYSICAL REVIEW B 92, 155140 (2015)[33] Y. Takada, Phys. Rev. Lett. 87, 226402 (2001).
[34] H. Maebashi and Y. Takada, Phys. Rev. B 84, 245134 (2011).
[35] R. M. Dreizler and E. K. U. Gross, Density Functional Theory(Springer, Berlin, 1990), Chap. 4.
[36] G. Stefanucci and S. Kurth, Phys. Rev. Lett. 107, 216401 (2011).
[37] O. Gunnarsson and B. I. Lundqvist, Phys. Rev. B 13, 4274
(1976).
[38] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev.
Mod. Phys. 73, 33 (2001).
[39] C. J. Umrigar and X. Gonze, Phys. Rev. A 50, 3827 (1994).
[40] R. N. Barnett and K. B. Whaley, Phys. Rev. A 47, 4082 (1993).
[41] Cornell-Holland Ab initio Materials Package (CHAMP) written by C. J. Umrigar, C. Filippi, and J. Toulouse,
http://www.physics.cornell.edu/cyrus/champ.html[42] P. W. Anderson, Phys. Rev. Lett. 18, 1049 (1967).
[43] Z.-Q. Ma, J. Phys. A: Math. Gen. 39, R625 (2006).
[44] J. Friedel, Philos. Mag. 43, 153 (1952).
[45] V. U. Nazarov, C. S. Kim, and Y. Takada, Phys. Rev. B 72,
233205 (2005).
[46] P. Jena and K. S. Singwi, Phys. Rev. B 17, 3518 (1978).
[47] A. E. Carlsson and N. W. Ashcroft, Phys. Rev. B 25, 3474
(1982).
[48] P. Nozières, J. Low Temp. Phys. 17, 31 (1974).
[49] Y. Takada and T. Cui, J. Phys. Soc. Jpn. 72, 2671 (2003).
[50] M. Shimomoto and Y. Takada, J. Phys. Soc. Jpn. 78, 034706(2009).
[51] S. A. Bonev and N. W. Ashcroft, Phys. Rev. B 64, 224112(2001).
[52] T. Dietl, Nat. Mater. 9, 965 (2010).
[53] S. Doniach, Phys. B (Amsterdam) 91, 231 (1977).
[54] Q. Xu and A. Van der Ven, Phys. Rev. B 76, 064207 (2007).
[55] R. Quijano, R. de Coss, and D. J. Singh, Phys. Rev. B 80, 184103(2009).
[56] D. Chattaraj, S. C. Parida, S. Dash, and C. Majumder, Int. J.
Hydrogen Energy 39, 9681 (2014).
[57] N. W. Ashcroft, Phys. Rev. Lett. 21, 1748 (1968).
[58] N. W. Ashcroft, Phys. Rev. Lett. 92, 187002 (2004).
[59] P. Cudazzo, G. Profeta, A. Sanna, A. Floris, A. Continenza, S.
Massidda, and E. K. U. Gross, Phys. Rev. B 81, 134506 (2010).
[60] J. M. McMahon and D. M. Ceperley, Phys. Rev. B 84, 144515(2011).
[61] A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov, andS. I. Shyein, Nature (London) 525, 73 (2015).
[62] Y. Li, J. Hao, H. Liu, Y. Li, and Y. Ma, J. Chem. Phys. 140,
174712 (2014).
[63] I. Errea, M. Calandra, C. J. Pickard, J. Nelson, R. J. Needs, Y.
Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, Phys. Rev. Lett. 114,
157004 (2015).
[64] D. Duan., Y. Liu, F. Tian, D. Li, X. Huang, Z. Zhao, H. Yu, B.
Liu, W. Tian, and T. Cui, Sci. Rep. 4, 6968 (2014).

155140-11