Many-particle effects in adsorbed magnetic atoms with easy-axis anisotropy: the case of Fe on CuN/Cu(100) surface

We study the effects of the exchange interaction between an adsorbed magnetic atom with easy-axis magnetic anisotropy and the conduction-band electrons from the substrate. We model the system using an anisotropic Kondo model and we compute the impurity spectral function which is related to the differential conductance (dI/dV) spectra measured using a scanning tunneling microscope. To make contact with the known experimental results for iron atoms on the CuN/Cu(100) surface [Hirjibehedin et al., Science {\bf 317}, 1199 (2007)], we calculated the spectral functions in the presence of an external magnetic field of varying strength applied along all three spatial directions. It is possible to establish an upper bound on the coupling constant J: in the range of the magnetic fields for which the experimental results are currently known (up to 7T), the low-energy features in the calculated spectra agree well with the measured dI/dV spectra if the exchange coupling constant J is at most half as large as that for cobalt atoms on the same surface. We show that for even higher magnetic field (between 8 and 9T) applied along the ``hollow direction'', the impurity energy states cross, giving rise to a Kondo effect which takes the form of a zero-bias resonance. The paper introduces an approach for calculating the expectation values of global spin operators and all components of the impurity magnetic susceptibility tensor in numerical renormalization group (NRG) calculations with no spin symmetry. An appendix contains a density-functional-theory (DFT) study of the Cu and Fe adsorbates on CuN/Cu(100) surface: we compare magnetic moments, as well as orbital energies, occupancies, centers, and spreads by calculating the maximally localized Wannier orbitals of the adsorbates.


Introduction
Recent advances in the experimental techniques, such as low-temperature scanning tunneling microscopy (STM) and the X-ray magnetic circular dichroism (XMCD), made it possible to study magnetic properties of single adsorbed atoms on various surfaces [1,2,3,4,5,6,7]. Particularly interesting are magnetic impurities on noble metal surfaces where many-particle effects such as the Kondo effect play an important role [8,9]. Due to the reduced symmetry in the surface region, such magnetic adsorbates are strongly anisotropic [1,10]. A simple parametrization of the leading magnetic anisotropy terms takes the form of where the S x , S y , and S z are the cartesian components of the quantum-mechanical spin-S operator, while the parameters D and E are known as the longitudinal and transverse magnetic anisotropy, respectively. If D < 0 (the "easy-axis" case), the spin tends to maximize the absolute value of the z component, thus the S z = ±S states dominate in the ground state doublet, which may be split for integer S even in the absence of an external magnetic field due to the transverse magnetic anisotropy E. If D > 0 (the "easy-plane" or "hard-axis" case), the ground-state is a doublet consisting mostly of the S z = ±1/2 states for half-integer spin, or a singlet consisting mostly of the S z = 0 state for integer spin. Similar models are also used to describe molecular magnets [11,12,13,14,15]. When adsorbed on metallic surfaces, the coupling of the spin to the conduction-band electrons can introduce new phenomena, the Kondo effect being the most interesting due to its unique signature expected in the tunneling spectra. For half-integer spins with easy-plane anisotropy, the level degeneracy in the ground state actually permits the occurrence of the Kondo effect at zero magnetic field, as indeed observed experimentally in the system of cobalt atoms adsorbed on the CuN ultra-thin layers deposited on the Cu(100) substrate surface [16]. For magnetic atoms with easy-axis anisotropy, however, the Kondo effect has not yet been observed. Its absence is either due to the zero-field level splitting caused by the transverse magnetic anisotropy E, which in the case of integer S plays the role of an effective magnetic field which suppresses the Kondo effect, or (for half-integer S or for E ≪ |D|) due to the fact that the spin-flip scattering events of the conduction-band electrons can only change the impurity spin component by 1.
Thus the two impurity ground-state levels only couple via higher-order processes and the Kondo screening is strongly suppressed, i.e., the Kondo scale is extremely small [12]. Nevertheless, interesting many-particle effects in fact do occur in the case of easy-axis anisotropy, thus we study this class of the problems in the following parts of the paper, focusing in particular on the case of spin 2 as appropriate for iron atoms adsorbed on the CuN/Cu(100) substrate [10].
The absence of the Kondo effect at zero magnetic field does not imply, that the exchange coupling of an easy-axis magnetic impurity to the substrate is of no consequence. In this work we show that the value of the exchange coupling constant J, defined through the Kondo coupling Hamiltonian where s is the spin-density of the conduction-band electrons at the position of the impurity, affects the impurity spectral function in a measurable way. The most drastic effect is the occurrence of a Kondo effect induced by a magnetic field, similar to the singlet-triplet Kondo effect in quantum dot systems [17,18,19,20,21,19,22]: for some high enough external magnetic field applied along the x-axis, the impurity level crossing leads to a situation where the conditions for the occurrence of a Kondo effect are met. We will show that the presently available experimental results (measured up to 7 T, Ref. [10]) just stop short of the region where this Kondo resonance should be observed (between 8 T and 9 T). For this reason, the presently known experimental results can only be used to establish an upper bound on the value of J, but by extending the magnetic field range into the 9 T range, the Kondo resonance could be probed and the exchange coupling J might be determined even quantitatively. The paper is structured as follows. In Sec. 2 we provide the full definition of the model under study and describe the numerical technique which is used to compute the impurity spectral function. We also comment on the relation between the impurity spectral function and the experimentally measured dI/dV spectrum. In Sec. 3 we discuss the impurity level structure and the predicted dI/dV spectra in the limit of a decoupled impurity, J → 0. Such spectra do not encompass any many-particle effects; however, they account accurately for the high-energy spin excitations which are observed experimentally. In Sec. 4 the results of the numerical calculations for the full manyparticle problem are presented and discussed in relation with the experimental results. Appendix A presents a comparative density-functional-theory (DFT) study of Co and Fe adatoms on the CuN/Cu(100) surface, while Appendix B contains the description of a method for computing the thermodynamic impurity spin susceptibility in problems having no symmetries in the spin space with the numerical renormalization group (NRG) method.

Model and method
We model the system using the Hamiltonian The anisotropy term H aniso is given by Eq. (1), and the coupling term H K by Eq. (2). The operator c kσ corresponds to a conduction-band electron with momentum k, spin σ ∈ {↑, ↓}, and energy ǫ k . In terms of these operators, the spin density of the conductionband electrons at the position of the impurity (assumed to be at the origin of the coordinate system) is given by N is the number of the conduction-band states, and σ αβ is the vector of Pauli matrices. Finally, g is the gyromagnetic factor (assumed to be isotropic), and µ B ≈ 0.0579meV/T is the Bohr magneton. For an iron adatom on the CuN/Cu(100) surface, the following parameters have been established: spin S = 2, g = 2.11 ± 0.05, D = −1.55 ± 0.01 meV, E = 0.31 ± 0.01 meV (Ref. [10]). For a complete description one should, furthermore, know the density of states (DOS) of the conduction band, ρ(ω) = (1/N) δ(ω − ǫ k ), and the value of the Kondo coupling J. In the interesting energy range (on the order of 10 meV around the Fermi level) the energy dependence of the DOS may be neglected and we can use a constant value ρ = ρ(ω = 0). The low-energy behavior of the problem then depends solely on the dimensionless combination ρJ. We note that the model does not include any orbital moment degrees of freedom. This approximation is based on the fact that the orbital moment of an impurity adsorbed on the surface is typically strongly quenched due to the surface field effects (for the particular case of Fe and Co adsorbates on the surface of CuN/Cu(100), the total suppression of orbital degeneracy follows from the DFT results tabulated in Appendix A). We solve the problem using the numerical renormalization group [23,24,25] (for the details on the application of the NRG to the problems of this class see also Refs. [26,27]). In this work, the focus of our interest will be the impurity spectral function defined as [28] A σ (ω) = − 1 π ImT σ (ω + iδ), where T σ (ω) is the T-matrix for the scattering of electrons of spin σ on the spin-S impurity, given by with the operators Here S ± = S x ± iS y . The notation A; B ω denotes the Laplace transform of the correlator A; B t = −iθ(t) {A(t), B(0)} + between two (fermionic) operators A and B, i.e., A; B ω = ∞ 0 e iωt A; B t dt with Im ω > 0; {A, B} + = AB + BA is the anti-commutator.
The T-matrix itself is defined through where G kkσ (ω) = c kσ ; c † kσ is the conduction-band electron Green's function and G 0 kkσ its unperturbed counterpart. In the STM experiments, the spin components are not resolved, thus we will mostly plot the spin-averaged spectral functions, A(ω) = [A ↑ (ω) + A ↓ (ω)]/2.

Decoupled impurity limit
The spin excitation spectrum of a magnetic atom adsorbed on a thin insulating layer may for the most part be accounted for by fully neglecting its exchange coupling to the substrate electrons. According to this picture, the tunneling differential conductance will increase stepwise when the bias voltage is increased beyond the values (E i − E gs )/e, where E i are the energies of the excited states and E gs is the ground-state energy [2,3,10]. At each of these voltages, an additional inelastic scattering (spin-flip) channel opens, typically resulting in an increased conductance [29,30,31,32,33,34,35,36,37]. In fact, the experiments have shown that the conductance step heights are described to a good approximation by [2] | i|S i.e., by the transition matrix elements for the spin operator. This empirical observation has recently received theoretical support [38,39]. The expression suggests that the impurity is well thermalized with the substrate: the impurity relaxes to the ground-state on a time-scale which is short compared with the mean time between the tunneling events. Using this simple approach, one can compute dI/dV spectra which compare very well with the experimental results (as long as the many-particle effects are not important) if thermal broadening is taken into account. Here we show the results obtained by this procedure without adding any artificial thermal smearing, see Fig. 1a. The plots should be compared with Fig. 2A,B,F and Fig. S1,A,C in Ref. [10]. Presenting the results in the zero-temperature limit uncovers more details, in particular the crossing of the the levels in the case of a field applied along the x-axis (this corresponds to the "hollow axis" in Ref. [10]). As is well known, the experimental results are in very good agreement with the simulated dI/dV spectra [10]. One comment, however, is in order. In experiments, the step height of the first step in the dI/dV spectra is observed to decrease as the field strength along the x-axis is increased. This can be explained, on one hand, by the thermal broadening of the spectra which leads to reduced heights of spectral features of the width smaller than the width of the broadening kernel. Another explanation is based on a possibility of a slight misalignment of the magnetic field from the direction of the actual x-axis of the magnetic atom. A simple calculation (Fig. 1b) shows that even a small misalignment of 1 • (which is quite likely to occur) can suppress the step height of the first excitation (at T = 0) without significantly modifying those at the other magnetic excitation energies.
We now comment on the important observation that the energy difference between the two lowest levels decreases as the field strength along the x-axis is increased, see Fig. 1c, top. A simple calculation neglecting any possible many-particle effects indicates that the levels become degenerate at For this field strength, a magnetic-field-induced Kondo effect can occur if the Kondo exchange coupling is antiferromagnetic, J > 0. We note that the aforementioned misalignment of the field will also tend to suppress the Kondo effect at B x ≈ 8.8 T, since a component of the field along the z-axis will mix the levels, thereby leading to an avoided crossing of the levels (Fig. 1c, bottom). In experiments seeking to detect the field-induced Kondo effect, it will be thus essential to not only fine-tune the magnitude of the magnetic field, but also its exact direction.

Impurity spectral functions
In Fig. 2 we plot the impurity spectral functions A(ω) computed using the numerical renormalization group (NRG). The calculation is performed for zero temperature, thus the widths of spectral features are due to intrinsic effects (the NRG spectral function overbroadening artifacts are small in these calculations; for discussion see Ref. [40]). Compared with the simulated dI/dV curves, all spectral features except for those near the Fermi level are almost completely smeared out. This is at odds with the experimental results, where the steps are observed to be sharp apart from the thermal broadening. While the calculated spectral function does encompass some inelastic scattering effects (see Fig. 3 in Ref. [27]), the tunneling current includes additional inelastic scattering contributions which are beyond the scope of our calculations that are constrained to the case of the magnetic impurity being in thermal equilibrium with the substrate electron reservoir at all times and do not include any additional interactions between the tunneling electron and the impurity local moment. Nevertheless, the results obtained from a NRG calculation provide a reliable description of the behavior of the system in the immediate vicinity of the Fermi level, which is probed by tunneling electrons of sufficiently low-energy so that the system is not driven significantly out of the equilibrium.  Figure 2. Impurity spectral functions A(ω) computed using the NRG. The horizontal energy axis has been transformed to show the equivalent bias voltages, while the vertical axis corresponds to the dI/dV differential conductance in the simple approximation where the LDOS is taken to be proportional to the spectral function.
The agreement between the low-energy features (interval [−1 mV : +1 mV]) in the calculated and measured spectra [10] in the 0 to 7 T range is very good if thermal broadening is taken into account. The calculation was performed for ρJ = 0.05. This value is half as large as the exchange coupling for Co impurities on the same surface, ρJ = 0.1, which has been established in Refs. [26,27]) by determining the value of ρJ that allows to reproduce the experimental differential conductance spectra for the Co/CuN/Cu(100) adsorbate system. Due to similarities between Co and Fe atoms the exchange coupling constants should not be too different, but very likely of the same order of the magnitude and of the same sign (Co and Fe are neighbours in the periodic system; see Appendix A for a more rigorous comparative study using the density functional theory). In Fig. 3a we compare the low-energy part of the spectral function for a fixed magnetic field of 7 T in the x-axis direction for a range of values of J. For large J, the dip in the spectral density turns into a (Kondo) resonance which becomes fully developed for ρJ ≥ 0.1. For ρJ in the interval 0.06 to 0.1, we find an intermediate structure: a dip with protrusions at the step edges. At high enough temperature, the protrusions will be washed out and a single antiresonance would be detected. Taking all these elements into account, the estimate of the upper bound ρJ = 0.05 ± 0.02 therefore seems reasonable.  The magnetic field at which the Kondo effect occurs for integer spin Fe is not given exactly by Eq. (12), but it rather depends on the value of J, since the different impurity levels are renormalized differently due to the anisotropy. In the J → 0 limit, the Kondo effect occurs at B x ≈ 8.8 T, while for finite J the magnetic field where the Kondo peak may be observed is reduced. In addition, the parameter J controls the width of the Kondo resonance (and, of course, the Kondo temperature itself), thus for smaller J the interval of the magnetic fields where the Kondo effect may be observed is accordingly narrower and the Kondo scale lower.
The magnetic-field-induced Kondo effect exhibits a number of rather peculiar features. We first note that it is very difficult to observe the fully developed Kondo resonance unless J is rather large. For the estimated exchange coupling of ρJ = 0.05, where the Kondo effect occurs for B * x = 8.535 T, the Kondo temperature is actually extremely low, on the order of 0.1 mK, see Fig. 4. Even at dilution refrigerator temperatures, it thus appears unlikely that a fully developed Kondo resonance will be observed. This does not, however, preclude the experimental detection of the Kondo physics at play in this system, since the evolution of the tunneling spectrum as the magnetic field is tuned should exhibit a characteristic merging of the two "halfresonances".
The results for the thermodynamic properties also reveal that the zz component of the magnetic susceptibility multiplied by the temperature [i.e., the effective moment T χ zz (T )] exhibits a plateau at very high value of 3.5 (in suitable units), which is rather near the limiting value for a decoupled S z = ±2 doublet, which is 4. This plateau corresponds to a ln 2 plateau in the impurity entropy, i.e., to an effective two-state system. On the scale of T K , the moment is then screened. Note that the temperature dependence of T χ zz (T ) is different (Kondo-like, see the arrow in Fig. 4) for B x = B * x as compared with those for other values of the magnetic field, which merely follow the exp(−∆E/k B T ) rule, where ∆E is the difference between the two lowest impurity level energies.
It is equally interesting to note that the spectral functions change discontinuously as we cross the transition value B * x , see Fig. 5a; this corresponds to the system being in different ground states on either side of the transition. In the presence of a small magnetic field component along the z-axis, the level-crossing is replaced by a rapid cross-over. Any such rapid change of the spectral function over an extended energy interval would be an indicator of the proximity to the Kondo regime, even at higher temperatures.
The magnetic-field-induced Kondo effect is not the only observable consequence of the exchange coupling with the substrate. The value (as well as the sign) of J is actually reflected in the spectral function for any value of the magnetic field. In Fig. 6 we compare, for example, the impurity spectral function at zero field for a range of coupling constants. The exchange coupling strongly affects the features around the Fermi level, but also those at higher energies, although to a lesser degree.
The results which we have established for the S = 2 Kondo model also apply more generally for other integer spins: at some B * x there will be level crossing which can induce a Kondo resonance at the Fermi level, but no such crossing occurs for the magnetic field along the y and z axis. The same is true for half-integer spin models, but in this case there is, in addition, a two-fold level degeneracy also for zero magnetic field, thus for half-integer-spin magnetic impurities one could observe both the zero-field Kondo effect and the magnetic-field-induced Kondo effect.

Conclusion
We have studied thermodynamic and spectral properties of the Kondo model with an easy-axis anisotropy. By comparing the calculated spectral functions with the measured  We plot the impurity contribution to the magnetic susceptibility in various directions, χ imp,αα , where α ∈ {x, y, z} is a direction in space, and the impurity contribution to the entropy, S imp . The Kondo temperature is estimated from the data for S imp to be T K ∼ 10 −5 meV ∼ 0.1 mK. The arrow points to the curve with exhibits Kondo-like temperature dependence.
ones in the magnetic field range up to 7 T we have estimated the upper bound for the exchange coupling of a Fe adatom on CuN/Cu(100) surface with the bulk conductionband electrons to be ρJ = 0.05 ± 0.02. By extending the range of the magnetic fields in the experiments up to 9 T, it should be possible to obtain an improved upper bound on J or perhaps even observe a magnetic-field-induced Kondo resonance if J is high enough. We note here that, strictly speaking, the quantity J scales under the renormalizationgroup flow, thus its actual value depends on the choice of the bandwidth and the corresponding density of states ρ; a more rigorous treatment of the problem would take into account the full density-of-states dependence of the substrate and, furthermore, the momentum dependence of the Kondo coupling J(k, k ′ ), see for example Ref [41]. Nevertheless, a simple characterization in terms of the dimensionless combination ρJ is sufficient to compare the over-all effect of the Kondo screening for different adsorbed atoms. For quantitative estimates, it will be important also to reduce the temperature at which the experiments are performed as far as possible. Irrespective if the Kondo resonance around 8.5 T will eventually be observed or not, the real motivation for  performing such experiments actually comes from the need to estimate the exchange coupling to the bulk electrons, since this parameter also determines the spin relaxation time; thus it plays an important role in the possible applications of high-spin states of magnetic impurities for storing and processing information.
Appendix A. DFT study of Co and Fe adatoms on CuN/Cu(100)  Table A2. The properties of the MLWFs of Co on CuN/Cu(100). We tabulate orbital energy ǫ, occupancy n, orbital spread s, and orbital center along the z-axis. The orbital labels correspond to the initial projections used in the calculation of the corresponding maximally localized Wannier orbitals, although the orbitals are significantly modified from the free-atom orbitals due to the strong hybridization with the neighboring N atoms and the Cu atom just below the adsorbate. The position of the Co adatom in the z-direction is z(Co) =-0.148Å; the origin for the z-axis corresponds to the top-most layer of Cu atoms prior to relaxation.

Co
Minority spin ↓ Majority spin ↑ To further substantiate our claim of the similarities between Co and Fe, which are neighbours in the periodic system differing by a single electron, we have performed a detailed comparative study using the density functional theory (DFT) [42,43]. We used the PWSCF code (from the Quantum Espresso package, http://www.quantum-espresso.org/, Ref. [44]) with a plane-wave basis set, ultra-soft pseudopotentials, and Perdew-Burke-Ernzerhof σ-GGA exchange-correlation functional [45,46,44]. The kinetic-energy cutoff was 35 Ry and the density cut-off 400 Ry. We used an 8x8x4 Monkhorst-Pack mesh of k-points in the self-consistent-field band-structure calculation (but 6x6x1 for the initial relaxation calculation). The cold-smearing by 0.035 Ry has been performed using the Marzari-Vanderbilt scheme [47]. We used the Table A3. The properties of the MLWFs of Fe on CuN/Cu(100). The position of the Fe adatom in the z-direction is z(Fe) =-0.158Å.

Fe
Minority spin ↓ Majority spin ↑ experimentally-determined lattice constant for Cu, a = 0.361 nm. The substrate was modeled using a 2-by-2 supercell in the lateral direction consisting of 4 Cu layers (i.e., 37 atoms in total including the N atoms and the adsorbate). The Cu atoms in the bottommost layer were fixed, while all others were allowed to relax. We note that similar calculations have already been performed for Fe and Mn adatoms [10] and recently also for Mn chains [48] on the same surface.
To gain more insight into the orbital structure of the impurities, we have also calculated the maximally localized Wannier functions (MLWF) [49,50] on the impurity site using the Wannier90 package (http://www.wannier.org/, Ref. [51]). The MLWFs provide a very compact and accurate local representation of the electronic structure: they give direct insight into the nature of the chemical bonding. A 4x4x2 uniform grid of k-points was used here. The initial projections consisted of d and s states on the adatom and the Cu atoms, and of sp3 hybridized states on the N atom; in addition, s states were added in the interstitial sites of the lattice to describe the delocalized electrons (1 per Cu atom). The energy window was chosen so as to encompass all d and s levels of the adatom. Table A1 presents an overview of the results: the magnetic moment and the occupancies of the adatom levels. The results are expected -the two adatoms differ by one electron in the d orbital. More details can be found in Tables A2 and A3 where we show the orbital energies (defined as the average energy of the spectral function of a particular MLWF), occupancies (defined as the integral of the spectral function of MLWF up to the Fermi energy), the centers of the MLWFs, as well as their spreads. It is manifest that the difference consists of one additional electron in the "minority" spin orbitals for the Co adatom. In fact, the spectral functions shown in Fig. A1 reveal that the "majority" spin spectral functions are essentially the same in both cases, while the "minority" spin spectral functions are to a good approximation just rigidly shifted.  Figure A1. Partial density of states (spectral functions) for the MLWF on Co (left) and Fe (right) adatoms. We plot both spin projections: "majority" spin (above), "minority" spin (below). The spectral functions were computed on a 16x16x3 grid of k-points with a 0.1 eV binning resolution.

Appendix B. Calculation of the magnetic susceptibility in anisotropic models
In general, the magnetic susceptibility is calculated as [52,25]: where S α is the α ∈ {x, y, z} component of the spin and β = 1/k B T . The impurity contribution to the total magnetic susceptibility, χ imp,α , is then obtained by subtracting the susceptibility of a system without the impurity, χ (0) α (T ), from the magnetic susceptibility of the full system (i.e., the continuum of the conductance-band electrons plus the impurity), χ full,α (T ): In the problems with U(1) spin symmetry, S z is a conserved quantum number and it is used to classify the invariant subspaces of the Hilbert space, thus the calculation of the expectation values in (B.1) is particularly simple for α = z: Sz;n ) and similarly for S 2 z (T ). Here Z (N ) is the partition function at shell N, defined as Z (N ) = Tr N (e −βH ) = Sz,n exp(−βE (N ) Sz;n ), S z the conserved spin quantum number, n a multi-index composed of all the remaining quantum numbers used to classify the states, and E (N ) Sz;n the energy of the state with given quantum numbers. For other directions α, and for more general problems without any symmetry in the spin space, a full dynamical calculation of the correlator S α (τ )S α (0) has to be performed.
The calculation of χ α requires an iteration of the operator matrix elements for all components of the total spin operator S α . The approach is similar to the calculation of the matrix elements of an arbitrary local impurity operator, but it is extended to allow for total ("global") spin operators S α . The matrix elements at a given step N + 1 are calculated from those at the previous step by using the unitary basis transformation resulting from the diagonalization of the Hamiltonian matrices. We follow the notation of Ref. [25]; see Eqs. (39,40,43) in the cited work: where the matrix elements M N ww are those that correspond to the operator S α . The imaginary-time integral of the correlator C(τ ) = S α (τ )S α (0) may be expressed using the Lehmann representation as Note that for the spin operator along the y direction, the NRG calculation has to be performed using complex numbers. Where comparisons can be made (e.g., in isotropic models) the magnetic susceptibility calculated in the traditional way and using the proposed technique differ by no more than a few per mil in calculations performed for Λ = 2, and even less for higher Λ. We observed, however, that more states need to be kept in this approach than in the standard NRG thermodynamic calculations. The calculation of χ is made possible by the property of the energy-scale separation of total spin operators, much like the entire NRG iteration is based on the energy-scale separation property of the Hamiltonian.
The approach described in this section can be immediately generalized to the calculation of the out-of-diagonal components of the magnetic susceptibility tensor.