{'fraction_non_alphanumeric': 0.048400235326062364, 'fraction_numerical': 0.019550165180793774, 'mean_word_length': 4.0387641089955535, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Emergence of a Kondo singlet state with the Kondo temperature well beyond 1,000K in the proton-embedded electron gas: Possible route to high-Tc superconductivity

Hydrogen in metals has attracted much attention for a long time from both basic scientific and technological points 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 is well 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 distribution n(r) by the combination 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-like spin-singlet resonance, the emergence of which is confirmed by the presence of an anomalous Friedel oscillation characteristic to the Kondo singlet state with the Kondo temperature T_K well beyond 1,000K. 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 long-sought novel high-T_K system, opening an unexpected route to room-temperature superconductivity in a Kondo lattice made of protons.


I. INTRODUCTION
Physics in the heavy fermion superconductors has been understood by the concept of quantum criticality in a system of 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 T c and the Kondo temperature T K 4-7 . More specifically, T c is of the order of 0.1T K , as shown in Fig. 1 plotted based on the information available in the literature, from which we can conceive an idea that high-T c superconductivity will be obtained if we can discover a Kondo system with very high T K . In fact, the recentlydiscovered plutonium compounds such as PuCoGa 5 with T c = 18.5K and T K ≈ 260K 8,9 may be regarded as the successful realization of this idea. Thus we should make further pursuit of this idea by searching for a new class of Kondo systems with T K higher than 1,000K. Theoretically, this search can be done by the first-principles quantitative determination of T K for the composite system of an impurity atom embedded in a metal.
As first suggested by Debye and Hückel 10 , an atomic nucleus charge +Ze in a metal is screened by accumulation of metallic electrons which is regarded as a charge resonance and well described by the linear response theory. This concept of 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 is an odd integer, such as Z = 1 (case of a proton), in which a spin-polarized bound state might appear at the impurity atom. With the above basic scientific issues in mind, we have concerned 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 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][13][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 by quantum Monte Carlo methods [15][16][17] . Nonlinear metallic screening is another extensively examined topic in this system [18][19][20] , but no serious attention has been paid so far to the spin resonance effect.
The homogeneous EG with the average density n 0 is specified by a single parameter r s , defined by r s = (3/4πn 0 ) 1/3 in units of the Bohr radius a B . (We use atomic units hereafter.) Its characteristic energy is the Fermi energy ε F , given by k 2 F /2 (= 1.84r −2 s ) with k F (= 1.92r −1 s ) the Fermi momentum ( Fig. 2(a)). Hydrogen, on the other hand, has two typical energies ( Fig. 2(b)), the 1s level ε 1s (= 0.5) and the electron affinity ε A (= 0.0278). Ratio of ε F to ε 1s or ε A determines the character of the ground state in the proton-embedded EG; for high n 0 corresponding to ε F ≫ ε 1s (or r s ≪ 2), the 1s level is buried in the continuum of EG and thus no electrons are bound to a bare proton H + , leading to a charge resonance (CR) state in which H + is dielectrically screened by accumulation of itinerant electrons near the Fermi level µ (Fig. 2(c)). For low n 0 with ε F ≪ ε A (or r s ≫ 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 large that the localized electrons in H − begin to spill out into EG. Then a crucial question is whether this state at ε F ≈ ε A is the same as that in Fig. 2(c) or not. . Three possible ground states are schematically illustrated in (c)-(e) corresponding, respectively, to a bare proton H + screened by metallic electrons with the screening length ≈ k −1 F (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 ≫ k −1 F . These temporary states are superposed with interchanging the roles of spins to make a resonance state ΨSSR.
Intuitively, for ε A ε F ε 1s , we can imagine a spinpolarized state made of a single electron with either up or down spin captured by H + , but in view of the concept of spin 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-like spinsinglet resonance (SSR) state Ψ SSR (Fig. 2(e)). Because there is no clear distinction between conduction and localized electrons, this SSR state is composed of only itinerant electrons near µ without a local spin moment, similar to CR, but an important difference exists in the screening length; for CR, it is the Thomas-Fermi length λ TF ≈ k −1 F , but for SSR, the Kondo-screening length ξ K is much longer than k −1 F , leading to an anomalous Friedel oscillation [23][24][25][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 T K from first principles, but this confirmation is not an easy task due to the existence of various difficulties, as we shall explain below in some detail.
The Hamiltonian for IAM, H A , is written as 22 in second quantization with use of the annihilation operator c kσ for a conduction electron with wave vector k, spin σ and one-body band energy ε k , while c dσ is an operator to destroy a localized d electron with spin σ at the energy level E d . The d electrons not only interact to each other at the localized site with the strength U but also hybridize with the conduction electrons with the strength V dk . Based on H A , Ψ SSR is given as 27,28 where Φ 0 is the Slater determinant made of conductionelectron orbitals and the parameters, a 0 and a k , are to be determined variationally. As Eq. (2) clearly shows, Ψ SSR is not descibed 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 indispensable for the construction of this Ψ SSR . In our work, we are not allowed to employ this simple model H A . Instead, we have to start with the firstprinciples Hamiltonian H, described in atomic units as in first quantization. Here we have considered a neutral atom of atomic number Z at the origin of coordinates immersed into the jellium sphere of radius R and average density n 0 . The number of electrons contained in the jellium sphere is 4πR 3 n 0 /3 = (R/r s ) 3 , so that the total electron number N is equal to Z + (R/r s ) 3 , satisfying global neutrality, from which we obtain R = (N−Z) 1/3 r s . In Eq. (3), the external potential working on an electron v ext (r) is composed of the potential from the nucleus and that from the positive background, written as with θ(x) the Heaviside function and C Z N represents the Coulomb self-energy stemming from both nucleusbackground and intra-background interactions, given by In solving Eq. (3), we impose the fixed boundary condition to make the wave function vanish at |r i | = R. Irrespective of whether we include the constant term C Z N in H or not, there is no problem of divergence in this finite-N system, but C Z N is needed in order to achieve the mutual global cancellation in energies between the electronbackground attractive potentials and the repulsive potentials of both electron-electron and intra-background interactions for the bulk (N → ∞) system 29 .
By comparing H in Eq. (3) with H A , we readily see their differences; first, there is no predetermined localized state in H and thus a (virtual) localized state, if any, must be determined in the first place in constructing Ψ SSR with related parameters such as E d , U , and V dk , if necessary. Second and more importantly, the long-range Coulomb interaction, which is completely neglected in H A , works among metallic electrons in H, making solution of the problem quite difficult. In fact, even in the homogeneous EG without the embedded atom, the problem is sufficiently complicated in both variational [30][31][32] and Green's-function 33,34 approaches. Incidentally this long-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 another important aim of this paper, we can never neglect this long-range Coulomb interaction, making all theoretical and computational techniques developed so far for H A useless to H.
Usually, the first-principles Hamiltonian is solved by either diffusion Monte Caro (DMC) simulations or DFTbased 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 can easily treat the bulk system, but the ground-state physical quantities are calculated in terms of a single Slater determinant made 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 we always have to resort to some approximation to the exchange-correlation 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) and the corresponding ground-state energy E 0 by projecting the real interacting many-body system to a fictitious auxiliary non-interacting system in which n(r) can be calculated with use of a single Slater determinant made of KS orbitals, even if we know 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 non-degenerate 35 , as is the case for SSR. Of course, information obtained only through n(r) and E 0 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 and SSR states and determining T K 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 spinresolved electron distribution n σ (r), so that we can check a possibility of the spin-polarized (i.e., non-spinsinglet) ground state by detecting the difference between n ↑ (r) and n ↓ (r). (ii) We assess the results in LSDA at finite-N systems in comparison with those in fixed-node DMC 38 with taking N up to 170, a much larger size than those in previous variational Monte Carlo (VMC) calculations 16,17 . (iii) We obtain n(r) in the 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 can assume that they will be free from any restrictions incurred from the fixed-node approximation, the only approximation adopted in DMC simulations, because the fixed-node positions in DMC are prescribed by R (and consequently by N for given r s due to R = (N −Z) 1/3 r s ) 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, a very important inherent property of the anomalous Friedel oscillation, to which we call serious attention for the first time in quantitatively determining ξ K and consequently T K from first principles.
In accordance with those strategies, we have investigated n(r) in finite-N systems in both LSDA and DMC to find spin-unpolarized ground states with a strong size effect and a series of magic numbers (10, 60, 170, · · · ) of N at which convergent results are easily obtained. These features can be explained in terms of the emergence of SSR with its long-range nature of ξ K . Its emergence is also signaled in LSDA at N → ∞ for r s 2 by the appearance of a strange shallow bound KS orbital below the conduction band bottom with an unusually long binding radius. In spite of the large size effect, for each N , n(r) in LSDA agrees very well with that in DMC.

FIG. 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 sharp transition. 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.
By summarizing the results thus calculated, we have obtained a ground-state diagram in (r s ,N −1 ) space, shown in 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-range change of n(r) over ξ K in SSR in the form of an envelope density. We also find that hydrogen is most stably embedded in EG in the SSR region, especially, optimally firmly at r s ≈ 4 with the Kondo temperature T K ≈2,100K, indicating our success in discovering a long-sought high-T K system. Thus we may expect that superconductivity occurs at a temperature as high as about 0.1T K in a metallic hydrogen alloy at ambient pressure in which a macroscopic number of protons are regularly embedded in a metal in this density region to form a Kondo lattice.
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 their implications and future directions. Finally in Sec. V, we give a summary of this paper.

A. LSDA in the finite-N system
Let us consider the neutral system of a single nucleus with atomic number Z embedded in the N -electron jellium sphere of radius R. Its Hamiltonian H is given in Eq. (3). In LSDA to DFT, the KS equation is written as where ε iσ and φ iσ are the energy level and the normalized wave function for KS orbital i and spin σ, respectively, where v xc With use of the lowest-N σ KS orbitals, n σ (r) is given by and n(r) is the sum of n ↑ (r) and n ↓ (r). The spin density n σ (r) and consequently N σ with N = σ N σ should be determined by the self-consistent solution of Eqs. (6)-(9), together with the fixed boundary condition at |r| = R = (N − Z) 1/3 r s . By using those converged quantities, we can calculate E 0 (N, Z) the ground-state energy including the constant term C Z N by With use of the lowest-N KS orbitals thus obtained, we can define the Slater determinant Φ 0 (r 1 ,· · ·, r N ), with which the trial many-body ground-state wave function Φ(r 1 ,· · ·, r N ) for the VMC calculation can be constructed in the Slater-Jastrow type as 40 Φ(r 1 ,· · ·, r N ) = exp[J(r 1 ,· · ·, r N )]Φ 0 (r 1 ,· · ·, r N ), (12) where the Jastrow function J(r 1 ,· · ·, r N ) contains the terms to describe electron-nucleus correlation u 1 (r i ), two-electron correlation u 2 (r i − r j ), and three-body nucleus-two-electron correlation u 3 (r i , r j , r i −r j ) as The actual choice of the forms for u 1 (r i ), u 2 (r i −r j ), and u 3 (r i , r j , r i −r j ) as well as their optimization is done by adopting the CHAMP-code package 41 as it is. Then the expectation value A of an operator A is given by By putting A =n(r) = i δ(r − r i ) in Eq. (14), we obtain n(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 where E R is the reference energy to be adjusted to E 0 in the course of DMC simulations by removing the τdependence from the asymptotic form of Ψ(τ ) at τ → ∞.
Note that the formal solution to Eq. (15) is written as where {Ψ n } is the normalized mutually-orthogonal complete set of eigen functions for H with the corresponding set of eigen energies {E n }. Then, as long as Ψ 0 |Φ = 0, the asymptotic τ -independent wave function Ψ is reduced to the true ground-state wave function Ψ 0 , apart from the normalization 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 the presence 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 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 depend on N in the fixed boundary condition, leading to the hope that unphysical nodeposition dependent effects will be removed by extracting the N -independent results. In performing actual fixednode 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 second-order approximation to the exact expectation value 40 , which amounts to where n(r) VMC = n(r) VMC .

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 states in KS orbitals are continuum ones, for which we may write ε i = k 2 /2 with momentum k and φ iσ (r) = R klσ (r)Y lm (r/r) with the spherical harmonics Y lm (r/r) and the radial wave function R klσ (r) satisfying the following boundary condition at r (≡ |r|) → ∞ : (18) apart from a normalization factor, where j l (kr) and n l (kr) are the spherical Bessel functions and δ lσ (k) is the phase shift of angular 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 number of the bound states in KS orbitals N BS is given by lσ (2l + 1)δ lσ (0)/π, while Z = lσ (2l+1)δ lσ (k F )/π by the Friedel sum rule 44 .

E. Embedding energy
The embedding energy δE is defined as the difference of the ground-state energies between the atom-embedded EG and the 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 E Z a , we can obtain δE through Since δE is of the order O(1) in comparison with E 0 of the order 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) − n 0 ] and δ ′ lσ (k), the derivative of δ lσ (k) with respect to k, into the following form with ensuring the convergence of integrals: 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 density of states, δ ′ lσ (k)/π.

III. RESULTS FOR PROTON IMMERSION
A. LSDA for the bulk system In line with the previous result 45 , the ground state in the proton-embedded (Z = 1) bulk EG is perfectly spin-unpolarized, i.e., n ↑ (r) = n ↓ (r) at every r, at least for r s < 15 and our results on n(r) and δE in LSDA are in good agreement with those in previous calculations [12][13][14][15]18,19,46 . In Fig. 4, the obtained n(r) normalized by n(0) is plotted as a function of 2k F r/π for r s = 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 r s < 1.97, we obtain δ lσ (0) = 0 and thus N BS is zero, leading to the typical n(r) in CR with H + screened by metallic electrons in a short range. For r s ≥ 1.97, on the other hand, δ s (0) = π and thus N BS = 2, seemingly implying the sudden appearance of H − ion confined in EG at r s = 1.97 12 . This must be true, if n L (r) the localizedelectron distribution is about the same as n(r) for |r| smaller than the H − ion range, as is the case for r s 12.5, in which δ s (k) is distinctive and almost a universal function of k/k F , as seen in Fig. 4(b). Note that the deep dip in n(r) just outside the ion region, as seen for r s = 14 in Fig. 4(a), is a typical electron profile describing the repulsive action of localized closed-shell electrons to exclude metallic electrons from the ion region by the Pauli exclusion principle. For 1.97 < r s 12.5, however, the "localized" electrons behave much differently; at r s = 4, for example, they extend long up to |r| ≈ 82 and concomitantly n L (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, by remembering that KS orbitals in DFT are not necessarily tied with real physical entities but just introduced for mathematical convenience to correctly reproduce n(r), we can assume that this "local- In (a), the dotted 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 this situation.
ized" orbital does not represent a real localized state but just describes the long-range change of n(r) in the form of an envelope density over ξ K in SSR. Notice that δ s (k F ) in this density region lies between 0.7×(π/2) and 1.3×(π/2), which means that δ s (k F ) is close to π/2, a value expected in the Kondo resonance state in IAM 22 .

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 r s = 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 , though the results in VMC do not match so well, assuring the importance to perform DMC for taking the expectation values in accordance with Eq. (17).
As for N dependence or the size effect, we find that

C. Cusp theorem and the on-top density
According to the cusp theorem 47 , n(r) near a nucleus of atomic number Z behaves rigorously in the manner as where n cusp (r) exhibits strictly a linear change with r in semilog plots (the dashed lines in Fig. 6(a) and (b)). By exploiting this linear behavior, we can rather easily and accurately determine the on-top density n(0) in DMC by looking 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 the rapid increase in energy scale determined by the r −1 -Coulomb potential.
For r s 11, n(0) thus obtained is about the same as the exact one for an isolated H − ion 39 , as shown in Fig. 6(b) for r s = 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 at r 0.2 with the increase of r s ; for r s 12.5, the majority of data points deviate upward, while opposite is the case for r s 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 is consistent with the formation of the closed-shell ion, leading to the conclusion that the transition to the H − ion confinement state occurs at r s ≈ 12, 5, the same r s as that in LSDA.
In Fig. 7, we plot n(0) by changing r s and N and find that the results do not depend on N for r s 3, but they do for r s 3, i.e., in the CR-SSR transition region. The transition is signaled by a jump in n(0) in both LSDA and DMC, although its magnitude decreases with increasing N and eventually at N → ∞ no jump is seen in LSDA even at r s = 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 find that this jump becomes much enhanced at the magic  : (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, 170, respectively, reflecting the competition between the −r −1 and r 2 potentials in Eq. (4) and the outermost s 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 the actual SSR state in Kondo physics is situated at µ, extending very long over the range ξ K with an s-wave character. Then the jump 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 r s = 3.85, 2.65, and 2.87, respectively, in LSDA. In DMC, convergent results are also easily obtained for the same series of N and the transition occurs, for example, at r s ≈ 2.65 for N = 58 → 60, just as in LSDA. The next N in this 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, a change in the stacking sequence occurs in LSDA. Those transition points, along with the second transition point obtained in DMC by the data shown in Fig. 6(b), provide the ground-state diagram in (r s ,N −1 ) space in Fig. 3.

D. Embedding energy
In calculating δE by use of Eq. (19) at some finite N , we need E 0 (N −1, 0) in addition to E 0 (N, 1), but for an even number of N , the ground state of the (N−1)-electron system is necessarily a spin-polarized one. Nevertheless, for smoothly connecting to the spin-unpolarized ground state at N → ∞, it is better to calculate in the spinunpolarized situation. In LSDA, by a spin-symmetrized sum of KS orbitals by fractionally occupying the orbitals at the Fermi level, we can obtain this needed spinunpolarized ground-state energy E 0 (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 E 0 (N − 2, 0) in DMC, starting with the Slater determinant Φ 0 (r 1 , · · · , r N −2 ) without the outermost s orbital. Then we estimate E 0 (N −1, 0) by with the coefficient α N−1 determined through E 0 (N−1, 0) and E 0 (N−2, 0) in LSDA. We have deduced this approximation scheme by considering that the leading term in E 0 (N−1, 0) is, in general, in proportion to (N−1) due to 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 E 0 (N − 2, 0), we calculate E 0 (N, 0) in DMC, with which we estimate E 0 (N−1, 0) by a similar strategy.
In Fig. 8, δE is given as a function of r s , 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][13][14][15][16][17][18][19] , though 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 to play a key role in hydrogen storage in metals and hydrogen is expected to reside at a site with r s ≈ 4 in an inhomogeneous metal.

E. Anomalous Friedel oscillation
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 n N,Z (r) in the finite-N system corresponding to E 0 (N, Z) and then we estimate n(r) through the cancellation of the sphere-boundary effect by subtracting n N −1,0 (r) from it as n(r) ≈ n N,1 (r) − n N −1,0 (r) + n 0 .
In LSDA, we can employ Eq. (23) by using n N −1,0 (r) which is obtained in the spin-unpolarized situation simultaneously with E 0 (N−1, 0), but in DMC, only n N −2,0 (r) associated with E 0 (N − 2, 0) is available. Thus we estimate n N −1,0 (r) from this n N −2,0 (r) as where λ is a parameter to correct the difference in the sphere radius, defined as λ = [(N − 2)/(N − 1)] 1/3 , and the r-dependent exponent β(r) is determined with use of the data for the electron distribution in LSDA as An example of n(r) obtained through the above procedure at r s = 4 is given in Fig. 9, in which we find good convergence for r 7 in changing N by the comparison between N = 60 and 170, illustrating that this procedure doubles the size-convergent 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 and size 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 the outermost 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][24][25][26] Here only the s-wave contribution, which indeed dominates others in the present case, is considered and δ (0) s 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 can rewrite Eq. (26) into with the amplitude A(r) and the phase δ(r), given by If δ (0) s is in the range (0, π/2), δ(r) increases gradually from δ (0) s for small r to δ (0) s + π/2 for large r in accordance with the change of F (r/ξ K ). On the other hand, if δ (0) s is in the range (π/2, π), δ(r) decreases gradually from δ (0) s for small r to δ (0) s −π/2 for large r. With taking care of such a gradual change in δ(r), we can determine r i the ith zero of δn(r) [≡ n(r)−n 0 ] by 2k F r i +δ i = iπ+π/2 with δ i ≡ δ(r i ) and i = 1, 2, 3, · · · . Then the half period of the Friedel oscillation ∆ i in SSR is given by 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 effect is an important consequence of the presence of SSR, but this modulation effect has not been well recognized even in the Kondophysics community due to the fact that the effect is totally absent for δ (0) s being equal to a multiple of π/2 (see the denominator 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 2k F r/π at r s = 4 to check whether the modulation effect exists or not in our first-principles calculations. In DMC, the first four zeros of δn(r) in units of ∆ (0) , 2k F r i /π, are given by 0.994 ± 0.005, 1.970 ± 0.008, 2.949 ± 0.013, and 3.927 ± 0.019 with the errors estimated by the distribution of data points around δn(r) = 0 at each r i . 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 r i and r i+1 . Those results, at least definitely those for i = 1 and 2, assure the existence of the shortening of the Friedel-oscillation period, confirming the emergence of SSR.
Once the data for {r i } are known, we can independently calculate the modulation factors through Eq. (29) as a function of ξ K and δ (0) s with using F (x) in Fig. 10(a). The results are shown in Fig. 10(c), from which we find that our DMC data for the modulation factors (actually the shortening factors in this case) for i = 1 − 3 agree very well with those obtained at ξ K ≈ 35/k F = 73 and δ (0) s = 0.86 × (π/2) (see the horizontal and vertical double-dotted-dashed lines in Fig. 10(c)), indicating that by quantitatively analyzing first few modulation factors, we can determine both ξ K and δ (0) s uniquely, even if only the data 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 v F , the Kondo temperature T K is estimated by T K = v F /ξ K ≈ 0.0066 Hartree = 2, 100K with about 10% errors, revealing that this is an astonishingly high-T K system. With such a high T K , 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 .
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 from the "localized" KS state obtained at Z = 1, δnL(r)r 3 , is also given and seen as playing a role of an envelope density.
The shortening effect in the Friedel-oscillation period is also found in LSDA, as seen in Fig 11, in which δn(r)r 3 is plotted for both Z = 1 and −1 in the bulk system at r s = 4. For large enough r outside the SSR binding radius ξ K , the Friedel-oscillation 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 in which the anomalous Friedel oscillation is absent. In LSDA, the modulation of the Friedel oscillation is brought about by the contribution from the "localized" density n L (r). Incidentally, the oscillation-period shortening factors in LSDA are much smaller than those in DMC; (δ 2 −δ 1 )/π = 0.0133 and all others are less than 0.004, indicating δ (0) s ≈ 0.014×(π/2) in LSDA. (We note that δ s (k F ) in Fig. 4(b) is equal to δ (0) s +π/2 ≈ 0.507π.) We can easily understand the reason for this large difference in δ (0) s 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 non-interacting system is qualitatively different from the true correlated SSR wave function. Due to δ (0) s ≈ 0 in LSDA, we cannot employ the diagram in Fig. 10(c) to determine ξ K very accurately. Therefore we estimate ξ K from the extent of n L (r), which is 82, giving T K to be about 1,900K. In relation to n L (r), ε BS the binding energy for this "localized" KS state in LSDA is given as 0.0115 Hartree = 3,600K, which is about twice as large as T K . Thus, although ε BS has no direct physical meaning, this quantity seems to be a good measure for the magnitude of T K . In our calculations in LSDA, the values for ε BS are 97 K, 2,200K, 3,500K, 2,200K, and 720K for r s = 2.07, 3, 6, 8, and 10, respectively, indicating that in the majority of the SSR region, namely, for 3 r s 8, we may expect T K to be well beyond 1,000K.

IV. DISCUSSION
Five comments on the present work are in order: (i) In SSR, we find a good semi-quantitative 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 calculation of n(r).
(ii) From our present perspective, we may regard our previous study in LSDA on the spin-polarized ground states for second-period atoms in Periodic Table 45 as a successful extension to multi-channel Kondo systems 1 , in which the Hund's-rule coupling plays a crucial role in producing the spin-polarized ground states.
(iii) In a short term, the immediate next target of research is a hydrogen molecule H 2 immersed in EG 12 to pursue a new concept in chemical bonding 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 r s > 3.2. In a longer term, we can expect fruitful research on new aspects in the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction 1 and dilute magnetic semiconductors 52 .
(iv) If the number of immersed protons is increased up to a macroscopic level and those protons are arranged into a lattice in a metal with the SSR-density region, we may obtain a Kondo lattice with high T K . Furthermore, if the lattice constant is so arranged as to tune the RKKY interaction to be about the same as T K for realizing the quantum-critical situation in the Doniach phase diagram 53 , then as indicated in Fig. 1, we may obtain high-T c superconductivity in the spin-fluctuation mechanism in a metallic hydrogen alloy. In this regard, the transition-metal-hydride system such as TiH 2 and ZrH 2 54-56 , usually used for batteries, might be a promising candidate, although it seems that the metallic electron densities in the metal hydrides so far synthesized are too high for our purpose.
(v) The basic reason for about hundred-times increase of T K and consequently T c 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 T c over 100K has been discussed in solid hydrogen at pressures of about 500GPa in the conventional phonon mechanism 57-60 . Very recently, stimulated by the experimental result 61 , similar discussions are made on H 2 S 62,63 as well as H 3 S 64 at pressures of about 200GPa. Note that the metallic densities in those systems are found to be in the CR region, namely, r s < 2 (or typically r s ≈ 1.4). Therefore, our present proposal has nothing to do with those conventional theories for superconductivity. We also emphasize that high metallic densities realized at such very high pressures are not needed in our proposal, indicating that high-T c superconductivity is expected to occur even at ambient pressure, contrary to the case of solid H 2 or sulfur hydrides.

V. SUMMARY
In conclusion, by employing both DMC and LSDA, we have revealed the emergence of SSR in the protonembedded electron gas by the confirmation of an anomalous Friedel oscillation characteristic to the Kondo-like spin-singlet state with quantitatively determining T K 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, depending on the metallic electron density and Z, SSR takes the place of CR. This work also provides a first concrete example to show how the KS orbitals behave to represent the exact n(r) in strongly-correlated electron systems such as those describing Kondo physics, even though they are not always physically relevant. Finally, this work discloses a long-sought high-T K system, opening an unexpected route to room-temperature superconductivity in a Kondo lattice made of protons and motivating people engaged in hydrogen-based physics, chemistry, and technology to synthesize hydrogen alloys exhibiting high-T c superconductivity at ambient pressure.