Valence fluctuations and Kondo resonance in Co adatom on Cu2N/Cu(100) surface: DFT + ED study

We report density functional theory plus exact diagonalization of the multi-orbital Anderson impurity model calculations including the spin–orbit coupling for the Co adatom on the top of Cu2N/Cu(001) surface. It is found that experimentally observed zero-bias peak in differential conductance can occur for quasi-degenerate many-body solution representing a mixed valence state of the Co adatom with the non-integer d manifold occupation of 7.4. Estimated 5 meV inelastic step energy is in a quantitative agreement with the experimental scanning tunnelling spectroscopy value of 6 meV. We find the non-zero spin and orbital magnetic moments in the applied magnetic field, in a reasonable agreement with the x-ray magnetic circular dichroism data. The presence of large orbital magnetic moment was not considered in previous spin-only theoretical interpretations, and can unravel complexity of the physics behind the experimental spectra for Co and other adatoms on Cu2N/Cu(001) surface.


Introduction
Experimental and theoretical studies [1,2] of magnetic atoms adsorbed on non-magnetic surfaces provide a fundamental insights into the quantum many-body phenomena at the nanoscale. They imprint non-trivial signatures in scanning tunneling microscopy (STM) and x-ray magnetic circular dichroism (XMCD) measurements, and can serve as a prototype for potential applications as qubits for quantum information processing [3].
Surprising resistivity increase rather than expected decrease below 8 K observed by de Haas and van der Berg in the gold samples [4], and a resistivity minimum followed by a strong increase with reducing temperature in the dilutely doped non-magnetic metals [5] is a hallmark of the effect that is explained by Kondo [6,7] by the presence of magnetic impurities. This Kondo effect emerges from the exchange interaction between the sea of conduction electrons and the magnetic moment formed by localized d (or f) electrons of impurity. The most distinctive feature of the Kondo effect is the formation of so-called Abrikosov-Suhl [8][9][10] (or Kondo) resonance, the sharp resonance in the density of states (DOS) near the Fermi level (E f ). It is a commonly accepted viewpoint that this resonance has the many-body origin, and is not explained by the single-electron theory, like the density functional theory (DFT) [11].
Direct experimental observation of the Kondo resonance in dilutely doped metals using the photoemission is difficult as it requires measuring of the spectra very near E F with sufficient energy resolution. The more favorable approach is achieved by scanning tunneling spectroscopy (STS) measurements [12] for the magnetic impurities adsorbed on noble metal surfaces where the differential conductance G = dI/dV is used to probe the local impurity DOS near E f . The so-called 'zero-bias' anomalies are experimentally observed and associated with the Kondo resonance. Very recently, alternative interpretation was proposed [13] based on the spin-polarized time-dependent DFT in conjunction with many-body perturbation theory.
In order to reduce the strong hybridization between magnetic adatom and the underlying metallic substrate the ultra-thin insulating layer was introduced [15] where the Co adatom was deposited on the 'insulating' Cu 2 N layer on the top of Cu(100) substrate. Presumably, the direct tunneling from the tip into the substrate is suppressed, and the Fano interference [14] plays no role. The zero-bias anomaly in G takes a form of a Lorentzian resonance [15]. Simultaneously, the steps in differential conductance at finite bias voltages are observed signaling the spin-flip excitations caused by the magnetic anisotropy (MA) caused by the crystal-field (CF) and spin-orbit coupling (SOC) acting on the d-shell of the Co adatom. This interpretation was supported by the numerical renormalization group calculations for the S = 3/2 Kondo impurity model [16].
Theoretical electronic structure of the Co atom on the top of Cu 2 N/Cu(001) was considered by the one-crossing-approximation (OCA) [17] to the multiorbital single impurity Anderson model [18] (SIAM) augmented by the ad hoc MA term, and assuming the localized spin S = 3/2. The OCA model describes very well experimentally observed STS spectra [15,17]. Nevertheless, the model parameters choice remains purely phenomenological.
The XMCD experiments [19] for Co@Cu 2 N/Cu(001) show large orbital contributions to the magnetic moment of Co adatom. This orbital moment is not accounted in S = 3/2 models [15,17]. In this work, we approach Co@Cu 2 N/Cu(001) making use of the combination of DFT with the exact diagonalization of multiorbital SIAM (DFT + ED). Previously, there were already attempts to explore this direction of the research for Co in bulk Cu, and Co on the Cu surface [20], and the Co atom on the top of Cu 2 N/Cu(001) [21] neglecting the SOC. In this work, we extend the DFT + ED approach by including the SOC, in order to elucidate the spin and orbital state of the Co-adatom.

Methods
Consider the effective SIAM (see book of Hewson [18]) The operator d † mσ creates an electron with the magnetic quantum number m and spin σ = {↑, ↓} in the impurity d shell with energy d . The operator b † kσ creates an electron in the conduction bands of the substrate with the energy k . The hybridization V mk parameters describe the coupling of substrate to the impurity orbitals. The last term in equation (1) represents the Coulomb interaction. Parameters of the Coulomb interaction are taken as external parameters of the model and their particular choice will be discussed below.
The equation (1) can be rewritten [22] using the energy-dependent states |b( , m, . The Hamiltonian equation (1) transforms into: where n b The impurity-level position d which yield the desired n d , and the bath energies km are measured from the chemical potential μ, that was set to zero. The SOC ξ parameter specifies the strength of the SOC, whereas Δ CF matrix describes CF acting on the impurity.

DFT
The DFT calculation were performed on a supercell of three Cu(100) layers, a Cu 2 N-(2 × 2) monolayer, and the Co adatom followed by four empty Cu layers modeling the vacuum. Figure 1 shows ball model of the Co@[Cu 8 N 4 /3Cu 8 ] supercell employed for the adsorbate atop of Cu. The structure relaxation was performed employing the VASP method [23] together with the generalized gradient approximation (GGA) to spin-polarized DFT without SOC. The adatom-substrate distance as well as the atomic positions within the Cu 2 N layer, and one Cu(100) layer underneath were allowed to relax.
The occupation n d , spin M S , and orbital M L magnetic moments for d-states of Co adatom calculated with scalar-relativistic VASP (GGA), the full-potential linearized augmented plane wave method (FLAPW) [24,25] and local spin-density approximation (LSDA), FLAPW with SOC (LSDA(SOC)) [26], and the rotationally-invariant relativistic LSDA(SOC)+U method [27] are given in table 1. In these calculations magnetization is directed along the z-axis perpendicular to the surface. The SOC is included in a self-consistent second-variational procedure. The values for the Coulomb U = 4 eV and exchange J = 0.9 eV for Co were used [20] which are in the ballpark of commonly accepted U and J for transitional 3d-metals, were used in LSDA + U calculations.
The number of Co-d-electrons (n d occupation) in the Wigner-Seitz sphere (R WS , VASP) and the 'muffin-tin' sphere (MT, FLAPW) are shown in table 1 together with the values of spin and orbital magnetic moments on the Co atom d shell. Since the Co-atom R WS = 2.5 a.u. is slightly bigger than the R MT = 2.2 a.u., the n d occupation obtained in VASP slightly exceeds their FLAPW value. The corresponding M S are in reasonable agreement. When SOC is allowed, the non-zero M L is obtained which is more that a twice bigger than the orbital moment for the bulk Co (0.07 μ B ) obtained in LSDA. Coulomb-U does not significantly increase the M S , but has a strong effect on the M L due to additional orbital polarisation induced by non-spherically symmetric contributions from LDA + U functional. All of the listed static mean-field results suggest non-integer occupation of Co with n d > 7 typically found in the 3d metallic systems.
Following the recipes of the dynamical mean-field theory (DMFT) [28,29] we make use of the DFT (LDA) local Green's function G 0 (z) calculated with help of the FLAPW [24,25], in order to define the parameters for the equation (3). Here, the energy z is counted from the Fermi energy E F , and the index γ ≡ lmσ marks the d-orbitals in the MT-sphere of the Co adatom. Note that we make use of the LDA + U projectors from the Kohn-Sham DFT solutions into the local MT-sphere as described in reference [30], and the non-spin-polarised LDA is used to extract the hybridization function.
With the specific choice of the Cartesian reference frame (see figure 1), the local Green's function G 0 (z) becomes diagonal in the basis of cubic harmonics m = {xz, yz, xy, x 2 − y 2 , 3z 2 − r 2 }. Moreover, it is convenient to use the imaginary energy axis over the Matsubara frequencies iω n . The corresponding non-interacting Green's function of the equation (3) will then become Thus, the hybridization function equation (5) can be evaluated making use of the local Green's function G 0 (z).

Exact diagonalization
The discrete bath model is built by finding bath energies and amplitudes which reproduce the continuous hybridization function as closely as possible.
The fitting is done by minimizing the residual function, using the limited-memory, bounded Broyden-Fletcher-Goldfarb-Shanno method [31,32], with the parameters km and V km as variables. The factor 1 ω γ n with γ = 0.5 is used to attenuate the significance of the higher frequencies. The fitted bath parameters are shown in table 2. These parameters are used to build the SIAM Hamiltonian equation (3). The corresponding LDA Δ m together with fittedΔ m are shown in figure 2.
The SOC parameter ξ = 0.079 eV is taken from LDA calculations. The Slater integrals F 0 = 4.00 eV, F 2 = 7.75 eV, and F 4 = 4.85 eV are used for the Coulomb interaction. They correspond to the same values of U and J used in the LDA + U calculations (see table 1).
The model equation (3) is solved using exact diagonalization in a reduced many-body state space [33]. In order to minimise the computational effort a reduced many-body basis [34] is introduced with the cutoff on the number of excitations N exc for each N-electron sector of Hilbert space H N ,  Table 3. The occupation N , fluctuation ( N 2 − N 2 ) 1 2 , n d , spin, orbital and total moments of the impurity d-shell as a function of μ. Grand-canonical averages are at low temperature where n is a number of d-shell electrons, p bath electrons above Fermi level and q bath holes under Fermi level. N < b is the number of bath orbitals located below the Fermi level (in the present case, N < b = 20). The Hamiltonian equation (3) was diagonalized using the Lanczos method, with 100 iterations and by calculating the 20 lowest eigenstates. The excitation cutoff N exc = 4 was used.

Results and discussion
The total number of electrons N, and the d-shell occupation are controlled by the d parameter. It has a meaning of the chemical potential μ = − d in equation (3). In DMFT it is quite common to use μ = V dc , the spherically-symmetric double-counting which has a meaning of the mean-field Coulomb energy of the d-shell, and to use standard (AMF) V dc = (U/2n d + 2l 2(2l+1) (U − J) n d ) [35] form, or the fully localized limit V dc = (U − J)/2(n d − 1) [36]. Since precise definition of n d depends on the choice of the localized basis, we adopt a strategy of reference [20], and consider a value of μ as a parameter.
The the expectation values of total number of electrons (d-shell + bath) N , charge fluctuation ( N 2 − N 2 ) 1 2 , the d-shell occupation n d , as well as the expectation values of spin (S), orbital (L) and total spin-orbital (J) moments as grand-canonical averages at low temperature k B T = β −1 = (1/500) eV (20 K) are shown in table 3. Note that spin S, orbital L and total J moment 'quantum numbers' are found using X 2 = X(X + 1) for X = S, L, J where, The spectral densities of d-states (DOS), are shown in figure 3. It is seen that for μ = 25 eV, and the occupation N is integer, and there are no charge fluctuations. Corresponding DOS has no peak near E F (see figure 3(A)). With an increase of μ to 26 eV, charge fluctuations become non-zero. Simultaneously, DOS has a peak in a close vicinity of E F (see figure 3(B)). Note that the peak at E F occurs only in a presence of SOC. Further increase of μ to 27 eV suppresses the charge fluctuations, and moves the DOS peak position away from a vicinity of E F as shown in figure 3(C). Now we analyze the GS and low-energy excitations. Eigenstates |Ω for different μ are shown in table 4. Each |Ω corresponds to an integer N occupation (d-shell + bath) sinceN commutes with Hamiltonian equation (3). The probabilities to find the atomic eigenstates |n with integer occupation d n , P n = n|Ω Ω|n are shown in table 4, and define the occupation n d = n P n n d n .
For μ = 25 eV the |Ω (N = 27) doublet is a GS for which we obtain P d 7 = 0.58 and P d 8 = 0.23. This is predominantly d 7 state with non-integer valence (n d = 7.07) due to hybridization with the substrate. Similarly, for μ = 27 the |Ω (N = 28) singlet is a GS with P d 7 = 0.32 and P d 8 = 0.55. This is predominantly d 8 state with non-integer valence n d of 7.70. For both, μ = 25 eV and μ = 27 eV, the DOS (see figures 3(A) and (C)) does not have a peak at E F . The first DOS peak is shifted by ∼0.1 eV above E F for μ = 25 eV which corresponds to transition from |Ω (N = 27) to |Ω (N = 28). For the μ = 27 eV, there is a peak at ∼0.1 eV above E F (transition from |Ω (N = 28) to |Ω (N = 29)), and a peak at ∼0. 15 eV below E f due to |Ω (N = 28)-to-|Ω (N = 27) transition. Corresponding DOS very near E F has a 'Kondo'-like resonance. We note a similarity to the intermediate  valence [37] in the rare-earth compound SmB 6 , where the multi-orbital 'Kondo-like' singlet ground state is formed due to the valence fluctuations between f 5 and f 6 states [38] resulting in a 'Kondo'-like resonance at E F . The low temperature STM experiments [15,17] have observed the peak in differential conductance G(V) of Co@Cu 2 N/Cu(001) at zero bias. This peak can be identified as the elastic contribution to G(V), and is consistent with the DOS peak at E F for μ = 26 eV shown in figure 3(B). In reference [15] this peak was interpreted as a Kondo-resonance for S = 3 2 spin with positive MA D > 0. This interpretation was supported by the OCA calculations [17], in which the effects of the CF and the SOC were included as a model MA term equation (9). Our calculations demonstrate that this peak can be explained differently by charge fluctuations between |Ω (N = 27)d 7 -like and |Ω (N = 28)d 8 -like states.
The expectation values of the z-axis projections of the total Ω|J z |Ω , orbital Ω|L z |Ω , and spin Ω|S z |Ω angular momenta for GS and low-energy excitation energies for different values of μ are shown in table 5. The states in table 5 approximately correspond to the eigenstates of the effective Hamiltonian equation (9) with the uniaxial D and transversal E magnetic anisotropies,    For the half-integer J = 3/2 (see figure 5(A)) the GS quadruplet splits into two |m = ± 1 2 and |m = ± 3 2 doublets [39]. For D > 0, the |m = ± 1 2 becomes a GS. The in-plane anisotropy E mixes the GS with exited doublet. However the double degeneracy of resulting GS and exited states remains. From the data shown in table 5, we estimate the uniaxial anisotropy D = 9 meV, and the in-plane anisotropy E = 0.0 meV.
For the integer effective J = 1 (see figure 5(B)) the anisotropy D lifts the degeneracy of (|m = −1 , |m = 0 , |m = 1 ) [39]. In case of D > 0, the |m = 0 state becomes GS, and the states m = ±1 become an excited doublet. The in-plane anisotropy E allows for a quantum tunneling between these doublet states, so that they become a linear combinations |χ ± ∼ |m = −1 ± |m = 1 with corresponding energies D ± E. From the data shown in table 5 for μ = 26 eV, we obtain D = 4.85 meV and E = 0.05 meV.
The low-energy inelastic contributions to the tunneling can be easily estimated making use of the data in table 5. Considering the μ = 26 eV as a non-coherent statistical mixture of effective total moments J = 3 2 (probability α = 0.15) and J = 1 (probability 1 − α = 0.85) states. For J = 3 2 we obtain a single inelastic step in differential conductance G(V) at Δ 0 = 9 meV shown in figure 5(A). For J = 1, we get two steps, Δ 1 = 4.8 meV and Δ 2 = 4.9 meV (see figure 5(B)) which are practically coincide within the experimental resolution.
We calculated inelastic tunneling spectrum making use of reference [40]. The inelastic contribution to the differential conductance G(V) at T = 2.5 K are shown in figure 5(C). When mixed with the 15/85 percentage they would give two-step conductance with the main step at 5 meV (close to the experimental 6 meV step [15]), and a shallow second step at 9 meV. Either this second step cannot be resolved in the experiments of references [15,17] or (most likely) it is not there. This contradicts to the experimental results showing the only one inelastic step in the conductance. Also, we note that that reference [17] claims an increase of the step location up to 10 meV depending on the position of the Co adatom with respect to the edge of the Cu 2 N island on Cu(001) surface.
This partial contradiction between our theory and the experimental observations brings about the important issue: to what extend the discrete ED-model with finite (and integer) N-number of particles describes the experimentally measured open quantum system with an infinite N. The reason to introduce the mixed-state solution is that the energy minimum (figure 4) corresponds to the non-integer N, and no eigenstate of equation (3) exists. In the quantum mechanics, the mixed state solution can occur when the Hamiltonian does not include all the information needed for the complete description of the physical system. It is likely that we deal with such a case, and our finite bath model equation (3) needs to be extended to the continuum bath with N = ∞. Nevertheless, our calculations allow to speculate that when going to N = ∞, the mixed state will transform to the coherent superposition of J = 3/2-like and J = 1-like states for which the ISTS will produce the single step only. Now we make a comparison with available XMCD results [19]. The XMCD spectra are obtained in the external magnetic field of 5 T at low 8 K temperature. When the B = 10 T (0.6 meV) along the z-axis (comparable with experimental value) is used in equation (3), we obtain nearly paramagnetic solution with tiny spin and orbital magnetic moments, mainly due to the spin-flip terms included in the full Coulomb interaction matrix [21]. In order to obtain the magnetic solution with the spin magnetic moment corresponding to the experimental XMCD data, we used larger B-field of 20 meV (≈300 T) [21]. The results are shown in table 6. For μ = 26 eV, which corresponds to the paramagnetic ground state, we obtain smaller n d occupation of 7.3. Slight increase of the μ, in order to keep the same n d = 7.4 as in paramagnetic calculations, yields the spin, orbital and dipole [41] moments, and their R LS ratio in a reasonable quantitative agreement with the XMCD experimental data (see table 6). The calculated spectral DOS for B = 20 meV case is shown in figure 6 in comparison with the paramagnetic DOS (B = 0). It is seen that the DOS peak near E f for B = 0 is split in two by the magnetic field, in qualitative agreement with STM spectra [15].

Summary
We report the many-body calculations of the Co adatom on the top of Cu 2 N/Cu(001) within the multi-orbital Anderson impurity model. The exact diagonalization is performed for discrete bath model containing forty bath orbitals at low temperature (β = 500 eV −1 ) including SOC. DFT calculations were used to define the input parameters for AIM.
Considering the electronic structure as a function of the localized d-level position d in equation (3) (= −μ, the AIM chemical potential), we found that the experimental 'Kondo' peak can occur for quasi-degenerate solution of AIM only. We propose a scenario that the Co@Cu 2 N/Cu(001) is in the mixed valence state fluctuating between d 7 , J = 3/2-like and d 8 , J = 1-like states with non-zero spin and orbital angular moments, and the non-integer d manifold occupation of 7.4. The calculated peak very near E F is in agreement with experimental STM spectra [15,17]. The estimated 5 meV energy position of the inelastic step is in a quantitative agreement with the experimental value of 6 meV [15].
Moreover, we estimate non-zero spin M s = 1.58 μ B and orbital M L = 0.71 μ B magnetic moments for the Co d-states in the applied magnetic field. They are in reasonable agreement with experimental XMCD results [19] M s = 1.4 ± 0.1 μ B , and M L = 0.9 ± 0.2 μ B . The presence of this large M L was not considered in previously reported spin-only theoretical interpretations for Co@Cu 2 N/Cu(001) [16,17]. It can unravel complexity of the physics behind the experimental STM spectra for Co adatom on Cu 2 N/Cu(001) surface. The results presented here can be important to explore the role of valence fluctuations and non-zero orbital magnetic moments for other adatoms, and push for revision of Kondo physics in the multiorbital systems.