M Smogunov

E Fabrizio

Tosatti

P P Baruselli

R Requist

A Smogunov

M Fabrizio

E Tosatti

Co adatoms on Cu surfaces: Ballistic conductance and Kondo temperature

Keywords: number(s): 73.63.Rt, 73.23.Ad, 72.15.Qm, 74.55.+v

Co adatoms on Cu surfaces: ballistic conductance and

Kondo temperature

I. INTRODUCTION

Since the observation of zero-bias anomalies (ZBA's) for Ce adatoms on silver [1] and Co adatoms on gold [2] by scanning tunneling microscopy (STM), the Kondo effect [3] of magnetic adatoms has become a subject of extremely high interest. STM measurements of the Kondo effect offer the possibility of achieving exquisite external control over a paradigmatic strongly correlated system. Several of the established Kondo systems have been reviewed in Ref. [4]. Despite the apparent simplicity of these systems, a full theoretical description is still lacking. Indeed, great efforts using both ab initio and many-body approaches have been made, but several open issues still exist.

The Kondo effect in magnetic adatoms has been successfully treated within an Anderson model approach [5] (mostly with a single impurity orbital of d z 2 symmetry; only recently has the whole d shell [6,7] been taken into account), but the role of the tip is still a subject of debate. STM measurements have conventionally been performed in the tunneling regime where the tip does not affect the results. Recent works [8][9][10][11] have looked beyond tunneling measurements to explore the contact regime, where the geometric details of the tip and its position above the adsorbate can affect the Kondo ZBA. Ab initio calculations including tip-induced perturbations are needed to describe, for instance, the observed progression of the Kondo temperature as a function of tip height above a Co adatom on Cu(100) [8][9][10].

STM conductance is usually calculated in the Tersoff-Hamann model [12,13], in which the tip is described by a single s orbital and current flows thanks to the coupling between the tip and nearby metallic states. However, the situation is more complicated in the presence of an adsorbate. The interference between tunneling directly into the surface and tunneling via the adsorbate causes a ZBA in the STM conductance

G(V ) = G back + G q 2 + 2qV /(k B T K ) -1 [q 2 + 1][(V /k B T K ) 2 + 1] (1)
with a characteristic Fano line shape [14], the fingerprint of the Kondo effect; here we introduce a 1/(q 2 + 1) factor in such a way that G/G back represents the signal to background ratio, which is experimentally found to be on the order of 10-30%. The parameter q describes the shape of the ZBA (q = 0 corresponds to a minimum, q = ±∞ to a maximum, while intermediate values give rise to asymmetric line shapes), while T K , the Kondo temperature, is proportional to its width (k B is the Boltzmann constant). The full STM-adatom geometry can alternatively be viewed as a nanocontact and the electrical conductance calculated within the Landauer-Büttiker formalism [START_REF] Datta | Electronic Transport in Mesoscopic Systems[END_REF]. However, when Kondo correlations are present the ballistic scattering matrix cannot be obtained from density functional theory (DFT) in the standard semilocal approximations because the latter do not include the many-body correlations responsible for ZBA's. DFT calculations are nevertheless indispensable in singling out the relevant adsorption geometries and electronic degrees of freedom.

To model Kondo ZBA's from first principles, recent works [START_REF] Lucignano | [END_REF][START_REF] Requist | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Baruselli | [END_REF] have developed a scheme to quantitatively join DFT and many-body calculations via an intermediate Anderson impurity model (AIM). The model parameters are determined by matching the mean-field scattering matrix of the AIM to the ballistic scattering matrix of a spin-polarized DFT calculation. Observables are then obtained by solving the AIM with the many-body numerical renormalization-group (NRG) method. This DFT + NRG scheme has successfully predicted the Kondo ZBA of nitric oxide adsorbed on the Au(111) surface, albeit underestimating the experimental Kondo temperature [START_REF] Requist | Proc. Natl. Acad. Sci. USA[END_REF]. Predictions have also been made for Co impurities and vacancies in carbon nanotubes [START_REF] Baruselli | [END_REF][19][20]. Different approaches to adding dynamical many-body correlations to first-principles transport calculations have been proposed [7,21].

Two main approaches have been adopted to calculate the Fano parameter q in Eq. ( 1).

The first is the "two-path model" [8,[22][23][24][25][26], in which the tip is coupled to the adatom d z 2 orbital via hopping t pd and to the surface via t pk , where k denotes conduction states of the surface and p denotes the states of the tip, giving the expression

q = t pd + k t pk V kd G k k t pk V kd G k → t 2 + t 1 V d G(0) t 1 V d G(0) , ( 2
)
where V kd are matrix elements between the d orbital and the surface and G k is the Green's function of the clean surface. By considering t pk , t pd , and V kd to be energy and momentum independent, and calling them, respectively, t 1 , t 2 , and V d , the right-hand side of Eq. ( 2) is obtained, where only the surface Green's function G(0) at the Fermi energy appears.

The second approach consists in neglecting the coupling of the tip to d states, i.e., t 2 ≡ t pd = 0 [27,28], so that just the density of states of the metal is probed by STM via t 1 ≡ t pk , assumed for simplicity to be momentum independent; this leads to the Fano parameter

q = k t pk V kd G k k t pk V kd G k → G(0) G(0) , (3)
which is zero for a particle-hole symmetric band if the momentum dependence of the matrix elements can be neglected. This second approach is justified for Co/Au by the fact that the experimentally observed Fano resonance is independent of tip height [24]. If t 2 were not negligible, q would be expected to have non-negligible z dependence because t 1 and t 2 generally have different dependence on the tip position R = (R,ϕ,z) (z being the height above the adatom and R the lateral displacement in the angular direction ϕ). The R dependence of q has been observed, but it is mainly a consequence of probing variations of the surface Green's function at different positions.

In this paper, we use our DFT + NRG scheme to study how the Kondo temperature and Fano line shape are affected by the location of the tip in two experimentally wellcharacterized cases: single Co impurities on Cu(100) and Cu(111) surfaces [8][9][10][11]. Experimentally, it is found that on the Cu(100) surface Co adatoms show a ZBA with T K = 88 K and q = 1.13; upon moving the tip laterally q decreases down to 0.6 [8]. When the tip approaches the Co adatom, T K increases to 700 K and q increases to ∼70 in one experiment [10] and T K =150 K and q ∼ 2 in another [9]. On the Cu(111) surface, T K = 54 K and q = 0.18; upon moving the tip laterally, q decreases to approximately zero [8]; when the tip approaches the adatom, both T K and q are unaffected [11].

Our main results can be summarized as follows. Spinpolarized DFT calculations show that on both surfaces the spin state of Co is S = 1, each of two magnetic d orbitals contributing approximately one Bohr magneton. On the (100) surface these orbitals are inequivalent; the d z 2 orbital is found to have a much higher Kondo temperature than the d x 2 -y 2 orbital. The effect of the tip is to increase the hybridization of orbital d z 2 by pushing the adatom into the surface, as well as itself providing another source of hybridization. On the (111) surface, the two magnetic orbitals are degenerate, and their hybridization does not increase as the STM tip gets closer, due to symmetry and structural reasons. Consequently, the Kondo temperature does not vary appreciably all the way from the tunneling regime to the contact regime. On both surfaces a precise determination of the Kondo temperature, which depends exponentially on the parameters of the AIM, is beyond the capabilities of our method. Both the Kondo temperature and the Fano parameter q are affected by numerous fine details, comprising the electronic structure of the surface and adatom; the details of the tip-adatom-surface nanocontact, in particular, how strongly symmetry is broken by the tip; and possibly nonequilibrium and multiorbital (beyond 2) effects.

The paper is structured as follows. In Sec. II we describe our method; in Secs. III and IV we present our results for Co/Cu(100) and Co/Cu(111), respectively; and in Sec. V we discuss the conclusions of our work.

II. OVERVIEW OF THE METHOD

We have employed, with a few simplifications, the method presented in Refs. [START_REF] Lucignano | [END_REF][START_REF] Requist | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Baruselli | [END_REF][19][20], to which we refer for further details.

First we perform a self-consistent, fully relaxed calculation of the electronic properties of the scattering region (as shown in Figs. 1 and3) by density functional theory. This is constituted by a 3 × 3 Cu supercell in the xy plane with a Co coverage of 1/9; in the z direction we use eight Cu layers plus a "pyramid" of five Cu atoms to simulate the STM tip for the (100) surface; for the (111) surface, we use seven layers and a four-atom pyramid. At this coverage the interaction between periodic replicas of the adatom is small. The calculations are carried out with the standard planewave package QUANTUM ESPRESSO [29] using the generalized gradient approximation (GGA) to the exchange-correlation functional in the parametrization of Perdew, Burke, and Ernzerhof [30]. The plane-wave cutoffs are 30 and 300 Ry for the wave functions and charge density, respectively. Integration over the two-dimensional Brillouin zone is accomplished using a 6 × 6 grid of k points and a smearing parameter of 10 mRy.

After obtaining the self-consistent electronic structure, the conductance in the z direction is calculated using the PWCOND routine [31]. Scattering eigenchannels and eigenvalues depend on the transverse momentum k xy ; hence, we restrict our conductance calculations to the single most representative k xy points: B = π L ( 1 2 , 1 2 ) for the (100) surface, and K = π L ( 2 3 ,0) for the (111) surface. This procedure introduces small systematic errors in the estimation of parameters but has the advantage of keeping the computational effort low. We verified that the error with respect to a more accurate calculation with 5 × 5 k xy points is less than a few percent. In the above expressions, L = 7.77 Å is the length of the supercell in the x and y directions, set to three times the equilibrium nearest-neighbor distance for bulk Cu for our pseudopotential, 2.59 Å, slightly larger than the experimental value 2.56 Å.

With the knowledge of the scattering eigenvalues t n , it is possible to compute the energy-dependent Fano factor [32]:

F = n t n (1 -t n ) n t n , (4)
which is experimentally accessible through noise measurements [33,34].

In the final step, an AIM is built in such a way that it reproduces the DFT scattering phase shifts as closely as possible when solved in the Hartree-Fock (HF) approximation. The Kondo temperature can be estimated after the AIM is solved by numerical renormalization group (NRG) [35,36]. In contrast to Refs. [START_REF] Lucignano | [END_REF][START_REF] Baruselli | [END_REF] where a phase-shift analysis was performed to determine the Fano parameter q, here, due to the additional complication of the dependence on transverse momentum and lack of even or odd symmetry along the z direction, the line shape is inferred directly at the DFT level by looking at the energy-dependent transmission eigenvalues, and fitting them with a Fano line shape.

This procedure is repeated for a series of tip-surface distances to show how the Kondo temperature and line shape vary in going from the tunneling to contact regime.

III. Co/Cu(100)

In this section we present results for the Co/Cu(100) system, which is found experimentally to have T K = 88 K and q = 1.13 in the tunneling regime [8]. These values are modified in the contact regime, where T K grows up to 700 K and q grows up to ∼70, [10]. Reference [9] found that q does not go beyond 2-a discrepancy which is probably due to the different nanocontact geometries in different experimental conditions.

A. DFT results

We find that the Co adatom adsorbs in the hollow position, with four nearest-neighbor Cu atoms. In this configuration the symmetry group is C 4v , and Co 3d orbitals are split into three singlets (d z 2 with symmetry A1, d x 2 -y 2 with symmetry B1, d xy with symmetry B2) and a doublet (d xz and d yz with symmetry E); the 4s orbital has A1 symmetry. In our scattering geometry the tip is built so as not to lower the C 4v symmetry of the adatom plus surface system. We discuss possible consequences of this approximation later, since symmetry is not preserved in real experiments.

GGA calculations, both with and without the tip, show the presence of two magnetic orbitals, d z 2 and d x 2 -y 2 , in agreement with Ref. [37], while the 4s orbital is highly hybridized, almost spin unpolarized and singly occupied; the electronic configuration is thus 3d 8 4s 1 , and the magnetic moment is close to 2μ B , in agreement with, for example, Ref. [38]. Table I reports some structural and electronic data for this system at different tip-surface separations. Structural data are in good agreement with Ref. [39]. Table I covers the approximately-known range of experimental tip-surface separations. The tip-surface separation could not be further reduced because already at the smallest value reported in Table I, d tip-sur = 4.12 Å, z 2 and U z 2 are affected by large errors. The source of error is the breakdown of our procedure for estimating parameters when the occupation of the d z 2 orbital grows significantly greater than 1 and part of the polarization is transferred to the d xz and d yz orbitals. This is expected since with increased coordination Co approaches the configuration it has as a bulk impurity, where the polarization is shared by all d orbitals [40].

Figure 2 shows representative density of states and transmission eigenvalues for d tip-sur = 5.58 Å. The interference between the s and d z 2 orbitals in the down-spin A1 channel gives a Fano line shape with q ∼ 1 centered at an energy around ↓ d z 2 ∼ 0.5 eV; all other channels give a negligible contribution to the total conductance.

B. Anderson model

In building an effective AIM, only the d z 2 and d x 2 -y 2 magnetic orbitals are retained, each coupled to a linear combination of conduction states with the same symmetry, A1 and B1, respectively. Within this approximation, the adatomsurface Hamiltonian is

H ad-sur = kσ i = A1,B1 ki c † kiσ c kiσ + kσ i = A1,B1 V ki (c † kiσ d iσ + d † iσ c kiσ ) + i=A1,B1 ( i n i + U i n ↑ i n ↓ i ) + U 12 n z 2 n x 2 -y 2 -J S z 2 • S x 2 -y 2 , ( 5
)
where we have introduced the single-particle energies ki and fermionic operators c kiσ and c † kiσ associated with conduction states with momentum k, symmetry i, and spin σ , the on-site energies i and fermionic operators d iσ and d † iσ associated with impurity orbitals (d A1 ≡ d z 2 , d B1 ≡ d x 2 -y 2 ), the hopping elements V ki between a conduction state and d orbital, the on-site Hubbard repulsion U i for d z 2 and d x 2 -y 2 orbitals, the TABLE I. Parameters for Co/Cu(100) at different tip-surface separations d tip-sur : the distance d tip-Co between the tip and the Co adatom, d Co-sur between the Co adatom and the surface, d Co-NN between the Co adatom and its four nearest neighbors, the total magnetization m of the unit cell in μ B , the on-site energies z 2 and x 2 -y 2 , the Hubbard repulsions U z 2 and U x 2 -y 2 , the hybridizations z 2 and x 2 -y 2 , the interorbital Hubbard repulsion U 12 , the Hund exchange J , the calculated Kondo temperatures T K,z 2 and T K,x 2 -y 2 in K, the Fano parameter q z 2 , the DFT conductance g = G/G 0 = g ↑ + g ↓ , and Fano factor F at the Fermi energy; distances are in Å and energies are in eV. interorbital Hubbard repulsion U 12 between d z 2 and d x 2 -y 2 , and the Hund exchange J > 0 between d z 2 and d x 2 -y 2 . We denote

d tip-sur d tip-Co d Co-sur d Co-NN m z 2 U z 2 z 2 x 2 -y 2 U x 2 -y 2 x 2 -y 2 U 12 J T K,z 2 T K,x 2 -y 2 q z 2 g F 7.
n σ i = c † iσ c iσ , n i = σ n σ i , S i = 1 2 μν d † iμ σ μν d iν ,
where σ is the vector of Pauli matrices σ = (σ x ,σ y ,σ z ). A generalized Tersoff-Hamann model is used for the interaction of the tip with the surface and the d orbitals on the adatom:

H tip = pσ i = A1,B1 pi c † piσ c piσ + pσ i = A1,B1 t 2i (c † piσ d iσ + d † iσ c piσ ) + pkσ i = A1,B1 t 1i (c † piσ c kiσ + c † kiσ c piσ ), (6)
where pi and c piσ denote the single-particle energies and destruction operators associated with conduction states of the tip with momentum p, symmetry i, and spin σ . The usual approach only takes into account the A1 symmetry channel, i.e., the tip-surface hopping t 1A1 and t 2A1 , which is motivated by the fact the apex atom of the tip (Cu in this case) has a single s orbital at the Fermi energy. In contrast, we allow for conductance through channels with different symmetry. Since our geometry preserves the symmetry of the Co d states, the total conductance g tot ≡ G tot /G 0 , where G 0 = e 2 /h is the quantum of conductance, is the sum of the conductance from A1 and B1 channels (we ignore all other symmetry channels as they do not contribute to the ZBA):

g tot = i=A1,A2,B1,B2,E g i g A1 + g B1 . ( 7
)
We are aware that symmetry is not preserved in real experiments, because the tip cannot be expected to have the ideal pyramidlike shape we have assumed; however, it is a good starting point for studying the problem. Should symmetry be broken, channels with different symmetries would start interfering, leading to a modified line shape. However, we believe the disagreement we find with the experimentally determined q is only partially due to this approximation: since a single Kondo temperature is relevant (see next subsection), the main effect of interference would be to slightly modify the hopping parameters from the tip to the surface and the adsorbate, only weakly affecting our estimate of q.

We now introduce the hybridization functions s i ( ) due to the coupling with the surface, t i ( ) due to the coupling with the tip, and the total i ( ):

s i ( ) = π k δ( -ki )V 2 ki → πρ si V 2 i , ( 8
) t i ( ) = π k δ( -ki )t 2 ki → πρ ti t 2 2i , ( 9
) i ( ) = s i ( ) + t i ( ), i = d z 2 ,d x 2 -y 2 , ( 10
)
where the expressions πρ s V 2 i and πρ t t 2 2i are valid as long as we assume energy independent quantities (ρ si and ρ ti denote the density of states of symmetry i at the Fermi energy for the surface and the tip, respectively). Our method only allows us to infer the total i ( ), which we assume to be energy independent, and simply call i . Also, we assume that the density of states of conduction electrons is flat, and extends from -D i to D i :

k δ( -k ) = θ (D i -| |) 2D i , i = A1,B1, (11)
and we take D A1 = D B1 = 7eV as estimated from the density of states at the Fermi energy of s electrons for Cu atoms on the clean surface.

In principle, the 4s orbital of the Co atom should also be taken into account:

H s = s σ s † σ s σ + kσ V kA1s (s † σ c kA1σ + c † kA1σ s σ ) + t s pA1σ (c † pA1σ s σ + s † σ c pA1σ ) -S s • i J si S i , (12)
with on-site energy s , destruction operator s σ for spin σ , hybridization matrix elements V kA1s with A1 conduction states, coupling t s to A1 states of the tip, and exchange coupling J si with d i states. However, due do its large hybridization ( s ∼ 2eV), the 4s orbital can be taken as part of the A1 conduction band, thus effectively enhancing t 1A1 in Eq. ( 6) for the d z 2 orbital.

In practice, since we are mainly interested in the Kondo temperature, instead of solving Eqs. ( 5) and ( 6) together by NRG, we always solve Eq. ( 5) alone but replace the hybridization s i due to the surface with the total hybridization i = s i + t i due to the surface plus tip. As explained in Sec. II, the line shape is approximated during the DFT step without resorting to the model Hamiltonian in Eqs. ( 5) and (6).

The parameters in Eq. ( 5) are then fixed by trying to reproduce as closely as possible the GGA results within the HF approximation of the AIM in the wide-band limit [START_REF] Lucignano | [END_REF], which gives

↑ i = i + U i n ↓ i + j U ij n j - j J ij m j /4, ( 13
) ↓ i = i + U i n ↑ i + j U ij n j + j J ij m j /4, ( 14
)
where we sum over all j = i atomic orbitals (n σ i , with σ = ↑,↓, is the occupation of orbital i in the spin channel σ ;

n i = n ↑ i + n ↓ i ; m i = n ↑ i -n ↓ i ).
J = 2( ↓ xy - ↑ xy ) m . ( 15
)
Hubbard repulsions U i are taken from the splitting of magnetic orbitals, once J is known:

↓ i - ↑ i = U i m i + J 2 (m -m i ), i = d z 2 ,d x 2 -y 2 . ( 16
)
The interorbital Hubbard repulsion U 12 is approximated by [START_REF] Baruselli | [END_REF]]

U 12 = U ave -5 4 J, ( 17
)
where U ave is the average of the Hubbard repulsion over the d z 2 and d x 2 -y 2 orbitals. Finally, the on-site energies i are fixed from the orbital energies σ i together with the knowledge of U i and U 12 :

i = ↑ i + ↓ i 2 - U i 2 n i -U 12 n j , i,j = d z 2 ,d x 2 -y 2 ,j = i. ( 18
) Here, the numerical values of σ i , n σ i , n i , and m i are taken from DFT. The results of this procedure are reported in Table I.

C. Kondo temperature

In order to estimate the Kondo temperature, Eq. ( 5) for the d z 2 and d x 2 -y 2 orbitals is solved by a two-channel NRG calculation. Each orbital, with on-site energy i and Hubbard repulsion U i , is coupled to its own Wilson chain through the full broadening i = s i + t i , which takes into account the interaction with both the surface and tip; the two channels are coupled by Hund exchange coupling J and interorbital Hubbard repulsion U 12 . The Kondo temperature is obtained by computing the spectral function for both impurity levels, and taking the half width of the zero energy resonance.

GGA predicts the down-spin orbitals to lie just above the Fermi energy (see Fig. 2). When translated into an AIM, this means that ↓ i ∼ i + U i 0, which leads to high particle-hole asymmetry. As a consequence, NRG predicts the magnetic orbitals to be almost in the mixed-valence regime, which explains why the values of T K in Table I are higher than the experimental ones. However, we believe this high particle-hole asymmetry is a spurious effect, which could be amended by resorting to some more sophisticated method, such as GGA + U [41]. For example, we found that a small value of U ∼ 0.5 eV in the GGA + U approach is enough to reproduce the experimental T K in the tunneling regime.

Our value of z 2 0.18 eV in the tunneling regime is comparable, but somewhat lower than other values found in the literature: 0.20 eV in Ref. [28] and 0.24 eV in Ref. [9]. Even though our method of computing the hybridization from the broadening of scattering states is in principle more accurate than a simple estimation from the density of states after a self-consistent calculation, since it involves the interaction of impurity levels with a continuum of states, hence requiring no artificial broadening, the approximation of using a single k xy point can easily lead to a ∼10% error, as a consequence of the interaction among periodic replicas of the impurity in the xy plane, which could be alternatively reduced by the use of larger supercells. In addition, the way the nuclear relaxation is performed-either taking into account magnetic effects or ignoring them-is found to be another source of uncertainty. For example, performing a spin-unpolarized relaxation increases z 2 from 0.183 to 0.208 eV with respect to the standard spin-polarized calculations we use in this paper. This differs with earlier work [39], where magnetic and nonmagnetic calculations were found to yield similar relaxed atomic coordinates due to a cancellation between adatomsubstrate and adatom-tip interactions; see also Ref. [11]. For comparison, we find the relaxed Co-surface distance without the tip to be 1.59 Å when magnetism is taken into account, and 1.49 Å otherwise; in Ref. [39] this value was found to be 1.51 Å, regardless of magnetism.

As such, a precise evaluation of the Kondo temperature remains a challenge. However, we stress that we obtain the correct growth of T K,z 2 in passing from the tunneling regime to the contact regime, as reported in Table I and Fig. 5. In addition, in the intermediate regime, i.e., for g ∼ 0.3, T K,z 2 decreases slightly, which is probably true experimentally [see Fig. 3(a) of Ref. [10]]. This is a consequence of the fact that the Co adatom is attracted by the tip, and therefore pulled away from the surface (see Table I, column d Co-sur ; this agrees well with Ref. [39]). This causes a decrease in the hybridization z 2 , which is compensated only in the contact regime, when the tip is close enough to the Co atom to be considered as an additional nearest-neighbor atom. At this point, z 2 and T K,z 2 are greatly enhanced. However, this is only true for the d z 2 orbital, which is probably the one producing the experimental ZBA. In the case of d x 2 -y 2 , the tip greatly decreases T K,x 2 -y 2 when approaching the adatom, because the decreased hybridization of the orbital with the surface due to the increased surface-adatom distance is not compensated by an additional hybridization with the tip, due to symmetry mismatch; should symmetry be broken, things would not be different, since the d x 2 -y 2 orbital lies flat on the surface, thus hardly coupling with the STM tip no matter how this approaches the adatom. This is true until the tip really "pushes" the atom closer to the surface, eventually enhancing T K,x 2 -y 2 too. We stress that, in contrast to Refs. [9,10], we attribute the change of the Kondo temperature mainly to changes in the hybridization, rather than in the on-site energy and Hubbard repulsion.

D. Line shape

As remarked, in our symmetry-preserving geometry, the total line shape is the sum of the A1 and B1 conductances in Eq. (7). However, after solving Eq. ( 5) with the parameters shown in Table I, we find that T K,z 2 T K,x 2 -y 2 ; moreover, g A1 g B1 because the d x 2 -y 2 orbital is flat and only couples to the second layer of the tip for symmetry reasons. This suggests that most of the experimental Kondo signal is due to the d z 2 orbital and, in what follows, we shall assume g tot = g A1 and ignore g B1 , together with all other symmetry channels, which do not carry ZBA's and contribute very little to the conductance. In any case, we find that the Fano parameter associated with the orbital d x 2 -y 2 should always be much larger than 1, and this would show up as an additional weak anomaly superimposed to the standard one.

Applying the phase-shift analysis described in Ref. [START_REF] Lucignano | [END_REF] turned out to be too cumbersome for this system, due to its intrinsic three-dimensional character, so instead we estimated q from the shape of the energy dependent DFT transmission coefficient, which shows interference at energies around ↓ d ∼ 0.5 eV. If we assume that the hopping parameters in Eq. ( 6) are weakly energy dependent, the DFT line shape can be a good approximation to the ZBA, which involves interference at the Fermi energy. Unfortunately, when we do so, the agreement with experiment is not always good. In the tunneling regime, we obtain q ∼ 1, which nicely matches experiments [8]. However, when going into the contact regime, we find a decrease of q. In this regime the conductance G is close to the unitary value G 0 , thus the interference between the s and d z 2 orbitals can only be destructive, leading to a dip in the conductance, while in experiments the opposite is observed (q increases slightly [9] or strongly [10]).

A possible reason for this disagreement in the contact regime is that the junction is formed in a different way than we have modeled it. Other sources of error in our calculations are the inclusion of nonequilibrium effects and our assumption of energy independent parameters; see also the conclusions in Sec. V.

IV. Co/Cu(111)

In this section we report our results for Co on the Cu(111) surface, which has T K = 54 K and q = 0.18 in the tunneling regime [8]. These values remain almost unchanged in passing to the contact regime [11].

A. DFT results

Once again, we find that Co adsorbs in the hollow position, this time with three nearest-neighbor Cu atoms in a configuration with C 3v symmetry. The d orbitals split into a singlet with symmetry A 1 (d z 2 ) and two doublets with symmetry E (d α1,2 , d β1,2 ); the s orbital has A1 symmetry. The electronic configuration is 3d 8 4s 1 , and the total magnetic moment is close to 2μ B , just like on the (100) surface. We again model the tip to preserve the symmetry (C 3v ).

If we take Cartesian axes as in Fig. 3, the doublets can be written in the following way:

d α1 = cos θd xz + sin θd xy , ( 19
)
d α2 = cos θd yz + sin θd x 2 -y 2 , ( 20
)
d β1 = -sin θd xz + cos θd xy , ( 21
)
d β2 = -sin θd yz + cos θd x 2 -y 2 , ( 22
)
where the d α1,β1 orbitals are odd with respect to the symmetry operator P x : x → -x, while the d α2,β2 orbitals are even. From DFT calculations, it turns out that the d α1,2 doublet is magnetic, while d β1,2 is fully occupied; moreover we find θ = 0.70 rad.

In Table II we show some structural and electronic data for different tip-surface separations. In Fig. 4 we show density of states and transmission eigenvalues for the shortest distance, d tip-sur = 4.33 Å. At K, the degeneracy between d α1 and d α2 is weakly broken, so two different peaks appear in the density of states and in transmission eigenvalues at ∼-0.5 eV (d β orbitals) and ∼0.5 eV (d α orbitals) for down-spin electrons in the E channels.

B. Anderson model

The full atomic Hamiltonian is

H = i=z 2 ,α1,α2,β1,β2 ( i n i + U i n i↑ n i↓ ) - j<i J ij S i • S j + j<i U ij n i n j , (23)
with on-site energies α1 = α2 ≡ α and β1 = β2 ≡ β , Hubbard repulsion U α1 = U α2 ≡ U α and U β1 = U β2 ≡ U β , interorbital Hubbard repulsion U ij , and Hund exchange coupling J ij . After dropping fully occupied orbitals, which means keeping only d α1 and d α2 , we introduce metallic states, to get i and t ↓ i (eigenvalues close to 1 correspond to the A1 channel while two lower eigenvalues are from the E channels), and the shot-noise Fano factor. Energies are with respect to the Fermi energy.

H ad-sur = kσ i = α1,α2 k c † kiσ c kiσ + kσ i = α1,α2 V kα (c † kiσ d iσ + d † iσ c kiσ ) + i=α1,α2 ( α n i + U α n ↑ i n ↓ i ) + U 12 n α1 n α2 -J S α1 • S α2 . ( 24
In addition to the d orbitals, the s orbital is half filled and highly hybridized, exactly as on the (100) surface, with the same Hamiltonian [Eq. ( 12)], but it is irrelevant when dealing with Kondo physics: it only contributes to the conductance in the A1 channel which shows no ZBA, and does not interfere with the d α orbitals, which have different symmetry E. TABLE II. Parameters for Co/Cu(111) at different tip-surface separations d tip-sur : the distance d tip-Co between the tip and the Co adatom, d Co-sur between the Co adatom and the surface, d Co-NN between the Co adatom and its three nearest neighbors, the total magnetization m of the unit cell in units of μ B , the on-site energy α , the Hubbard repulsion U α , the hybridization α , the interorbital Hubbard repulsion U 12 , the Hund exchange J , the calculated Kondo temperature T K in K, the Fano parameter q, the DFT conductance g = g ↑ + g ↓ , and the Fano factor F at the Fermi energy; distances are in Å and energies are in eV.

d tip-sur d tip-Co d Co-sur d Co-NN m α U α α U 12 J T K q g F 7.
H tip = pσ i = α1,α2 pi c † piσ c piσ + pkσ i = α1,α2 t 1i (c † piσ c ki + c † kiσ c piσ ) + pσ i = α1,α2 t 2i (c † piσ d iσ + d † iσ c piσ ), (25)
is used and the hybridization functions

s α ( ) = π k δ( -k )V 2 kα → πρ s V 2 α , ( 26
) t α ( ) = π k δ( -k )t 2 kα → πρ t t 2 2α , ( 27
) α ( ) = s α ( ) + t α ( ) (28)
are introduced and approximated as energy independent quantities s α , t α , and α . The parameters are fixed in the same way as for the Co/Cu(100) case; we get D α = 5 eV from the density of states of surface s orbitals; J is fixed from the splitting of β orbitals, which have m β = 0.06. Numerical values are reported in Table II. The value α ∼ 0.16 eV is slightly below the value 0.18 eV of Ref. [23].

C. Kondo temperature

The Kondo temperature is obtained by solving Eq. (24) with NRG for the model parameters in Table II. We used the total hybridization α , including contributions from both the surface and tip.

Since the d α1 and d α2 orbitals are degenerate, there is a single Kondo temperature as reported in Table II. Spin-orbit effects will lift this degeneracy, leading to two different Kondo temperatures, but the effect should be small. As for the previous case, we overestimate the Kondo temperature due to the excessive particle-hole asymmetry which comes from GGA, but in this case the disagreement is not too bad.

When the tip is brought down to the surface, we find that the Kondo temperature first decreases slightly (even though the difference is below the accuracy of our method), because the adatom is pulled farther away from the surface, then weakly increases. The s orbital of the tip is of different symmetry than the α orbitals, so t α = 0 to a first approximation. The effect of the tip-induced relaxation of the Co adatom and its neighbors is found to be negligible, in contrast to the (100) case, so s α is basically unaffected by the tip. This is in good agreement with a similar analysis in Ref. [11] and with experimental results that show a constant T K as a function of the tip position [11]. At a distance d tip-sur = 4.33 Å, the Co atom is pushed towards the surface, and at this point the Kondo temperature starts to increase, but this regime is probably not reached in experiments.

In Fig. 5 we show the evolution of the conductance g and Kondo temperature T K,i as a function of the tip-surface distance for both Co/Cu(100) and Co/Cu(111). The reported values of T K,i are to be seen mostly as upper estimates, since the use of GGA + U , as mentioned, would decrease the Kondo temperature. In fact, our values overestimate the experimental Kondo temperature. It must be stressed that most of the uncertainties are systematic, so they affect all the data in the same direction.

D. Line shape

When the tip is above the adatom, symmetry is preserved, and the total conductance is

g tot = i=A1,α,β g i 2g α + G A1 . ( 29
)
In the A1 channel, there is no ZBA, because only the s orbital is involved, the d z 2 orbital being completely filled. In the E channel, there is a ZBA due to d α orbitals, but the signal should be small because they do not couple to the s orbital of the tip. DFT predicts that these channels give a peak in the conductance (q 1, see Fig. 4), the coupling of the tip to the orbitals being much higher than to Cu states, according to the two-path model, but still much lower than the coupling to the Co s orbital. However, when compared to experiment, where a dip (q ∼ 0) in conductance is observed, with a strong signal, this is wrong. Moreover, on the basis of the DFT results, one would expect to see considerable changes in the line shape when moving the tip laterally in the xy plane, due to symmetry breaking, but this is not seen either. Finally, when the tip approaches the surface, the line shape is unaffected, remaining a minimum, again in contradiction with what one would expect from the DFT results.

This probably means that symmetry is unimportant, most likely because the tip breaks it. This implies that one can observe E orbitals even when the tip is above the adatom, and when moving the tip laterally there is no symmetry breaking and no significant change in the signal. However, this is not sufficient to explain the experimental results because one would still expect the ZBA to be a peak. Another interpretation might be that since the magnetic orbitals have E symmetry, while tunneling only happens through metal states of A1 symmetry, these acquire a dip in their density of states, leading to q ∼ 0, as suggested in Ref. [37]; however, metal states of A1 symmetry should carry no, or at most very weak, Kondo signal.

Instead, it is likely that surface states are responsible for the dip. According to DFT, the tip "sees" the magnetic orbital and not the metallic states, while it should be the other way around. This is simply because the Co adatom is closer to the tip than to the Cu surface atoms [on the (100) surface things are different, because the Co s orbital, which is effectively part of the conduction band, can interfere with the magnetic orbital]. However, surface states, if taken into account, might prove to be more prominent than the d orbitals, leading to a dip in conductance. This interpretation is also suggested by the "quantum mirage" [42] experiment, which shows how surface states can carry Kondo information even far from the impurity, with the same line shape as when the tip is above the adatom. Our DFT slab calculation can in principle include surface states, but to describe them correctly we would need a much larger supercell. Å for Ref. [11] for best fit with theory. Approximate error bars are shown for the theoretical Kondo temperatures: they take into account all the incertitudes of our method when fixing the parameters of the Anderson Hamiltonian Eqs. ( 5) and ( 24) from the DFT calculations; once the parameters are known, the Kondo temperature can be extracted with negligible uncertainty from the NRG run. It must be stressed that most of the incertitudes are systematic, so they affect all the data in the same direction, which turns out, when compared to experiments, to be an overestimation of the Kondo temperatures.

The role of surface states has been discussed in the literature. The total hybridization from the surface s = surf + bulk is the sum of the hybridization from surface states and bulk states. While it is generally agreed [8,23,27,[43][44][45], the only exception being Ref. [46], based on calculations and experimental hints (for example, the lack of appreciable changes in the Kondo temperature at step edges and defects [47], where surface states are deeply affected), that bulk is much larger than surf , by up to a factor of 100 [23], not much is known about the relative magnitude of tip-surface and tip-d orbitals coupling, which is what controls the line shape. It is only known that for a clean surface about two-thirds of the current flows into surface states [48]. Reference [49] argues that surface states can give an important contribution to the conductance even in the presence of adsorbates.

In any case, it must be stressed that the usual assumption that the magnetic orbital is of d z 2 character is found to be wrong in this case. This makes the usual Tersoff-Hamann approach fail, because due to symmetry mismatch there can be no conductance from the s state of the tip into the magnetic orbital, or to the linear combination of conduction states to which the magnetic orbital is coupled, which leads to the paradox that when the tip is directly above the adatom it should give almost no signal.

V. CONCLUSIONS

We have demonstrated the application of our DFT + NRG method [START_REF] Lucignano | [END_REF] to Co adatoms on Cu surfaces probed by an STM tip, trying to describe the line shape and the Kondo temperature both in the tunneling and contact regimes.

The calculations show the severe difficulty in predicting the details of experimental Kondo anomalies, especially as far as the line shape is concerned, while the trend in the Kondo temperature, if not its absolute value, is reliable. Several different issues might cause these discrepancies.

First of all, our description of the tip is surely oversimplified. Yet, experimental measurements are only weakly dependent on the choice of the tip, at least in the tunneling regime, so in principle this should not be a large source of error.

We believe that on the Cu(100) surface the essential remaining point is to correct the estimation of the parameters which control the line shape: the asymmetric ZBA in the tunneling regime can be understood in terms of the interference between the s and d z 2 orbitals of the Co adatom. Issues may arise in the contact regime due to the geometrical details of the contact, which may differ from our model. The momentum dependence of the hopping parameters such as V k , t 1ipk , and t 2ip , which we have ignored, might also play a role, as well as nonequilibrium effects. Also, in the contact regime, which is not far from a bulk impurity situation, all d orbitals start to become magnetized, thus actively entering conduction processes, and making our two-orbital model insufficient.

On the Cu(111) surface, in contrast, more work should be done to include surface states. This might reverse the sign of the ZBA (giving a dip instead of a peak). Also, one should seek to understand if the E symmetry of the magnetic orbital, instead of A1 as usually assumed, can affect the line shape, as suggested, e.g., in Ref. [37].

On both surfaces our estimate of the Kondo temperature would be improved by correcting the excessive particle-hole asymmetry brought about by plain GGA, for example, by using GGA + U [41]. Moreover, our many-body model may be improved. For example, one could take into account spinorbit effects, correlated hopping, double hopping, and other two-body interactions, or keep all d orbitals, considering that some of them are almost but not completely filled, especially in the contact regime. In addition, the energy dependence of the hybridization functions i ( ) is likely to have some impact. Finally, we note that the coexistence of many magnetic solutions (not only of lowest energy presented in this paper) realized at different Co adatoms (on both Cu surfaces) could be responsible for statistical spread in Kondo temperatures and line shapes.

We emphasize that, according to our GGA results, Co has spin S = 1 on both Cu surfaces, showing that the usual assumption of S = 1/2 is probably wrong. However, this mistake might not have a big impact on the final result if one of the two magnetic orbitals has a much lower Kondo temperature than the other one, as on the (100) surface, or if the two magnetic orbitals are degenerate, as on the (111) surface. Of course, a detailed quantitative approach cannot overlook this fact.

In any case, despite several works claiming that the Kondo physics of adatoms is fully understood, we believe that further effort is needed to completely understand, or at least describe satisfactorily from first principles, the Kondo behavior of Co adatoms on Cu surfaces and, more generally, of magnetic adatoms on metallic surfaces. Future improvements should include effects of spin-orbit and anisotropy, and the direct treatment of electrical dissipation beyond the Landauer linear response approach used here.

Note added in proof. Calculations of Kondo conductance anomalies for Co/Cu(100) using a different ab initio + manybody method have been published since the submission of this paper [50] with results in agreement with ours; in particular the Co d z 2 orbital was found to host the observed Kondo resonance, as here.

FIG. 1 .

 1 FIG. 1. (Color online) Unit cell of the scattering region for Co/Cu(100) for d tip-sur = 7.60 Å (a) in the yz plane and (b) in the xy plane. (c) Brillouin zone for the supercell in the xy plane.

FIG. 3 .

 3 FIG. 3. (Color online) Unit cell of the scattering region for Co/Cu(111) for d tip-sur = 7.81 Å (a) in the yz plane and (b) in the xy plane. (c) Brillouin zone for the supercell in the xy plane.

FIG. 4 .

 4 FIG. 4. (Color online) Co/Cu(111): (a) Spin-polarized density of states at Co 3d and 4s atomic orbitals constructed from scattering states at K = π L ( 2 3 ,0) (positive DOS means spin up; negative DOS means spin down); (b) t ↑ and t ↓ components of the DFT transmission, transmission eigenvalues t ↑ i and t ↓ i (eigenvalues close to 1 correspond to the A1 channel while two lower eigenvalues are from the E channels), and the shot-noise Fano factor. Energies are with respect to the Fermi energy.

FIG. 5 .

 5 FIG. 5. (Color online) (a) DFT conductance g and (b) Kondo temperature T K ,z 2 on the Cu(100) surface compared with experimental data from Refs. [10] (Choi et al.) and [9] (Néel et al.). (c) DFT conductance g and (d) Kondo temperature on the Cu(111) surface compared with experimental data from Ref. [11] (Vitali et al.). All data are shown as a function of the distance between the tip and the surface d tip-sur : since only relative distances are known experimentally, we have rigidly shifted experimental data by 5 Å for Ref. [10], 9 Å for Ref. [9], and 5.5Å for Ref.[11] for best fit with theory. Approximate error bars are shown for the theoretical Kondo temperatures: they take into account all the incertitudes of our method when fixing the parameters of the Anderson Hamiltonian Eqs. (5) and (24) from the DFT calculations; once the parameters are known, the Kondo temperature can be extracted with negligible uncertainty from the NRG run. It must be stressed that most of the incertitudes are systematic, so they affect all the data in the same direction, which turns out, when compared to experiments, to be an overestimation of the Kondo temperatures.

  Here, z 2 and U z 2 have large errors (see text).

	60	6.01 1.59	2.43 2.08 -4.83 3.06 0.183 -4.89 3.24	0.147 1.56 1.28 340	50	1.20 0.01 0.99
	5.58	3.91 1.67	2.47 2.10 -4.67 2.98 0.165 -4.80 3.26	0.130 1.50 1.29 290	10	1.19 0.31 0.68
	5.17	3.41 1.76	2.51 2.10 -4.54 2.91 0.170 -4.67 3.23	0.116 1.44 1.30 410	3	0.76 0.61 0.53
	4.73	2.84 1.89	2.56 2.08 -4.47 2.86 0.195 -4.60 3.23	0.101 1.40 1.32 600	0.5	0.09 1.06 0.42
	4.48	2.61 1.87	2.54 2.05 -4.58 2.88 0.211 -4.65 3.22	0.104 1.40 1.32 1000	1	0.03 1.27 0.38
	4.12	2.42 1.70	2.48 1.96 -5.47 a 3.36 a 0.232 -5.23 3.22	0.122 1.63 1.33 1100	25	0.01 1.36 0.34

a

  The linewidths i are taken from the down-spin density of states by fitting each impurity orbital with a Lorentzian; this actually gives i ( It is inferred from the energy splitting of the d xy orbital (which has a very low magnetization m xy ∼ 0.04μ B ), induced by the total magnetization m via exchange interactions:

↓ i ), but we assume i ( ↓ i ) i (0)

. J is assumed to be constant in the d shell.

ACKNOWLEDGMENTS

Work in Trieste was supported by PRIN-COFIN Contract

No. 2010LLKJBX 512004, and by the European Union ERC Advanced Grant No. 320796 -MODPHYSFRICT. We acknowledge useful discussions with R. Zitko.