Spin excitations and correlations in scanning tunneling spectroscopy

With low-temperature scanning tunneling microscopes (STM) scientist have developed a tool that has the ability to detect and manipulate individual magnetic spin systems on the atomic and molecular level by an externally controlled probe with sub-nm precision. These instruments have opened a new field of research envisioning not only a deeper understanding of the origin of molecular magnetism by studying the interactions between nanoscale spins but also of the many-particle effects between the localized spins and the itinerant electrons of the supporting substrate. The progress in this field is best demonstrated by the recently achieved ability to build stable magnetic bits by cleverly arranging only a handful of Fe atoms on either a thin insulating [1] or on a nonmagnetic Cu(111) surface [2]. While in these experiments the transition from quantum-mechanical to classical magnetic behavior is explored [3], other applications might be envisioned ranging from spin-based logic circuits [4] to entangled systems in which the quantum mechanical nature is crucial for computational purposes [5].

At the basis of all these experiments lies the spectroscopic capability of the STM. About ten years ago Heinrich and his co-workers showed in a hallmark experiment that it is possible to detect inelastic spin-flip excitations on Mn atoms adsorbed on patches of Al₂O₃ on a NiAl surface [6]. Since then many experiments have focused on transition metal atoms on a thin layer of Cu₂N on Cu(100) [1, 7–20]. In these experiments, patches of Cu₂N are formed by sputtering a clean Cu(100) surface with N⁺ ions and subsequent annealing. This leads to a monolayer of Cu₂N on the surface on which the desired metal atoms are deposited (figure 1(a)). Experimentally, every 3d transition metal atom adsorbed on this surface reveals its fingerprint when a spectrum is taken by placing the tunneling tip of an STM over the atom (figure 1(b)). The spectra are measured by varying the bias voltage applied between tip and sample and recording the differential conductance $\rm{d}I/\rm{d}V$. All spectra show characteristic features that can be classified into two groups: (1) step-like increases in the differential conductance, positioned symmetrically around zero bias, and (2) peaks of the differential conductance at zero bias.

Zoom In Zoom Out Reset image size

The number of distinguishable conductance steps varies for different spin systems and the precise form of the steps often shows some asymmetry with respect to the applied bias direction and can exhibit some overshoot of the conductance at the step-energy. As we will see in detail, the steps are due to the opening of additional conductance channels precisely governed by the magnetic properties of the spin system. These excitations can be present even at zero magnetic field due to the magneto-crystalline anisotropy. The anisotropy is caused by the reduction of the geometric symmetry at the surface and by spin–orbit coupling [23–25], which have the effect of lifting the inherent degeneracy of the spin states. For single atoms the maximum anisotropy is limited to a few tens of meV [26], whereby the magnetic excitation energies usually range from less than one meV up to a few meV requiring experimental setups operating at temperatures $T\leqslant 4$ K. Here, it is worth to note that one has to be aware that not all energetically low lying step-like increases in the differential conductance must originate from magnetic excitations. The tunneling electrons can also excite low-energy mechanical vibrations which can produce similar strong spectroscopic features [12, 27–31] and which can interact with the spin excitations [32, 33]. Thus, to clearly distinguish magnetic excitations their behavior in external applied magnetic fields is often crucial.

The peaks at zero bias are due to the Kondo effect in which itinerant electrons from the substrate coherently scatter with the localized magnetic moment of the adatom [34–36]. This effect has first been detected by scanning tunneling spectroscopy on single metal atoms adsorbed on noble metal substrates [37–46]. In these early experiments the metal atoms are relatively strongly coupled to the substrate leading to a characteristic Kondo temperature in the range of ${T}_{K}\approx 30-300$ K, as determined by the full-width at half-maximum of the peak, and in most cases to a strongly asymmetric Fano lineshape due to interference effects with a potential scattering channel [47–49]. A decoupling layer, such as Cu₂N, significantly reduces the Kondo temperature suggesting the possibility of influencing this many-body state with experimentally accessible magnetic fields and temperatures [9, 18, 50, 51]. This enables the study of the interplay between Kondo screening, the magnetic anisotropy, and nearby spins [10, 20].

Apart from Cu₂N, other decoupling substrates have been exploited to study spin excitations. For example, low-lying excitations associated with a spin S = 1 of the magnetic atom have been found for single Fe atoms embedded into the semiconducting InSb(110) surface [52, 53], which exhibits a two-dimensional electron gas confined to the surface. Other experiments used two-dimensional materials, such as graphene and hexagonal boron nitride (h-BN), as substrates for metal adatoms or magnetic molecules. For example, the Kondo state was observed for individual Co atoms adsorbed on a graphene sheet on top of a Ru(0001) surface with different Kondo temperatures and effective g-factors depending on the adsorption site with respect to the underlying metal [50]. In a different experiment Co and CoH_x (x = 1–3) complexes on graphene on top of Pt(111) have been investigated revealing a S = 1 for Co and CoH₃ and a high magnetic anisotropy of $\approx 8.1$ meV for the Co adsorbate [54]. On the highly corrugated and insulating h-BN adlayer on top of a Rh(111) substrate, CoH and CoH₂ complexes showed both, spin-flip excitations and the Kondo effect, pointing to two different spin states of S = 1 and $S=1/2$, respectively [51]. Interestingly, the CoH complex revealed a dependence of the effective anisotropy on the coupling to the underlying substrate, similarly as observed for Co atoms on large patches of Cu₂N [16].

Spin excitations can also be observed for spin systems adsorbed on bare metal substrates, even though lifetime, anisotropy energies, and intensities are in general reduced due to the strong coupling with the substrate. For example, spin-excitations of Co and Fe on Pt(111) have been measured [55, 56], where the measurements of reference [55] were performed at $T\approx 6$ K and showed an approximately ten-times higher apparent magnetic anisotropy as the latter measurement [56] that was performed at T = 0.3 K. This discrepancy stems from the strong temperature dependent broadening of the spectrum, which leads to this gross overestimation and illustrates the need for low-temperature measurements [57]. Additionally, spin excitations have been detected for Fe atoms adsorbed on Cu(111) [58] and on Ag(111) [59]. Apart from the 3d transition metal atoms, 4f lanthanide atoms also show low-energy magnetic signals. For example, the Kondo effect was measured on small Ce clusters on Cu₂N [12] and spin-flip signals were observed for Gd adsorbed on Pt(111) and Cu(111) [60] and Ho on Pt(111) [61], even though the intensities of the spin-flip excitations are unclear but presumably very small due to the strong localization of the 4f wavefunctions close to the nuclei and their small spacial extension [62].

In addition to these, metal-organic complexes, such as M-phthalocyanine with $M=$ Mn, Fe, Co, Ni, and Cu, have been studied on different surfaces [63–70]. Many of them showed Kondo screening with a lower Kondo temperature enabling the tuning by dehydrogenation [63, 70] or the splitting of the peak by an externally applied magnetic field [64, 68] or by coupling to a ferromagnetic substrate [71]. Interestingly, FePc adsorbed on Au(111) showed a clear Kondo signature [69], while the same molecule adsorbed on top of a CuO decoupling layer showed a double-step in the differential conductance pointing to S = 1 [65].

The electronic decoupling can also be created by using a substrate in the superconducting state in which the creation of Cooper pairs leads to a strongly reduced density of unpaired electrons around the Fermi energy. This results in a significant increase of the spin lifetime of the metal-organic complex [72]. Furthermore, such a substrate allows one to study the competition between superconducting phenomena and Kondo screening as observed on MnPc adsorbed on Pb(111) [67, 73].

Apart from metal atoms and metal-organic complexes spin-flip excitations were detected in complex molecules in which several spin centers are coupled forming a giant spin, as it has been shown at the example of the prototypical molecular magnet manganese-12-acetate-16 which has a ground state with a total spin of S = 10 [21] (figures 1(c) and (d)). Here, the h-BN decoupling layer is crucial as the magnetism is strongly quenched when the molecule is adsorbed upon a Au(111) surface. Furthermore, fully organic molecules including the charge-transfer complexes TTF-TCNQ (tetrathiafulvalene-tetracyanoquinodimethane) [32] and an organic radical molecule [22], have shown the Kondo effect, whereby for the latter molecule the temperature and magnetic field dependence was fully understood and simulated in a perturbative approach (figures 1(e) and (f)).

To describe the inelastic spin-flip excitation spectra, scattering models have been developed that treat the interaction of the tunneling electron with the localized spin system effectively as a one-electron second-order perturbation [74–77]. Additionally, similar models allowed rationalization of the change in the spectra when a spin-polarized tip is employed [78] and the dynamics at increased coupling between tip and sample [79–81]. To address the experimentally observed bias asymmetry, co-tunneling models in the strong Coulomb-blockade regime have been proposed [82] and many-electron effects using the non-crossing approximation have been included [83]. Furthermore, third-order scattering models similar to the ones used here were employed to well describe the observed bias overshooting at certain inelastic excitation steps [84, 85]. Nevertheless, these models were restricted to Kondo-like interactions and did not include potential scattering.

In this paper, I plan to review the straightforward model of the exchange-interaction between an isolated spin system and the tunneling electrons and use a perturbative tunneling approach to simulate experimental data with unprecedented accuracy. The basic idea of the model goes back 50 years to the hallmark discovery of Juan Kondo that higher order scattering of bulk electrons on a magnetic impurity leads to logarithmic divergences [34, 35]. We will see that such a model enables the determination of physical properties like the magnetocrystalline anisotropy, the coupling strength to the substrate, the state lifetimes, and the magnetization directly by comparing the differential conductance measurements with simulations. Furthermore, it allows one to grasp some of the correlations and entanglements that are formed by the many-particle interactions due to the large electron bath of the substrate. The model, even though inherently limited due to its perturbative approach, has the advantage that we can develop it straightforwardly using simple matrix algebra. Due to this simplicity it is computationally very fast and can thus help us to understand future experimental observations. In these calculations the magnetic anisotropies, gyromagnetic factors, and coupling strengths to the substrate enter as experimentally determined parameters. Note however, that in particular for transition metal atoms adsorbed on metallic substrates the magnetic anisotropy have been successfully determined from first principle by, for example, time-dependent density functional theory [56, 86–88] or density functional theory which includes spin–orbit coupling [89, 90].

The basic idea of the model is to treat the tunneling as a scattering event of an incoming electron with the localized spin system (figure 2(a)). Postponing the exact meaning for later, it is amazing that only two parameters govern these interactions: the Kondo scattering parameter J and the potential scattering parameter U. These two parameters are crucial for the understanding of the low-energy excitations. One is tempted to ask 'What is the origin of these?' To find an answer it is helpful to step back from the scattering picture and look at the single-impurity Anderson model, which describes one half-filled atomic orbital (for example from a 3d transition metal atom) interacting with the continuous states of a host metal (figure 2(b)) [91]. In this picture all energies are in the eV range and it is not per se clear how this relates to the scattering model with low-energy excitations in the meV range. These relations were found by Schrieffer and Wolf, who were able to relate the high energies of the Anderson model with the scattering parameters [92]:

$$\begin{eqnarray}&&J={\displaystyle \frac{2}{\pi }}{\Delta }^{2}\left({\displaystyle \frac{1}{{\epsilon }_{d}}}-{\displaystyle \frac{1}{{\epsilon }_{d}+{U}_{d}}}\right)\lt 0,\quad \quad \mathcal{U}={\displaystyle \frac{2}{\pi }}{\Delta }^{2}\left({\displaystyle \frac{1}{{\epsilon }_{d}}}+{\displaystyle \frac{1}{{\epsilon }_{d}+{U}_{d}}}\right).\end{eqnarray} \tag{ 1 }$$

Here it is crucial for the understanding of the forthcoming model that the Kondo spin–spin exchange scattering J is for $S=1/2$ systems always negative, which means that the localized spin and the substrate electrons couple antiferromagnetically, while the potential scattering $\mathcal{U}$ can have both signs, either positive or negative. However, this result can be different for higher spins.

Zoom In Zoom Out Reset image size

The manuscript is organized as follows: in section 2 the basic model in the zero-current approximation is introduced and discussed for the example spin $S=1/2$. In section 3 high-spin systems are discussed. Experimentally observed examples are given and compared to the model calculations. Section 4 discusses the limitations of the model, in particular, its inability to cover all correlation effects that occur when a system enters the strong Kondo regime using experimental data obtained on the Co/Cu₂N system as an example. In section 5 the initial model is extended by including rate equations and tested against experimental data. Here, we will observe that the tunneling current too can lead to the appearance of a non-equilibrium Kondo effect. Finally, section 6 summarizes the manuscript and outlines possible routes for future extensions.

To describe the experimental observations we use a simplified model in which the Hamiltonian of the total quantum mechanical system is divided into the ones of the subsystems of the two electron reservoirs in tip and sample, the one of the localized spin system, and an interaction Hamiltonian that enables the exchange of charge carriers between the reservoirs:

$$\begin{eqnarray}&&\hat{H}={\hat{H}}_{\rm{t}}+{\hat{H}}_{\rm{s}}+{\hat{H}}_{\rm{a}}+{\hat{H}}^{\prime }.\end{eqnarray} \tag{ 2 }$$

The Hamiltonian of the tip ${\hat{H}}_{\rm{t}}$ and sample ${\hat{H}}_{\rm{s}}$, respectively, can be described in the framework of second quantization using creation and annihilation operators ${\hat{a}}^{\dagger }$ and $\hat{a}$:

$$\begin{eqnarray}&&{\hat{H}}_{\rm{t}}={\displaystyle \sum }_{k,\sigma }{\epsilon }_{\rm{t}k\sigma }{\hat{a}}_{\rm{t}k\sigma }^{\dagger }{\hat{a}}_{\rm{t}k\sigma },\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\hat{H}}_{\rm{s}}={\displaystyle \sum }_{k,\sigma }{\epsilon }_{\rm{s}k\sigma }{\hat{a}}_{\rm{s}k\sigma }^{\dagger }{\hat{a}}_{\rm{s}k\sigma },\end{eqnarray} \tag{ 4 }$$

with ${\epsilon }_{\rm{t}k\sigma }$ and ${\epsilon }_{\rm{s}k\sigma }$ as the energy of the electrons with momentum k and spin σ in the tip and sample, respectively.

These two many-particle systems are the source and sink for the tunneling electrons in the STM experiment. Instead of using creation and annihilation operators in the momentum space k, we assume in the small energy range of interest, i.e. to some tens of meV around the Fermi energy, for tip and sample a continuous and energetically flat density of states $\rho (\epsilon )={{\displaystyle \sum }}_{\sigma ,k(\epsilon ),k\prime (\epsilon )}{\delta }_{k,k\prime }\left\langle{\hat{a}}_{k\sigma }^{\dagger }{\hat{a}}_{k\prime \sigma }\right\rangle\equiv \rho$, with $\langle \cdot \rangle$ as the time averaged expectation value.

In general, the electronic states $| \varphi \rangle$ in these electrodes might be spin-polarized in an arbitrary direction which we account for by the corresponding spin density matrices $\varrho =| \varphi \rangle \langle \varphi | ={\displaystyle \frac{1}{2}}(\hat{I}+\vec{n}\cdot \hat{\boldsymbol{\sigma }})$. Here, $\vec{n}$ describes the direction and $-1\leqslant | \vec{n}| \leqslant 1$ the amplitude of the polarization in the chosen coordinate system, $\hat{\boldsymbol{\sigma }}=({\hat{\sigma }}_{x},{\hat{\sigma }}_{y},{\hat{\sigma }}_{z})$ are the standard Pauli matrices, and $\hat{I}$ is the $(2\times 2)$ identity matrix. With this convention the spin polarization is identical to the relative imbalance between majority and minority spin densities $\eta =| \vec{n}| =\left|{\displaystyle \frac{{\rho }_{\uparrow }-{\rho }_{\downarrow }}{{\rho }_{\uparrow }+{\rho }_{\downarrow }}}\right|$, where $\uparrow$ and $\downarrow$ account for the two different spin directions along the chosen quantization axis. This description via density matrices allows arbitrary polarization directions in tip and sample which obey the quantum statistics (see appendix).

The impurity spin system may either contain only a single spin or a finite number of coupled spins. We describe this system by a model Hamiltonian ${\hat{H}}_{\rm{a}}$ that includes Zeeman and magnetic anisotropy energy and—in the case of more than one spin—the Heisenberg coupling terms ${\vec{\mathcal{J}}}_{{ij}}$ ¹ and the non-collinear Dzyaloshinskii–Moriya coupling ${\vec{\mathcal{D}}}_{{ij}}$ between individual spins:

$$\begin{eqnarray}\begin{array}{ccc}{\hat{H}}_{\rm{a}} & = & {\displaystyle \sum }_{i}{g}_{i}{\mu }_{B}\vec{B}\cdot {\hat{\bf{S}}}^{i}+{\displaystyle \sum }_{i}\left({D}_{i}{({\hat{S}}_{z}^{i})}^{2}+{E}_{i}\left[{({\hat{S}}_{x}^{i})}^{2}-{({\hat{S}}_{y}^{i})}^{2}\right]\right)\\ & & +{\displaystyle \sum }_{i,j}{\vec{\mathcal{J}}}_{{ij}}{\hat{\bf{S}}}^{i}\cdot {\hat{\bf{S}}}^{j}+{\displaystyle \sum }_{i,j}{\vec{\mathcal{D}}}_{{ij}}{\hat{\bf{S}}}^{i}\times {\hat{\bf{S}}}^{j}.\end{array}\end{eqnarray} \tag{ 5 }$$

In this equation ${\mu }_{B}$ is the Bohr magneton and g_i, D_i, and E_i are the gyromagnetic factor, the axial and the transversal magnetic anisotropy for the ith spin in the reference coordinate frame, respectively² . The externally applied magnetic field is $\vec{B}$ and the total spin operators ${\hat{\bf{S}}}^{i}=({\hat{S}}_{x}^{i},{\hat{S}}_{y}^{i},{\hat{S}}_{z}^{i})$ are built from operators of the form ${\hat{S}}_{x,y,z}^{i}={\hat{I}}_{1}\otimes \cdots \otimes {\hat{S}}_{x,y,z}^{i}\otimes \cdots \otimes {\hat{I}}_{n}$, which only act on the ith spin and where $\hat{I}$ denotes the identity matrix and ⨂ the Kronecker matrix product. The spin operators ${\hat{S}}_{x,y,z}^{i}$ can be easily constructed remembering that ${\hat{S}}_{z}| {\psi }_{m}\rangle ={\hslash }{m}_{z}| {\psi }_{m}\rangle$ with ${m}_{z}=S,S-1,\ldots ,-S$ as the magnetic quantum number along our chosen z-axis and ${\hat{S}}_{\pm }| {\psi }_{m}\rangle ={\hslash }\sqrt{(S\mp m)(m\pm S+1)}| {\psi }_{m\pm 1}\rangle$. These then enable calculation of ${\hat{S}}_{x}={\displaystyle \frac{1}{2}}({\hat{S}}_{+}+{\hat{S}}_{-})$ and ${\hat{S}}_{y}=-{\displaystyle \frac{1}{2}}\rm{i}({\hat{S}}_{+}-{\hat{S}}_{-})$, respectively. For simplicity we will set ${\hslash }=1$ from now on.

Diagonalizing the Hamiltonian of the localized spin system (equation (5)) leads to discrete energy eigenvalues ${\epsilon }_{n}$ and eigenstates $| \psi {\rangle }_{n}$. For a single spin there are $2S+1$ eigenstates, while for complex spin structures that consist of several spins the number of eigenstates increases quickly as ${\prod }_{i}(2{S}_{i}+1)$. Because we neglect direct interactions between the electron baths in tip, sample and the localized spin system, we describe the total state as a product of the continuous electron states $| \varphi \rangle$ and the discrete spin states $| \psi \rangle$, i.e. $| {\varphi }^{\rm{t},\rm{s}},\psi \rangle$ $\;=| {\varphi }^{\rm{t},\rm{s}}\rangle | \psi \rangle$. Note, however, that the coupling of the spin-system with the substrate will lead to an entanglement between sample electrons and the spin system (see section 4) and to a renormalization of the parameters in equation (5) [16, 51, 93, 94]. Here we assume that the renormalization is already included in the anisotropies and gyromagnetic factors and omit for the moment the entanglement.

The interaction Hamiltonian H' of equation (2) allows for tunneling of electrons from tip to sample or vice versa only via Kondo-like spin-flip or potential scattering processes with the impurity (figure 3(a)):

$$\begin{eqnarray}\begin{array}{ccc}{\hat{H}}^{\prime } & = & {\hat{V}}_{\rm{t}\to \rm{s}}+{\hat{V}}_{\rm{s}\to \rm{t}}+{\hat{V}}_{\rm{s}\to \rm{s}},\\ {\hat{V}}_{\rm{t}\to \rm{s}} & = & {\displaystyle \sum }_{i,{\lambda }^{\prime },\lambda }{T}_{\mathrm{ta}}^{i}{J}_{i}\;{\hat{a}}_{\rm{s}{\lambda }^{\prime }}^{\dagger }\;{\hat{a}}_{\rm{t}\lambda }\left({\displaystyle \frac{1}{2}}\hat{\boldsymbol{\sigma }}\cdot {\hat{\bf{S}}}^{i}+{U}_{i}{\hat{I}}_{i}\right),\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\hat{V}}_{\rm{s}\to \rm{t}}={\displaystyle \sum }_{i,{\lambda }^{\prime },\lambda }{T}_{\mathrm{ta}}^{i}{J}_{i}\;{\hat{a}}_{\rm{t}{\lambda }^{\prime }}^{\dagger }\;{\hat{a}}_{\rm{s}\lambda }\left({\displaystyle \frac{1}{2}}\hat{\boldsymbol{\sigma }}\cdot {\hat{\bf{S}}}^{i}+{U}_{i}{\hat{I}}_{i}\right),\end{eqnarray} \tag{ 7 }$$

with T_taⁱ as the coupling constants between tip and the ith adsorbate spin, and ${U}_{i}={\mathcal{U}}_{i}/{J}_{i}$ as the unitless ratio between Kondo and potential scattering. While we assume that the impurity is much more strongly coupled to the sample than to the tip, spin–spin scattering between the impurity and substrate electrons are additionally considered:

$$\begin{eqnarray}&&{\hat{V}}_{\rm{s}\to \rm{s}}={\displaystyle \frac{1}{2}}{\displaystyle \sum }_{i,{\lambda }^{\prime },\lambda }{J}_{i}\ {\hat{a}}_{\rm{s}{\lambda }^{\prime }}\;{\hat{a}}_{\rm{s}\lambda }^{\dagger }\;\hat{\boldsymbol{\sigma }}\cdot {\hat{\bf{S}}}^{i}.\end{eqnarray} \tag{ 8 }$$

Zoom In Zoom Out Reset image size

With the above assumptions the model is largely identical to the ones studied by Appelbaum and Anderson already in the 70s for mesoscopic tunnel junctions [95–97].

2.1. Current and conductance in second order

We will now briefly review the calculation of the tunneling current using Fermi's golden rule. In STM experiments the bias eV applied between tip and sample shifts the Fermi level ${\epsilon }_{F}$ of the tip with respect to the one of the sample. At positive electrons from occupied tip states can cross the tunnel-barrier and interact with the localized spin system under exchange of angular momentum and energy (figure 3(b)). To obey energy conservation, the energy difference ${\epsilon }_{{if}}={\epsilon }_{f}-{\epsilon }_{i}$ between initial and final state of the spin system has to be accounted for by the energy of the tunneling electron. The current flowing between tip and sample is then given by:

$$\begin{eqnarray}\begin{array}{ccc}{I}^{\rm{t}\to \rm{s}}({eV}) & = & {\displaystyle \frac{2\pi e}{{\hslash }}}{T}_{0}^{2}{\displaystyle \sum }_{i,f}{p}_{i}{{\displaystyle \int }}_{-\infty }^{+\infty }{\left|{\mathcal{M}}_{{if}}^{\rm{t}\to \rm{s}}\right|}^{2}f(\epsilon -{eV})\left[1-f(\epsilon -{\epsilon }_{{if}})\right]\rm{d}\epsilon ,\\ {I}^{\rm{s}\to \rm{t}}({eV}) & = & {\displaystyle \frac{2\pi e}{{\hslash }}}{T}_{0}^{2}{\displaystyle \sum }_{i,f}{p}_{i}{{\displaystyle \int }}_{-\infty }^{+\infty }{\left|{\mathcal{M}}_{{if}}^{\rm{s}\to \rm{t}}\right|}^{2}f(\epsilon +{eV})\left[1-f(\epsilon -{\epsilon }_{{if}})\right]\rm{d}\epsilon ,\\ I & = & {I}^{\rm{t}\to \rm{s}}-{I}^{\rm{s}\to \rm{t}},\end{array}\end{eqnarray} \tag{ 9 }$$

where we have dropped the additional summation that would account for current through the different sites in coupled spin systems, i.e. we restrict our calculation to single spin systems.

The unitless tunnel barrier transmission coefficient ${T}_{0}^{2}=| {T}_{\mathrm{ta}}J{\rho }_{\rm{s}}{| }^{2}$ contains all experimental parameters that will determine the strength of the tunneling current. In particular, T₀ strongly depends on the distance between tip and adsorbate, and on the spin independent local densities of states in tip and sample. For the moment we want to assume that ${T}_{0}^{2}\ll 1$ so that the influence of the tip on the spin system is negligible and the time between consecutive tunneling events is much longer than the relaxation time of the spin system, e.g. via processes such as described by equation (8). Such conditions are usually easy to provide in STM experiments where the tip–sample separation can be adjusted to give tunneling currents in the pA range or below so that the tip can be seen as a local probe of the spin system that leaves the spins always in thermal equilibrium with the substrate (zero-current approximation).

Under this assumption the average state occupation p_i(T) of the spin system is only governed by the effective temperature T and follows the Boltzmann distribution,

$$\begin{eqnarray}&&{p}_{i}(T)={\displaystyle \frac{\mathrm{exp}[-{\epsilon }_{i}/({k}_{B}T)]}{{{\displaystyle \sum }}_{i}\mathrm{exp}[-{\epsilon }_{i}/({k}_{B}T)]}},\end{eqnarray} \tag{ 10 }$$

with k_B as the Boltzmann constant. Furthermore, we assume a flat density of states in tip and sample so that the electron occupation is given by the Fermi–Dirac distribution $f(\epsilon )={[1+\mathrm{exp}(\epsilon /{k}_{B}T)]}^{-1}$.

While most STM experiments do not record the tunneling current directly, we are interested in its derivative with respect to the bias voltage. With the energy independent matrix elements the derivative of the current is easy to evaluate. Integrating the Fermi–Dirac distributions of equation (9) and calculating the derivative results in a temperature broadened step function [98]:

$$\begin{eqnarray}&&\Theta (x)={\displaystyle \frac{1+(x-1)\mathrm{exp}(x)}{\mathrm{exp}{(x)}^{2}}},\end{eqnarray} \tag{ 11 }$$

with $x=\epsilon /({k}_{B}T)$, so that equation (9) becomes:

$$\begin{eqnarray}&&{\displaystyle \frac{\partial I}{\partial V}}{({eV})}^{\rm{t}\to \rm{s}}={\displaystyle \frac{2\pi {e}^{2}}{{\hslash }}}{T}_{0}^{2}{\displaystyle \sum }_{i,f}{p}_{i}{\left|{\mathcal{M}}_{{if}}^{\rm{t}\to \rm{s}}\right|}^{2}\Theta ({eV}-{\epsilon }_{{if}}).\end{eqnarray} \tag{ 12 }$$

The main goal of this paper lies in deriving the transition matrix elements ${\mathcal{M}}_{{if}}$ to be able to solve the equation (12). First, we start with second order processes, neglecting scattering that involves electrons which originate and end in the substrate bath, i.e. we concentrate on tunneling processes which consist of only one scattering event as depicted in figure 3(b) and described by the transport Hamiltonian of equations (6) and (7):

$$\begin{eqnarray}&&{\mathcal{M}}_{{if}}^{(1)}=\left.\langle {\varphi }_{f},{\psi }_{f}\right|{\displaystyle \frac{1}{2}}\hat{\boldsymbol{\sigma }}\cdot \hat{\bf{S}}+U\hat{I}\left|{\varphi }_{i},{\psi }_{i}\right.\rangle={M}_{{if}}+U{\delta }_{{if}}.\end{eqnarray} \tag{ 13 }$$

This matrix is energy independent and connects the initial states in the electron baths $| {\varphi }_{i}\rangle$ and spin system $| {\psi }_{i}\rangle$ with their final states $| {\varphi }_{f}\rangle$ and $| {\psi }_{f}\rangle$. It contains the spin-exchange scattering term M_if and a potential scattering term that has non-zero matrix-elements only when initial and final spin state are the same (${\delta }_{{if}}$ is the delta distribution). Computing the absolute square of ${\mathcal{M}}_{{if}}^{(1)}$ results in three terms that contribute to the tunneling conductance:

$$\begin{eqnarray}&&{| {\mathcal{M}}_{{if}}^{(1)}| }^{2}={| {M}_{{if}}| }^{2}+{| U| }^{2}{\delta }_{{if}}+2\mathfrak{R}\left(U\times {M}_{{if}}\right){\delta }_{{if}},\end{eqnarray} \tag{ 14 }$$

where $\mathfrak{R}(A)$ denotes the real part of the matrix A. The first term consist of the operators

$$\begin{eqnarray}&&\hat{\boldsymbol{\sigma }}\cdot \hat{\bf{S}}={\hat{\sigma }}_{x}{\hat{S}}_{x}+{\hat{\sigma }}_{y}{\hat{S}}_{y}+{\hat{\sigma }}_{z}{\hat{S}}_{z}={\displaystyle \frac{1}{2}}{\hat{\sigma }}_{+}{\hat{S}}_{-}+{\displaystyle \frac{1}{2}}{\hat{\sigma }}_{-}{\hat{S}}_{+}+{\hat{\sigma }}_{z}{\hat{S}}_{z},\end{eqnarray} \tag{ 15 }$$

which connect the initial and final state and accounts for spin-exchange processes in which angular momentum between the tunneling electron and the localized spin system can be exchanged. The second term accounts for potential scattering between the tunneling electron and the localized spin system and does not change the spin state. The third term results from the interference between potential and spin-exchange scattering and depends, as we will see, strongly on the angular momentum of the localized spin and is the origin of magneto-resistive tunneling [13, 99].

Due to the product-state of the total quantum-mechanical system, it is worth mentioning that the matrix elements have to be independently evaluated for the localized spin and for the tunneling electron. For an electron tunneling between a spin-average tip and sample the non-zero matrix elements are:

$$\begin{eqnarray}&&\langle \uparrow | {\hat{s}}_{+}| \downarrow \rangle \;=\;+1,\end{eqnarray} \tag{ 16a }$$

$$\begin{eqnarray}&&\langle \downarrow | {\hat{s}}_{-}| \uparrow \rangle \;=\;+1,\end{eqnarray} \tag{ 16b }$$

$$\begin{eqnarray}&&\langle \uparrow | {\hat{s}}_{z}| \uparrow \rangle \;=\;+1/2,\end{eqnarray} \tag{ 16c }$$

$$\begin{eqnarray}&&\langle \downarrow | {\hat{s}}_{z}| \downarrow \rangle =-1/2,\end{eqnarray} \tag{ 16d }$$

$$\begin{eqnarray}&&\langle \uparrow | \hat{I}| \uparrow \rangle \;=\;+1,\end{eqnarray} \tag{ 16e }$$

$$\begin{eqnarray}&&\langle \downarrow | \hat{I}| \downarrow \rangle \;=\;+1,\end{eqnarray} \tag{ 16f }$$

where, for completeness, we have additionally included the probability amplitudes for interacting with the unity operator $\hat{I}$. Since $\hat{\sigma }=2\hat{s}$, these matrix elements enable us to rewrite the spin-exchange scattering term in equation (14) for a spin-average tip and sample to [8, 99]:

$$\begin{eqnarray}&&{| {M}_{{if}}| }^{2}={\displaystyle \frac{1}{2}}{\left|\left.\langle {\psi }_{f}\right|{\hat{S}}_{-}\left|{\psi }_{i}\right.\rangle\right|}^{2}+{\displaystyle \frac{1}{2}}{\left|\left.\langle {\psi }_{f}\right|{\hat{S}}_{+}\left|{\psi }_{i}\right.\rangle\right|}^{2}+{\left|\left.\langle {\psi }_{f}\right|{\hat{S}}_{z}\left|{\psi }_{i}\right.\rangle\right|}^{2}.\end{eqnarray} \tag{ 17 }$$

Note, that (16c) and (16d) do not cancel out because the two different spin directions in the initial (tip or sample) and final (sample or tip) bath are incoherent ensemble states and therefore cannot interfere with each other.

For a tip with spin polarization $\eta =\left|\frac{{\rho }_{\uparrow }-{\rho }_{\downarrow }}{{\rho }_{\uparrow }+{\rho }_{\downarrow }}\right|$ along the z-axis, the transition matrix elements are tunnel-direction dependent because either the initial or the final electron bath has now a spin imbalance in the density of states [99]:

$$\begin{eqnarray}\begin{array}{ccc}{| {M}_{{if}}^{\rm{t}\to \rm{s}}| }^{2} & = & {\displaystyle \frac{1+\eta }{2}}{\left|\left.\langle {\psi }_{f}\right|{\hat{S}}_{-}\left|{\psi }_{i}\right.\rangle\right|}^{2}+{\displaystyle \frac{1-\eta }{2}}{\left|\left.\langle {\psi }_{f}\right|{\hat{S}}_{+}\left|{\psi }_{i}\right.\rangle\right|}^{2}+{\left|\left.\langle {\psi }_{f}\right|{\hat{S}}_{z}\left|{\psi }_{i}\right.\rangle\right|}^{2},\\ {| {M}_{{if}}^{\rm{s}\to \rm{t}}| }^{2} & = & {\displaystyle \frac{1-\eta }{2}}{\left|\left.\langle {\psi }_{f}\right|{\hat{S}}_{-}\left|{\psi }_{i}\right.\rangle\right|}^{2}+{\displaystyle \frac{1+\eta }{2}}{\left|\left.\langle {\psi }_{f}\right|{\hat{S}}_{+}\left|{\psi }_{i}\right.\rangle\right|}^{2}+{\left|\left.\langle {\psi }_{f}\right|{\hat{S}}_{z}\left|{\psi }_{i}\right.\rangle\right|}^{2}.\end{array}\end{eqnarray} \tag{ 18 }$$

As an example figure 4 shows the spectrum for a single spin with $S=1/2$. When applying the magnetic field, the conductance increases step-like at bias energies above the Zeeman-energy $| {eV}| \gt g{\mu }_{B}| B|$. These steps are symmetrical around zero voltage and thermally smeared out by about $5{k}_{B}T$ [98]. When the electrons in the tip are spin-polarized the step heights at positive and negative bias are different. For a paramagnetic tip with a spin polarization in direction of the applied field ($\eta \gt 0$) the conductance step at positive bias decreases, while it increases by the same amount at negative bias (figure 4(b)).

Zoom In Zoom Out Reset image size

With a spin-polarized tip the interference between the potential and Kondo-scattering shows a very particular effect [13, 99]. Depending on the sign of U the zero-bias conductance is either increased or reduced and can become even zero at $\eta =1$ and $U=-{\displaystyle \frac{1}{2}}S$ (figures 4(c)–(d)). This magneto-resistive elastic tunneling [99] arises from the third term in equation (14) and scales with the expectation value $\langle {\hat{S}}_{z}\rangle$ of the impurity. At non-equilibrium the strength and sign of this term can change; this allows to read-out of the z-projection of the spin without exciting it, as has been shown in spin-pumping (see section 5) [13, 99] and pump-probe experiments [14].

2.2. Expansion to third order

We now turn our attention to the next higher order of interaction. We want to consider processes in which, during the transport of an electron from tip to sample (or vice versa), the spin system additionally interacts with an electron from the sample bath (figure 3(c)). These processes involve an intermediate state (m) and are expressed in second order Born approximation as:

$$\begin{eqnarray}&&{\mathcal{M}}_{{if}}^{(2)}=J{\rho }_{\rm{s}}\underset{m}{{\displaystyle \sum\!\!\!\!\!\!\!\!\int}}\left({\displaystyle \frac{{\tilde{M}}_{{mf}}{M}_{{im}}}{{\epsilon }_{i}-{\epsilon }_{m}}}+{\displaystyle \frac{{M}_{{mf}}{\tilde{M}}_{{im}}}{{\epsilon }_{m}-{\epsilon }_{f}}}\right).\end{eqnarray} \tag{ 19 }$$

Here, the tilde on the matrix $\tilde{M}$ tags scattering processes between the localized spin and sample electrons only. The order of the two different electron-spin interactions can be exchanged. Thus, we have also to account for processes in which first an electron that originates and ends in the sample scatters the system into the intermediate state and then the tunneling electron interacts with the system scattering it into its final state (right fraction in equation (19)). Here, it is of fundamental importance to have in mind that, contrary to the initial and final state of the total system, the intermediate state $| {\psi }_{m}\rangle$ can be virtual, i.e., must not necessarily obey angular momentum and energy conservation. The integro-summation symbol in equation (19) indicates that we perform both, a summation over all possible discrete intermediate states in the local spin system $| \psi \rangle$ and an integration over the continuous states of the intermediate electron states $| \varphi \rangle$ in the substrate. For the integration, we consider states in an energy range of $\pm {\omega }_{0}$ around E_F as possible scatterer leading to the following characteristic function for electron-like processes [97, 100]:

$$\begin{eqnarray}&&{F}_{e}(\epsilon ,T)=-{{\displaystyle \int }}_{-\infty }^{+\infty }\rm{d}{\epsilon }^{\prime\prime }{{\displaystyle \int }}_{-{\omega }_{0}}^{+{\omega }_{0}}\rm{d}{\epsilon }^{\prime }{\displaystyle \frac{1-f({\epsilon }^{\prime },T)}{{\epsilon }^{\prime }-{\epsilon }^{\prime\prime }}}f\prime ({\epsilon }^{\prime\prime }-\epsilon ,T),\end{eqnarray} \tag{ 20 }$$

and for hole-like processes:

$$\begin{eqnarray}&&{F}_{h}(\epsilon ,T)=-{{\displaystyle \int }}_{-\infty }^{+\infty }\rm{d}{\epsilon }^{\prime\prime }{{\displaystyle \int }}_{-{\omega }_{0}}^{+{\omega }_{0}}\rm{d}{\epsilon }^{\prime }{\displaystyle \frac{f({\epsilon }^{\prime },T)}{{\epsilon }^{\prime }-{\epsilon }^{\prime\prime }}}f\prime ({\epsilon }^{\prime\prime }-\epsilon ,T),\end{eqnarray} \tag{ 21 }$$

respectively. Here, the Fermi–Dirac distributions in the numerator ensure an unoccupied (occupied) final state and the integration over the derivative f ' accounts for the temperature broadening during the total process (figure 5(a)). The switching from virtual to real processes at $\epsilon =0$ leads to a logarithmic singularity that is broadened by the temperature. Thus, as long as the energy ${\epsilon }_{m}$ of the intermediate state and the temperature k_BT are small compared to the integration bandwidth ${\omega }_{0}$, the equations (20) and (21) (with a change of sign) can be rewritten as:

$$\begin{eqnarray}&&F({eV}-{\epsilon }_{m})\approx -{{\displaystyle \int }}_{-\infty }^{+\infty }\rm{d}{\epsilon }^{\prime }\mathrm{ln}\left({\displaystyle \frac{{\omega }_{0}+| {eV}-{\epsilon }_{m}| }{{eV}-{\epsilon }_{m}+\rm{i}{\Gamma }_{0}}}\right)\times {\Theta }^{\prime }({eV}-{\epsilon }_{m}+{\epsilon }^{\prime },T).\end{eqnarray} \tag{ 22 }$$

Here, ${\Theta }^{\prime }(\epsilon ,T)=\partial \Theta (\epsilon ,T)/\partial \epsilon$ is the derivative of the temperature broadened step function (equation (11)) and a small ${\Gamma }_{0}$ accounts for additional non-thermal (lifetime) broadening. Figure 5 shows the energy and temperature dependence of F which has first been experimentally observed on a radical molecule adsorbed on a Au(111) surface (see figure 1(e) and f) [22]. Compared to the Lorentzian or Frota functions [101] normally used for describing Kondo resonances in STM measurements [12], the behavior of this function is quite different, even though—if analyzed only in a small energy range of a few k_BT around zero—, it has a similar shape [22]. The function shows a relatively flat top whose width is determined only by the temperature. For energies ${k}_{B}T\lt | \epsilon | \lt {\omega }_{0}$ the amplitude decays logarithmically as $-\mathrm{ln}(| \epsilon /{\omega }_{0}| )$ while the maximum peak height is proportional to $-\mathrm{ln}({k}_{B}T/{\omega }_{0})$. Due to the restriction to states in an energy interval $\pm {\omega }_{0}$ the function remains analytical even at $\epsilon \to \pm \infty$ where F approaches zero. Note, that the precise value of the cut-off energy ${\omega }_{0}$ is not critical and mainly changes only the background offset. If not otherwise noted we use ${\omega }_{0}=20$ meV throughout the paper.

Zoom In Zoom Out Reset image size

When calculating the conductance we now have to consider both processes depicted in figures 3(b) and (c) and thus have to replace ${\mathcal{M}}_{1}$ with ${\mathcal{M}}_{1}+{\mathcal{M}}_{2}$ in equation (12) leading to:

$$\begin{eqnarray}\begin{array}{ccc}{\left|{\mathcal{M}}_{{if}}^{(1)}+{\mathcal{M}}_{{if}}^{(2)}\right|}^{2} & = & {\left|{{M}}_{{if}}^{(1)}\right|}^{2}\\ & & +J{\rho }_{\rm{s}}{\displaystyle \sum }_{m}\left[\left({M}_{{fi}}{\tilde{M}}_{{mf}}{M}_{{im}}+\bar{{M}_{{fi}}{\tilde{M}}_{{mf}}{M}_{{im}}}\right)F({\epsilon }_{{mi}})\right.\\ & & +\left.\left({M}_{{fi}}{M}_{{mf}}{\tilde{M}}_{{im}}+\bar{{M}_{{fi}}{M}_{{mf}}{\tilde{M}}_{{im}}}\right)F({\epsilon }_{{im}})\right]\end{array}\end{eqnarray} \tag{ 23a }$$

$$\begin{eqnarray}\begin{array}{ccc} & & +J{\rho }_{\rm{s}}U{\displaystyle \sum }_{m}\left[\left({I}_{{fi}}{\tilde{M}}_{{mf}}{M}_{{im}}+\bar{{I}_{{fi}}{\tilde{M}}_{{mf}}{M}_{{im}}}\right)F({\epsilon }_{{mi}})\right.\\ & & +\left.\left({I}_{{fi}}{M}_{{mf}}{\tilde{M}}_{{im}}+\bar{{I}_{{fi}}{M}_{{mf}}{\tilde{M}}_{{im}}}\right)F({\epsilon }_{{im}})\right]\end{array}\end{eqnarray} \tag{ 23b }$$

$$\begin{eqnarray}&&+O{\left(J{\rho }_{\rm{s}}\right)}^{2}.\end{eqnarray}$$

The evaluation of the matrix elements up to third order yields two new terms due to the interference between the processes described by ${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{2}$ which are absent in the second order perturbation calculation. The term (23a) was first identified by Kondo [34] and can lead to temperature broadened logarithmic features in the conductance at the energy of intermediate states. For a non-vanishing potential scattering amplitude $(U\ne 0)$ the term (23b) will, in addition, produce a bias-asymmetry in the differential conductance even without spin-polarized electrodes.

We start with the evaluation of the Kondo-like processes, that are described by the term (23a), using a spin $S=1/2$ system with only two states as an example. We assume that in thermal equilibrium only the ground state $| \uparrow \rangle =| {\psi }_{1}\rangle$ is occupied, i.e. that the Zeeman-splitting induced by the external magnetic field is large compared to the thermal energy: $g{\mu }_{B}| B| \gg {k}_{B}T$. The Feynman diagrams of figure 6 depict all possible interaction processes. First, we review the second order processes (figure 6(a)): a spin-up electron that tunnels from tip to the sample cannot flip the spin ($| {\psi }_{1}\rangle \to | {\psi }_{1}\rangle$), while a spin-down electron can either scatter with exchange of angular momentum ($| {\psi }_{1}\rangle \to | {\psi }_{2}\rangle$) leaving the system in the state $| \downarrow \rangle =| {\psi }_{2}\rangle$ or without exchange ($| {\psi }_{1}\rangle \to | {\psi }_{1}\rangle$). To third order, there are a total of six diagrams to be accounted for, which we can label by the states occupied in the initial (i), intermediate (m), and final (f) state. Here, we mark the exchange diagrams by appending a 'R'. This results in the processes (112), (121), (122), and the reversed order scattering events (112R),(121R), (122R).

Zoom In Zoom Out Reset image size

To evaluate the conductance due to these processes we calculate the spin-flip matrix elements of equation (15) for the electrons and the localized spin-system. For electrodes without spin-polarization this results in conductances of:

$$\begin{eqnarray}\begin{array}{ccc}{\displaystyle \frac{\partial I}{\partial V}}{({eV})}^{\rm{t}\to \rm{s}} & = & {\displaystyle \frac{4\pi {e}^{2}}{{\hslash }}}{T}_{0}^{2}{\displaystyle \frac{J{\rho }_{\rm{s}}}{4\rm{i}}}{\displaystyle \sum }_{\underset{\{x,y,z\}}{j,k,l=}}{\varepsilon }_{{jkl}}\mathfrak{R}\left[\left.\langle {\psi }_{i}\right|{\hat{S}}_{l}\left|{\psi }_{f}\right.\rangle\left.\langle {\psi }_{f}\right|{\hat{S}}_{k}\left|{\psi }_{m}\right.\rangle\left.\langle {\psi }_{m}\right|{\hat{S}}_{j}\left|{\psi }_{i}\right.\rangle\right]\\ & & \times F({eV}-{\epsilon }_{{im}},T)\times \Theta ({eV}-{\epsilon }_{{if}},T),\end{array}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}\begin{array}{ccc}{\displaystyle \frac{\partial I}{\partial V}}{({eV})}^{\rm{t}\stackrel{\rm{R}}{\to }\rm{s}} & = & {\displaystyle \frac{4\pi {e}^{2}}{{\hslash }}}{T}_{0}^{2}{\displaystyle \frac{J{\rho }_{\rm{s}}}{4\rm{i}}}{\displaystyle \sum }_{\underset{\{x,y,z\}}{j,k,l=}}{\varepsilon }_{{jkl}}\mathfrak{R}\left[\left.\langle {\psi }_{i}\right|{\hat{S}}_{l}\left|{\psi }_{f}\right.\rangle\left.\langle {\psi }_{f}\right|{\hat{S}}_{k}\left|{\psi }_{m}\right.\rangle\left.\langle {\psi }_{m}\right|{\hat{S}}_{j}\left|{\psi }_{i}\right.\rangle\right]\\ & & \times F({eV}-{\epsilon }_{{mi}},T)\times \Theta ({eV}-{\epsilon }_{{if}},T).\end{array}\end{eqnarray} \tag{ 25 }$$

Here, equation (24) accounts for the direct diagrams (normal order) and (25) for the exchange diagrams (reversed order), respectively. In these equations ${\varepsilon }_{{jkl}}$ is the usual Levi-Civita tensor of rank three, which is 1 (−1) if $\{{jkl}\}$ is an even (odd) permutation of $\{{xyz}\}$, and zero otherwise [77]. The step-function Θ ensures that energy conservation at the final state of the scattering process is obeyed, while the function F has its peak at the intermediate state energy and is mirrored at E_F for the reversed tunneling process. Any spin-polarization in tip or sample changes the transition matrix elements for the interacting electrons and equations (24) and (25) would become more complex (see appendix).

Comparing the different contributions to the conductance (figures 6(b) and 7(a)) it is remarkable that all third order scattering events start with a spin-down electron except the process 121R. Here, a spin-up electron tunneling from tip to sample non-trivially interacts with the system because a substrate electron has flipped the localized spin into the intermediate spin-down state before the interaction takes place.

Summing up the contributions to the conductance due to all second and third order scattering events results in spectra similar to those exemplary given in figure 7(b) (cf figures 1(e) and (f) [22]). At zero field the differential conductance shows a peak that splits with increasing magnetic field. While in this calculation we assume to be in the weak-coupling Kondo limit and thus neglect any correlation energy due to the formation of a Kondo singlet, the peak splits as soon as the Zeeman energy overcomes the thermal energy. Note, that the resulting split-peak at small fields can lead to the deduction of an erroneously high g-factor when just evaluating the peak positions due to the superposition of peak and step-structures.

Zoom In Zoom Out Reset image size

If we consider now in addition that the scattering process between tip and sample or vice versa can also occur via the potential interaction (term (23b)), we observe an asymmetric line-shape as soon as the two eigenstates are no longer degenerate and $\langle \hat{{S}_{z}}\rangle \ne 0$ (figures 7(c) and (d)). This direction-dependent asymmetry cannot originate from the scattering matrix elements that involve only the localized spin site (the order of excitations shall be the same) but derives from the matrix elements involving the interacting electrons. Physically, the reason lies in the asymmetry of the tunnel-junction for which we assume that only the sample is coupled to the spin-system and thus neglect any intermediate scattering process that originate and end in the tip.

As an example, we examine in detail the process (121) (see figure 6(b)) in both tunneling directions. For a current to flow from tip to sample via this process, first a $| \downarrow \rangle$ electron originating from the tip is scattered into a $| \uparrow \rangle$ state in the sample exciting the spin system from $| {\psi }_{1}\rangle \to | {\psi }_{2}\rangle$. Second, a $| \uparrow \rangle$ electron in the sample is scattered into $| \downarrow \rangle$ bringing the spin system back to $| {\psi }_{1}\rangle$. Third, the sample $| \downarrow \rangle$ is scattered back into the tip as $| \downarrow \rangle$ electron. The last step of this process can either take place via the ${\hat{\sigma }}_{z}$ or, in the case of potential scattering, via the ${\hat{\sigma }}_{I}$ term, respectively. For an electron tunneling in the reverse direction we first have a $| \downarrow \rangle$ from the sample being scattered into an electron $| \uparrow \rangle$ state of the tip simultaneously exciting the spin system from $| {\psi }_{1}\rangle \to | {\psi }_{2}\rangle$. Then, the second process flips a $| \uparrow \rangle$ sample electron into the $| \downarrow \rangle$ hole from the first process. Finally, a $| \uparrow \rangle$ tip electron traverses the junction and fills the $| \uparrow \rangle$ hole in the sample. Comparing both tunneling directions we see that it is the third process which differs in the initial and final state, i.e. the matrix elements are either $\langle \downarrow | {\hat{\sigma }}_{z}| \downarrow \rangle$ and $\langle \downarrow | {\hat{\sigma }}_{I}| \downarrow \rangle$ or $\langle \uparrow | {\hat{\sigma }}_{z}| \uparrow \rangle$ and $\langle \uparrow | {\hat{\sigma }}_{I}| \uparrow \rangle$. Rearranging the matrix elements so that all processes become electron-like, the prefactors for calculating the tunneling conductance due to these four processes are for the two without potential scattering:

$$\begin{eqnarray} &&{(121)}^{\rm{t}\to \rm{s}}:+J{\rho }_{\rm{s}}\mathop{\underbrace{\langle {\downarrow }^{\rm{t}}| {\hat{\sigma }}_{z}| {\downarrow }^{\rm{s}}\rangle }}\limits_{-}\mathop{\underbrace{\langle {\downarrow }^{\rm{s}}| {\hat{\sigma }}_{-}| {\uparrow }^{\rm{s}}\rangle }}\limits_{+}\mathop{\underbrace{\langle {\uparrow }^{\rm{s}}| {\hat{\sigma }}_{+}| {\downarrow }^{\rm{t}}\rangle }}\limits_{+}\;\gt 0,\\ &&{(121)}^{\rm{s}\to \rm{t}}:-J{\rho }_{\rm{s}}\mathop{\underbrace{\langle {\uparrow }^{\rm{t}}| {\hat{\sigma }}_{+}| {\downarrow }^{\rm{s}}\rangle }}\limits_{+}\mathop{\underbrace{\langle {\downarrow }^{\rm{s}}| {\hat{\sigma }}_{-}| {\uparrow }^{\rm{s}}\rangle }}\limits_{+}\mathop{\underbrace{\langle {\uparrow }^{\rm{s}}| {\hat{\sigma }}_{z}| {\uparrow }^{\rm{t}}\rangle }}\limits_{+}\;\gt 0, \end{eqnarray} \tag{ 26 }$$

and with potential scattering:

$$\begin{eqnarray} &&{(121u)}^{\rm{t}\to \rm{s}}:+J{\rho }_{\rm{s}}U\mathop{\underbrace{\langle {\downarrow }^{\rm{t}}| {\hat{\sigma }}_{I}| {\downarrow }^{\rm{s}}\rangle }}\limits_{+}\mathop{\underbrace{\langle {\downarrow }^{\rm{s}}| {\hat{\sigma }}_{-}| {\uparrow }^{\rm{s}}\rangle }}\limits_{+}\mathop{\underbrace{\langle {\uparrow }^{\rm{s}}| {\hat{\sigma }}_{+}| {\downarrow }^{\rm{t}}\rangle }}\limits_{+}\;\lt 0,\\ &&{(121u)}^{\rm{s}\to \rm{t}}:-J{\rho }_{\rm{s}}U\mathop{\underbrace{\langle {\uparrow }^{\rm{t}}| {\hat{\sigma }}_{+}| {\downarrow }^{\rm{s}}\rangle }}\limits_{+}\mathop{\underbrace{\langle {\downarrow }^{\rm{s}}| {\hat{\sigma }}_{-}| {\uparrow }^{\rm{s}}\rangle }}\limits_{+}\mathop{\underbrace{\langle {\uparrow }^{\rm{s}}| {\hat{\sigma }}_{I}| {\uparrow }^{\rm{t}}\rangle }}\limits_{+}\;\gt 0, \end{eqnarray} \tag{ 27 }$$

where we assumed $J{\rho }_{\rm{s}}\lt 0$ and $U\gt 0$. Note, that the preceding sign change at the tunneling direction from sample to tip ($\rm{s}\to \rm{t}$) is due to the rearrangement of the interaction order together with the switching from hole-like to electron-like scattering. The contribution to the conductance from the processes in which only Kondo-like spin–spin interactions take place is positive for both tunneling directions, while the conductance for processes that include potential scattering changes its sign when inverting the tunneling direction.

Up to here we have revised the relevant formalism to calculate the conductance in the third order of the matrix elements quite closely following the work established by Appelbaum, Anderson, and Kondo [34, 95–97]. The $S=1/2$ system we used for illustration contained only two states and was not influenced by any magnetic anisotropy or near-by spin systems. Recently, efforts have been made to expand this perturbative model to higher spin systems, which also include magnetic anisotropy and couplings to neighboring spins [81, 83–85], but the importance of the potential scattering was not taken into account up to now. In the following, the power of this easily accessible model will be used to evaluate and describe the spectral features on more complex systems.

3.1. The smallest ($S\;=$ 1) high-spin system

Regarding the spectra of the $S=1/2$ system (figure 7), one could get the impression that the coupling to the sample always results in some superimposed peak-like structures as soon as a spin flip transition is possible, and which scales with the substrate coupling strength $-J{\rho }_{\rm{s}}$. In this section we will show that the situation is more complex and depend not only on the transition matrix elements but also on the state energies.

To expand the complexity, we turn now to a magnetic system with S = 1 in which axial and transverse magnetic anisotropy have broken the zero-field degeneracy of the three eigenstates. Assuming easy-axis anisotropy ($D\lt 0$) the energetically lowest eigenstates have weights only for ${m}_{z}=\pm 1$. The additional transverse anisotropy splits the remaining two-fold degeneracy forming an antisymmetric ground and a symmetric first excited state (figure 8(a)). Thus, in the basis of the magnetic easy axis the three eigenstates can be written as: $| {\psi }_{1}\rangle ={\displaystyle \frac{1}{\sqrt{2}}}(| +1\rangle -| -1\rangle )$, $| {\psi }_{2}\rangle ={\displaystyle \frac{1}{\sqrt{2}}}(| +1\rangle +| -1\rangle )$, and $| {\psi }_{3}\rangle =| 0\rangle$. Such situations have been found for example for Fe(II)-phthalocyanine molecules adsorbed on the (2 × 1) oxygen terminated Cu(110) surface [65], for Fe atoms adsorbed on InSb(110) [53], or for CoH complexes adsorbed on the h-BN/Rh(111) surface [51].

Zoom In Zoom Out Reset image size

By taking only second order spin-flip processes into account it is easy to calculate (using equation (17)) that the transition probabilities from the groundstate to the two excited states at B = 0 are equal, leading to a symmetric double step structure in the spectrum (figure 8). Expanding the calculation to third order, only the processes which involve all states, i.e. (123) and (132), have non-zero matrix elements when the potential scattering U = 0. Processes like (112), (121), or (131) cannot be interlinked using any combination of the operators ${\hat{S}}_{+}$, ${\hat{S}}_{-}$, and ${\hat{S}}_{z}$ and are thus forbidden. Nevertheless, additional potential scattering enables transitions involving the operators ${\hat{S}}_{+}$, ${\hat{S}}_{-}$, and $\hat{I}$ or ${\hat{S}}_{z}$, ${\hat{S}}_{z}$, and $\hat{I}$ which makes the aforementioned processes possible resulting in an asymmetry of the spectrum with respect to tunneling direction.

It is remarkable that the remaining third order processes at U = 0 change the spectrum in a quite different fashion, even though they have the same strength. Process (123) produces its peaks at ${eV}=\pm {\epsilon }_{2}$, but due to energy conservation it contributes to the spectrum only at $| {eV}| \gt {\epsilon }_{3}$, which efficiently cuts off the peak. In contrast, the peak at higher energy due to process (132) is not as strongly cut off (figure 8(b)). Thus, the full spectrum shows a peak-like increase of the conductance at the energy of the second step but not at the energetically lower first transition step. Furthermore, the conductance has a curved form for tunneling voltages between ${\epsilon }_{2}$ and ${\epsilon }_{3}$ (figure 8(c)). Note, the overall general behavior would not alter when changing from easy-axis to easy-plane anisotropy ($D\gt 0$) as long as all degeneracies are lifted.

3.2. Single Mn and Fe atoms on Cu₂N

After the, with only three eigenstates, rather easy S = 1 example, we now apply the model to the experimentally and theoretically intensively studied single 3d transition metal atoms adsorbed on a monolayer of Cu₂N on Cu(100). When these atoms are placed on top of a Cu site they form strong covalent bonds with the neighboring N atoms [8]. This highly anisotropic adsorption geometry leads to three distinct symmetry axes that are perpendicular to each other: The direction out-of-surface and two in-surface directions along the Cu-N bonds and perpendicular to it, along the so called vacancy rows (figure 9(a) inset).

Zoom In Zoom Out Reset image size

Single Fe atoms adsorbed on this surface have been found to be in the S = 2 state with a magnetic easy-axis along the N rows (z-direction) and a magnetically hard-axis along the vacancy row (x-direction). In this coordinate system anisotropy values of $D=-1.55$ mV and E = 0.31 mV, and a gyromagnetic factor of g = 2.11 described the experimental data well using a spin Hamiltonian like equation (5) and a second order tunneling model [8, 75, 102]. The Hamiltonian has as solution five non-degenerate eigenstates and, due to $D\lt 0$, favors, at zero field, ground states with weights at high m_z values. Similar to the S = 1 system the transverse anisotropy breaks the degeneracies leading to a symmetric and antisymmetric solution with the main weights at $| \pm 2\rangle$ as ground and first excited state and weights in $| \pm 1\rangle$ for the second and third excited state (figure 9(b)). To visualize the anisotropy we plot the total energy necessary to rotate the ground state into arbitrary directions showing the favored easy axis (z) and the unfavored hard axis (x) (figure 9(b)) [89]. In second order, spin-flip scattering is allowed between the groundstate and the three lowest excited states but a transition to the highest state is forbidden because this would require an exchange of $\Delta m=\pm 2$.

Experimental $\rm{d}I/\rm{d}V$ measurements on this system show, in addition to the conductance steps, peak-like structures at the second and third step but not at the lowest one (figure 9(a)). Additionally, they show an asymmetry between positive and negative bias. To rationalize these observations we can follow a similar argumentation as in the S = 1 case: in third order, transitions like (121) are not possible without additional potential interaction and processes like (123) or (124) are strongly cut off due to the high energy difference between ${\epsilon }_{2}$ and ${\epsilon }_{3}$ or ${\epsilon }_{4}$. In contrast, the processes (132) and (142) scale with $J{\rho }_{\rm{s}}$ leading to the peak features in the differential conductance. The additional asymmetry hints at a non-negligible potential scattering. As the computed curves in figure 9(a) reveal, our model almost perfectly fits the magnetic field data without any adaption of the parameters³ . We mention that similar good agreement between experimental data and computation can be reached in the other two magnetic field directions (not shown). The coupling strength in these simulations is $J{\rho }_{\rm{s}}=-0.087$, close to the −0.1 found in a similar perturbative approach [84]. A potential scattering term of U = 0.35 is necessary to reproduce the asymmetry. This value is significantly smaller than the $U\approx 0.75$ found in experiments where the magneto-resistive elastic tunneling was probed [99]. Part of this discrepancy can be understood by an additional conductance term that does not coherently interact with the spin-system and which would lead to an overestimation of U in magneto-resistive measurements. Indeed we need a constant conductance offset of about 20%, which is added to the calculated conductance to reproduce the spectra.

Switching from an integer to a half-integer spin system we now discuss individual Mn atoms on Cu₂N, which have a spin of $S=5/2$ and only a small easy-axis anisotropy of $D\approx -40\ \mu \mathrm{eV}$ along the out-of-plane direction and a negligible transverse anisotropy [7, 8]. The easy-axis anisotropy prohibits the immediate formation of a Kondo state due to a Kramer's degenerate ground state doublet with ${m}_{z}=\pm 5/2$ (figure 9(c)). At zero field a typical spectrum shows only one step, which belongs to the transition between the $\pm 5/2$ and the $\pm 3/2$ states that have superimposed asymmetric peak structures (figure 9(d)). The fit to the model yields $J{\rho }_{\rm{s}}=-0.029$ and resembles a $S=1/2$ split-Kondo peak at small magnetic fields (see figure 7(b)). A different Mn atom investigated at B_z = 7 T shows a significantly reduced $J{\rho }_{\rm{s}}=-0.0091$. Interestingly, we find for both atoms a potential scattering value of $U\approx {\displaystyle \frac{1}{2}}S$, which allows one to describe the spectra without the need of any additional conductance offset. This high U value that is the origin of the bias asymmetry has been independently found in spin-pumping experiments [13] (see section 5) and by measuring the magneto-resistive elastic tunneling contribution [99]. The extraordinary agreement between model and experiment can be seen in measurements using different spin-polarized tips on the same atom (figure 9(e)). Here, the strong influence of the tip-polarization on the inelastic conductance at bias voltages $| V| \lt 1$ mV is evident while the differential conductance at $V\gt 1.5$ mV stays, for all tips, constant.

When a half-integer spin with $S\gt 1/2$ has easy-plane anisotropy $D\gt 0$, its ground state at zero field is a doublet with its main weights in ${m}_{z}=\pm 1/2$. This enables an effective scattering with the substrate electrons and leads at low enough temperatures to the formation of a Kondo state. Experimentally this has been observed for Co atoms on Cu₂N [9, 11, 16–18] which have been found to have the total spin $S=3/2$ and which enter the correlated Kondo state with a characteristic Kondo temperature of T_K = 2.6 K in experiments performed on small patches of Cu₂N at temperatures down to T = 550 mK [9]. Apart from $D\gt 0$ the system also has a small in-plane anisotropy ($E\ne 0$) which creates an easy axis (x) along the nitrogen row and a hard axis (z) along the vacancy rows (figure 10(a)).

Zoom In Zoom Out Reset image size

Interestingly, in this system the coupling to the substrate $J{\rho }_{\rm{s}}$ changes with the position of the Co atom on larger Cu₂N patches, concomitant with a change in the anisotropy energies which separates the $| \pm 1/2\rangle$ states from the $| \pm 3/2\rangle$ ones [16, 94]. For us this allows the study of the transition from the weak coupling to the strong coupling regime. In the case where the Co atom is relatively weakly coupled to the substrate ($J{\rho }_{\rm{s}}\approx -0.1$), the model can be consistently fitted to the experimental data even at different field strengths along B_x (figure 10(b)). We observe a zero-energy peak that splits at applied fields similar to the $S=1/2$ system of figure 7(d). But while for $S=1/2$ the field direction does not play a role, here the peak splitting depends strongly on the direction due to the magnetic anisotropy [9]. For this high-spin system we furthermore observe inelastic steps due to the transition to energetically higher excited states which are located at $| {eV}| =2| D|$ for B = 0 and whose additional peak structure is well described within the model.

For Co atoms adsorbed closer to the edges of the Cu₂N patches, the coupling to the substrate increases and the fit to the model worsens significantly (figure 10(c)). While the experimental data measured for non-zero fields are reasonable well described with $J{\rho }_{\rm{s}}\approx -0.25$, at B = 0 the experimentally detected peak at E_F is stronger than the peak created by the model. Additionally, the experimental peak-width appears to be smaller than the temperature broadened logarithmic function, which was similarly observed for a radical molecule with $S=1/2$ for temperatures presumably close to T_K [22]. Furthermore, we observe that the calculated spectrum no longer well describes the steps at $2| D|$. At this energy the third order contributions are less pronounced in the experimental data, indicating that we reach the limit of the perturbative approach.

The full description of a spin system in the strong coupling regime requires complex theoretical methods like numerical renormalization group theory [36, 103, 104] which are beyond the scope of this paper. Nevertheless, we can discuss some of the physical consequences within the outlined perturbative model. In contrast to the examples discussed in section 3, the two lowest degenerate ground states of the Co/Cu₂N system have weights in states that are separated by $\Delta m=\pm 1$. Thus, at zero field, electrons from the substrate can efficiently flip between these two states. The computation of the transition rate between the two states $| {\psi }_{i}\rangle$ and $| {\psi }_{f}\rangle$ of the spin system that have the energies ${\epsilon }_{i}$ and ${\epsilon }_{f}$ is ${\Gamma }_{{if}}:= {I}_{{if}}/e$, i.e. the transition probability per time, is analogous to the calculation of the current (equation (9)) given by:

$$\begin{eqnarray}&&{\Gamma }_{{if}}^{\rm{s}\to \rm{s}}={\displaystyle \frac{{I}_{{if}}^{\rm{s}\to \rm{s}}}{e}}={\displaystyle \frac{2\pi }{{\hslash }}}{(J{\rho }_{\rm{s}})}^{2}{{\displaystyle \int }}_{-\infty }^{\infty }\rm{d}\epsilon \;{| {\mathcal{M}}_{{if}}^{\rm{s}\to \rm{s}}| }^{2}f(\epsilon ,T)[1-f(\epsilon -{\epsilon }_{{if}},T)].\end{eqnarray} \tag{ 28 }$$

Here we have set the bias to eV = 0 and replaced the coupling constant T₀ between adsorbate and tip with $J{\rho }_{\rm{s}}$, the coupling between adsorbate and substrate. When we now consider only processes up to second order, the matrix ${\mathcal{M}}_{{if}}$ is independent of the energy and the integral can be evaluated to

$$\begin{eqnarray}&&{{\displaystyle \int }}_{-\infty }^{\infty }\rm{d}\epsilon \;f(\epsilon ,T)[1-f(\epsilon -{\epsilon }_{{if}},T)]={\displaystyle \frac{{\epsilon }_{{if}}}{\mathrm{exp}\left({\displaystyle \frac{{\epsilon }_{{if}}}{{k}_{B}T}}\right)-1}}.\end{eqnarray} \tag{ 29 }$$

Thus, the scattering rate between the two degenerate states (i.e. ${\epsilon }_{12}=0$) is directly proportional to the temperature (figure 11(a)):

$$\begin{eqnarray}&&{\Gamma }_{12}^{(2)\rm{s}\to \rm{s}}={\displaystyle \frac{2\pi }{{\hslash }}}{(J{\rho }_{\rm{s}})}^{2}{| {\tilde{M}}_{12}| }^{2}{k}_{B}T.\end{eqnarray} \tag{ 30 }$$

Zoom In Zoom Out Reset image size

These scattering processes, which we have discussed in section 2 and displayed in figure 6(a), will tend to change the spin polarization of the electronic states in the sample near the adsorbate to be correlated with the localized spin. Nevertheless, this local correlation will be quickly destroyed by decoherent scattering processes with the remaining electron bath, which we can assume to be large and dissipative. This decoherence rate is also proportional to the temperature, ${\Gamma }_{\mathrm{decoh}}\propto {k}_{B}T$ [105], but usually stronger so that no highly correlated state can form. This picture changes when we additionally consider third order scattering processes which yield, using equation (23), the probability:

$$\begin{eqnarray}\begin{array}{ccc}{\Gamma }_{{if}}^{(3)\rm{s}\to \rm{s}} & = & {\displaystyle \frac{4\pi }{{\hslash }}}{(J{\rho }_{\rm{s}})}^{3}{{\displaystyle \int }}_{-\infty }^{\infty }\rm{d}\epsilon \;{\displaystyle \sum }_{m}\mathfrak{R}({\tilde{M}}_{{fi}}{\tilde{M}}_{{mf}}{\tilde{M}}_{{im}})\\ & & \times \left[F({\epsilon }_{{mi}}-\epsilon ,T)+F({\epsilon }_{{im}}-\epsilon ,T)\right]f(\epsilon ,T)[1-f(\epsilon -{\epsilon }_{{if}},T)].\end{array}\end{eqnarray} \tag{ 31 }$$

Due to the growing intensity of $F(\epsilon =0)$ at reduced temperatures (figure 5), for temperatures $T\to 0$ the scattering ${\Gamma }^{(3)}$ decreases significantly more slowly than ${\Gamma }^{(2)}$ (see figure 11(b)) so that their ratio steadily increases:

$$\begin{eqnarray}&&{\displaystyle \frac{{\Gamma }^{(3)}}{{\Gamma }^{(2)}}}\approx J{\rho }_{\rm{s}}\mathrm{ln}\left({\displaystyle \frac{{k}_{B}T}{{\omega }_{0}}}\right).\end{eqnarray} \tag{ 32 }$$

In contrast to that, the decoherent scattering rate with the bath ${\Gamma }_{\mathrm{decoh}}$ lacks localized scattering centers and therefore has no significant third order contributions. Equation (32) leads to a characteristic temperature, the Kondo temperature T_K, where ${\Gamma }^{(3)}$ and higher order scattering terms become the dominant processes and perturbation theory breaks down [34, 36, 106, 107]:

$$\begin{eqnarray}&&{T}_{K}\approx {\displaystyle \frac{{\omega }_{0}}{{k}_{B}}}\mathrm{exp}\left({\displaystyle \frac{1}{J{\rho }_{\rm{s}}}}\right).\end{eqnarray} \tag{ 33 }$$

Below this temperature the assumption used up to now, i.e. that the electronic states in the sample are not influenced by the presence of the spin system, no longer applies. Using exact methods like the modified Bethe ansatz [108, 109] or numerical renormalization group theory [104] reveals that the sample electrons in the vicinity rather form an entangled state with the impurity, i.e.

$$\begin{eqnarray}&&{\Psi }^{\mathrm{total}}={\displaystyle \frac{1}{\sqrt{2}}}\left(\left|{\downarrow }^{\rm{s}}\right.\rangle\left|+{\displaystyle \frac{1}{2}}\right.\rangle-\left|{\uparrow }^{\rm{s}}\right.\rangle\left|-{\displaystyle \frac{1}{2}}\right.\rangle\right),\end{eqnarray} \tag{ 34 }$$

as illustrated in figure 11(c). This combined state is quite complicated because the electronic states in the sample are continuous in energy and extend spatially, and therefore strongly alter the excitation spectrum of the adsorbate [110, 111].

Figure 11(d) shows the impact of the formation of the correlated state on the experimentally detected spectrum of the Co/Cu₂N system. At temperatures $T\ll {T}_{K}$, the peak at the Fermi energy can no longer be well reproduced by a temperature broadened logarithmic function (which in any case must diverge for $T\to 0$). The peak is much better described by an asymmetric Lorentzian, i.e. a Fano line-shape [47], or the Frota function [101, 112, 113] which has a finite amplitude and a half-width of $\Delta ={k}_{B}{T}_{K}$, corresponding to the correlation energy of the Kondo state. Superimposed on this zero-energy peak we detect conductance steps at voltages that enable scattering of the spin system into the ${m}_{z}=\pm 3/2$ state. These steps are well described using only a second order perturbation spin-flip model [9] and do not show any third order logarithmic contributions. This means that the probability of the third order scattering channels must be closed due to the ground-state correlation between the localized spin and the substrate electrons. Here we see that the simple perturbative approach of the model is no longer valid and fails to capture the full physics in this strongly coupled Kondo regime. Nevertheless, the model still gives valuable information because it enables us to identify the conditions under which correlations can evolve and is due to its computational simplicity also applicable for larger and more complex spin systems where the calculation of the exact solution might not be feasible or very time consuming.

In the preceding section we discussed the appearance of correlations due to higher order scattering between the substrate electrons and the localized spin system. Now we want to evaluate how the tunneling of electrons between the biased tip and the sample influences the state populations and the observable spectroscopic features in local conductance measurements. This is equivalent to dropping the zero-current approximation we have assumed until now. We will see that the change of time-averaged system properties results in characteristic fingerprints in the $\rm{d}I/\rm{d}V$ signal which can become crucial for a deeper understanding of the system's response and dynamics.

To describe the transition rate between the localized spin states $| {\psi }_{i}\rangle$ and $| {\psi }_{f}\rangle$ due to the tunneling current flowing between tip and sample we can again use equation (9) and the relation ${\Gamma }_{{if}}^{\rm{t}\to \rm{s}}({eV}):= {I}_{{if}}^{\rm{t}\to \rm{s}}({eV})/e$, which leads, in first order Born approximation, to:

$$\begin{eqnarray}&&{\Gamma }_{{if}}^{(2)\rm{t}\to \rm{s}}({eV})={\displaystyle \frac{2\pi }{{\hslash }}}{T}_{0}^{2}{{\displaystyle \int }}_{-\infty }^{\infty }\rm{d}\epsilon \;{| {\mathcal{M}}_{{if}}^{\rm{t}\to \rm{s}}| }^{2}f(\epsilon -{eV},T)[1-f(\epsilon -{\epsilon }_{{if}},T)].\end{eqnarray} \tag{ 35 }$$

The integral can be solved using equation (29) resulting in:

$$\begin{eqnarray}&&{\Gamma }_{{if}}^{(2)\rm{t}\to \rm{s}}({eV})={\displaystyle \frac{2\pi }{{\hslash }}}{T}_{0}^{2}{| {\mathcal{M}}_{{if}}^{\rm{t}\to \rm{s}}| }^{2}{\displaystyle \frac{{\epsilon }_{{if}}-{eV}}{\mathrm{exp}\left({\displaystyle \frac{{\epsilon }_{{if}}-{eV}}{{k}_{B}T}}\right)-1}}.\end{eqnarray} \tag{ 36 }$$

As long as ${eV}-{\epsilon }_{{if}}\gg {k}_{B}T$, equation (36) can be further simplified to an equation that is linear in the energy difference:

$$\begin{eqnarray}&&{\Gamma }_{{if}}^{(2)\rm{t}\to \rm{s}}({eV})=(2\pi /{\hslash }){T}_{0}^{2}{| {\mathcal{M}}_{{if}}^{\rm{t}\to \rm{s}}| }^{2}({eV}-{\epsilon }_{{if}}).\end{eqnarray} \tag{ 37 }$$

This linear dependence of the scattering rate on the energy difference is equivalent to the assumptions that in second order perturbation the matrix elements and the coupling constant T₀ between tip and sample are energy independent. Under these assumptions the rate is given by the energy window between the available energetically hot electrons and the energy needed to change the localized state [79]. However, one should keep in mind that at large energy differences the intrinsic variation of the local density of electronic states in tip and sample $\rho {(\epsilon )}^{\rm{t},\ \rm{s}}$ might alter the results.

Furthermore, we additionally include third order contributions to the tunneling transition rates. This can be archived by integrating the interference terms between ${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{2}$ of equation (23) and leads to an additional contribution to the scattering rate of:

$$\begin{eqnarray}\begin{array}{ccc}{\Gamma }_{{if}}^{(3)\rm{t}\to \rm{s}}({eV}) & = & {\displaystyle \frac{4\pi }{{\hslash }}}{T}_{0}^{2}J{\rho }_{\rm{s}}{{\displaystyle \int }}_{-\infty }^{\infty }\rm{d}\epsilon \;{\displaystyle \sum }_{m}f(\epsilon -{eV},T)[1-f(\epsilon -{\epsilon }_{{if}},T)]\\ & & \times \left[\mathfrak{R}({M}_{{fi}}{\tilde{M}}_{{mf}}{M}_{{im}})F({\epsilon }_{{mi}}-\epsilon ,T)\right.\\ & & +\mathfrak{R}({M}_{{fi}}{M}_{{mf}}{\tilde{M}}_{{im}})F({\epsilon }_{{im}}-\epsilon ,T)\\ & & +U\mathfrak{R}({I}_{{fi}}{\tilde{M}}_{{mf}}{M}_{{im}})F({\epsilon }_{{mi}}-\epsilon ,T)\\ & & +\left.U\mathfrak{R}({I}_{{fi}}{M}_{{mf}}{\tilde{M}}_{{im}})F({\epsilon }_{{im}}-\epsilon ,T)\right].\end{array}\end{eqnarray} \tag{ 38 }$$

For tunneling currents flowing in the opposite direction, i.e. for electrons which are scattered at the spin system when tunneling from sample to tip, equations (36) and (38) can be adapted straightforwardly.

The total transition rates ${\Gamma }_{{ij}}^{\rm{t}\to \rm{s}}={\Gamma }_{{ij}}^{(2)\rm{t}\to \rm{s}}+{\Gamma }_{{ij}}^{(3)\rm{t}\to \rm{s}}$ and ${\Gamma }_{{ij}}^{\rm{s}\to \rm{t}}={\Gamma }_{{ij}}^{(2)\rm{s}\to \rm{t}}+{\Gamma }_{{ij}}^{(3)\rm{s}\to \rm{t}}$ have to be evaluated for all possible initial and final localized spin states and, together with the dissipative substrate to substrate scattering rates discussed in section 4 (equations (28) and (31)), can be brought together to the characteristic master equation for the state populations:

$$\begin{eqnarray}\begin{array}{ccc}{\displaystyle \frac{\rm{d}}{\rm{d}t}}{p}_{i}(t) & = & {\displaystyle \sum }_{j\ne i}{p}_{j}(t)({\Gamma }_{{ji}}^{\rm{t}\to \rm{s}}+{\Gamma }_{{ji}}^{\rm{s}\to \rm{t}}+{\Gamma }_{{ji}}^{\rm{s}\to \rm{s}})-{p}_{i}(t){\displaystyle \sum }_{j\ne i}({\Gamma }_{{ij}}^{\rm{t}\to \rm{s}}+{\Gamma }_{{ij}}^{\rm{s}\to \rm{t}}+{\Gamma }_{{ij}}^{\rm{s}\to \rm{s}})\\ & = & {{\Upsilon }}_{{ij}}(t){p}_{j}(t).\end{array}\end{eqnarray} \tag{ 39 }$$

This set of first order differential equations describes the time-development of an initial state of the spin system, which might be a superposition state $| \psi (t)\rangle ={{\displaystyle \sum }}_{i}{p}_{i}(t)\left|{\psi }_{i}\right.\rangle$, at the time t under the influence of all possible, bias dependent, scattering events. The first summation in equation (39) accounts for all the probabilities to scatter into the ith eigenstate from any other eigenstate. These probabilities have to be weighted with the time-dependent state probability p_j of the originating state. The second summation accounts for all scattering events which reduces the ith state probability by scattering it into other eigenstates. It is important to note, that we use three crucial premisses in this approach:

• (i)  
  We assume that the future of the localized quantum state depends only on the actual state and not on any interactions which happened in the past. This Markovian limit means that we assume that the timescales at which this approximation breaks down, in particular where higher than two-time correlations are dominant, can be neglected [114].
• (ii)  
  We do not account for any phase coherences between the eigenstates of the spin system. This means that the phase coherence time is short compared to the state's lifetime τ and the time between successive tunneling events. Under these assumptions all off-diagonal elements of the density matrix in the rotating frame with the eigenstates $\psi \prime$ are zero and thus the density matrix can be described as a state mixture $\chi ={{\displaystyle \sum }}_{i}{p}_{i}| {\psi }_{i}^{\prime }\rangle \langle {\psi }_{i}^{\prime }|$.
• (iii)  
  Because we treat the total system as the product state between the continuous electronic states in tip and sample and the discrete spin states, a correlated state between substrate electrons and the localized spin as discussed in section 4 (equation (34)) can not be developed directly within the limitations of our model.

In the following, we further want to limit our evaluation to steady-state conditions, which means that the change of any external parameter, like the tunneling bias voltage, occurs adiabatically slowly, much more slowly than any relaxation times in the system. The steady-state condition is reached, when all p_i are static, i.e. ${\displaystyle \frac{\rm{d}}{\rm{d}t}}{p}_{i}(t\to \infty )=0$. This is equivalent to finding of the algebraic kernel of the rate matrix Υ, which additionally has to be normalized to account for the conservation of probabilities ${{\displaystyle \sum }}_{i}{p}_{i}=1$:

$$\begin{eqnarray}&&{p}^{\mathrm{stat}}({eV})={\displaystyle \frac{\mathrm{ker}{\Upsilon }({eV})}{\parallel \mathrm{ker}{\Upsilon }({eV})\parallel }}.\end{eqnarray} \tag{ 40 }$$

From these ${p}^{\mathrm{stat}}({eV})$ the current I(eV) and the differential conductance $\rm{d}I({eV})/\rm{d}({eV})$ can be calculated:

$$\begin{eqnarray}&&I({eV})=e{\displaystyle \sum }_{i,j}{p}_{i}^{\mathrm{stat}}({eV})\left({\Gamma }_{{ij}}^{\rm{t}\to \rm{s}}({eV})-{\Gamma }_{{ij}}^{\rm{s}\to \rm{t}}({eV})\right).\end{eqnarray} \tag{ 41 }$$

To show the influence of a non-zero tip–sample coupling strength T₀ on the tunneling spectra we are returning to the spin S = 1 example from section 3.1. In figures 12(a) and (c) we show the development of the $\rm{d}I/\rm{d}V$ spectra at B = 0 and B_z = 10 T for different coupling strengths between sample and tip, respectively. The coupling constants correspond to a tunneling resistance between ${R}_{T}\approx 130$ and $\approx 2.6$ MΩ, which is equivalent to a tunneling current of $I\approx 150$ pA–8 nA at a bias of V = 20 mV. Note, that these are quite typical parameters for STM experiments.

Zoom In Zoom Out Reset image size

Without magnetic field and at small couplings we observe a spectrum which does not differ significantly from the one calculated in the zero-current approximation. The average occupation of the two excited states increases only slowly for voltages above the threshold for these transitions and reaches at $V=\pm 10$ mV only about 4% for p₂ and less than $0.4\%$ for p₃ (dashed lines in figure 12(b)). This situation changes when the tip–sample coupling is increased. At larger T₀ an additional feature appears in the spectrum at about ±4 mV, which is due to transitions from the state $| {\psi }_{2}\rangle$ to $| {\psi }_{3}\rangle$. These transitions are only possible because the probabilities p_i are driven far from thermal equilibrium and p₂ has already a significant weight at $| V| =4$ mV. Similar current induced pumping to higher excitation states has been observed for example for Mn-dimers adsorbed on Cu₂N [13], for small Fe clusters containing only a few atoms on Cu(111) [2], or for Fe-OEP-Cl (Fe-octaethylporphyrin-chloride) molecules adsorbed on Pb(111) [115]. Interestingly, the latter experiment was performed on a superconducting surface and with a superconducting Pb-tip which made it necessary to account for the gap around the Fermi energy in the quasi-particle density of states of tip and sample.

If we now add a magnetic field of B_z = 10 T to our simulation, the situation changes: the energy differences ${\epsilon }_{12}$ and ${\epsilon }_{23}$ between the first and second, and second and third eigenstate, respectively, are almost identical (see figure 12(c) inset) masking the additional step-like feature when pumping into energetically higher states. Nevertheless, also here a clear change of the spectrum occurs at high tip–sample couplings. We observe an increased relative conductance between ${\epsilon }_{2}$ and ${\epsilon }_{3}$ and a stronger third order peak at ${eV}=\pm {\epsilon }_{3}$. This feature is mainly due to the third order scattering process $(232)$ now possible, as apparently visible when taking the relative difference between the spectra at strongest and weakest coupling strength (dashed dotted line in figure 12(c)).

As long as the tip is not spin-polarized, the current induced pumping into higher states cannot lead to an inversion of the state occupancy. The state's lifetime ${\tau }_{i}$ for the eigenstate $| {\psi }_{i}\rangle$ is inversely proportional to the sum of all scattering processes which leaves the state:

$$\begin{eqnarray}&&{\tau }_{i}^{-1}={\displaystyle \sum }_{j\ne i}\left({\Gamma }_{{ij}}^{\rm{t}\to \rm{s}}({eV})+{\Gamma }_{{ij}}^{\rm{s}\to \rm{t}}({eV})+{\Gamma }_{{ij}}^{\rm{s}\to ,\rm{s}}({eV})\right).\end{eqnarray} \tag{ 42 }$$

While for spin-averaged electrodes the scattering matrix elements do not change when inverting the initial and final state, the energetically higher states must have always a shorter lifetime.

This behavior changes drastically when a spin-polarized tip is used. We have seen that, in particular for spin systems with a potential scattering $U=S/2$, the tunneling conductance and thereby the scattering rate depends strongly on the bias polarity (figure 4). Experimentally, this was first detected for Mn adatoms adsorbed on Cu₂N [13] and successfully discussed theoretically in a second order perturbative scattering model [79, 80]. Figure 13 shows the simulated spectra of a Mn spin system at high field which—due to the small magnetic anisotropy—leads to an almost equidistant energy difference between the five eigenstates. These eigenstates are well described by pure states with the magnetic quantum numbers ${m}_{z}=-5/2,-3/2,\ldots ,+5/2$.

Zoom In Zoom Out Reset image size

At low tip sample couplings the $\rm{d}I/\rm{d}V$ spectrum is similar to the ones simulated without rate-equations (figures 9(d) and (e)), but at higher coupling we observe a drastic reduction of the differential conductance at negative bias. This decrease is concomitant with the reduction of the average magnetization $\langle {m}_{z}\rangle =\mathrm{tr}\left({\hat{S}}_{z}{{\displaystyle \sum }}_{i}{\rm{p}}_{i}| {\psi }_{i}\rangle \langle {\psi }_{i}| \right)$. For the highest T₀ the average magnetization becomes even negative at negative bias showing clearly the inversion of the state populations. Interestingly, the strong bias asymmetry in the $\rm{d}I/\rm{d}V$ spectrum is only due to the potential scattering. Simulating the system with U = 0 results in a much less asymmetric spectrum, even though the bias dependence of the state populations and average magnetization are unaltered.

Finally, we have to remark that experimentally an effective interaction of the substrate conduction electrons with the Mn spin of ${G}_{S}=2.7\ \mu \rm{S}$ was found [13]. This is much higher than the spin-substrate scattering determined solely by the spectroscopically estimated $J{\rho }_{0}\approx -0.02$ which results to ${G}_{S}=2{e}^{2}/h\times$ ${(J{\rho }_{\rm{s}})}^{2}S$ $(S+1)\approx 0.3\ \mu$ S. This means that in this particular system roughly 90% of the spin relaxations with the substrate electrons must originate from scattering processes which do not directly leave their fingerprint in the observable differential conductance spectrum. In this context, ab inito density functional calculations have found an about 3.1 times higher coupling to the substrate via the 3d_xz and $3{d}_{{x}^{2}-{y}^{2}}$ orbitals than experimentally observed [116], while the approach discussed here neglect any orbital symmetry.

5.1. Current induced correlations

Comparing the dynamical model outlined above with experimental data that require the evaluation up to third order scattering is a very intriguing test of the capability, as well as the limitations, of this approach. Thus, we will now analyze spectroscopic data that have been measured for Fe adatoms on Cu₂N using spin-averaged and spin-polarized tips [99]. Figure 14 shows the experimental data obtained for two different magnetic field directions; along the main anisotropy axis (easy axis, B_z) which leads to a strong polarization of the eigenstates along the magnetic field (figure 14(a)) and, perpendicular to it, along the hard axis (B_x) which produces only a small polarization of the lowest eigenstates (figure 14(d)).

Zoom In Zoom Out Reset image size

In this experiment the spin polarization was deliberately changed by vertical manipulation ('picking-up') [13, 117] of a Mn atom onto the apex of the tip, changing its polarization from ${\eta }^{\rm{t}}\approx 0$ to ${\eta }^{\rm{t}}\approx 0.4$ and enabling the study of the identical atom with and without a polarized tip. The experimental data for the spin-average and for the spin-polarized measurements along the hard axis can be well simulated within one set of parameters (figure 14(c)). Surprisingly, the spectra are almost identical to simulations using only the zero-current approximation. Thus, the moderate coupling between tip and sample does not disturb significantly the state populations. This is also evident by plotting the average magnetic moment along the applied field (lower panels in figures 14(b) and (c)). The small magnetic field (compared to the anisotropy energy) of only B_x = 3 T along the magnetically hard-axis cannot polarize the Fe atom significantly (figure 14(d)). This means, that the spin imbalance produced by a spin-polarized current has only little influence on the $\rm{d}I/\rm{d}V$ spectrum. However, note that the small shift toward higher energies of the step at $| V| \approx 4$ mV, which has its origin in transitions between the states $| {\psi }_{1}\rangle$ and $| {\psi }_{3}\rangle$, might be due to an interaction between the spin system on the surface and the spin polarized electrode of the tip [118, 119].

The situation changes quite drastically when the magnetic field is applied along the easy axis and the spectrum is taken with the spin-polarized tip (figure 14(b)). A strong peak at $V\approx -5$ mV appears concomitant with the apparent disappearance of the step at precisely the same energy, that was clearly visible in the spectrum measured with the spin-averaging tip. Fortunately, the model describes this spectrum quite well, too and thus allows us to discuss the physical origin of these differences.

Comparing the spectrum with the one calculated in the zero-current-approximation reveals that the disappearance of the step is indeed due to the increased coupling between tip and sample. The origin lies in the finite U, quite similar to the differential conductance reduction observed at negative bias for the Mn system (figure 13). Indeed, we see that the spectrum calculated with the same set of parameters except U = 0 leads to a spectrum in which the step-like increase of the $\rm{d}I/\rm{d}V$ at $V\lt -4$ mV is maintained. Furthermore, we notice that the calculated magnetization $\langle {m}_{z}\rangle$ decreased abruptly below this tunneling bias voltage, much faster than for the spin-averaging tip, due to the increased population of the states $| {\psi }_{2}\rangle$ and $| {\psi }_{3}\rangle$ (figure 14(b) lower panel).

While the effects discussed above give a plausible explanation for the overall shape of the spectrum, there is still a noticeable discrepancy between experiment and simulation. In particular, the experimental $\rm{d}I/\rm{d}V$ data show a significantly stronger peak at $V\approx -5$ mV than one would expect based on the simulation, resembling of the effects of strong correlations discussed in section 4 (which are not covered by the model). Thus, it is worth to dissect the characteristic peak which has its origin in the third order process (131) in which the spin system is scattered from $| {\psi }_{1}\rangle \to | {\psi }_{3}\rangle \to | {\psi }_{1}\rangle$. With the spin-polarized tip this transition has a significantly higher probability to occur when the electrons tunnel from sample to tip. The process is illustrated in figures 14(d) and (f). It starts with an $| \downarrow \rangle$ electron in the sample which scatters on the spin system changing its grounstate from $| {\psi }_{1}\rangle$ to $| {\psi }_{3}\rangle$ and then tunnels, with increased probability, as $| \uparrow \rangle$ electron into the majority states of the tip. The spin-down hole in the sample will lead to a slightly higher spin-up electron density close to the Fe atom. Next, a $| \uparrow \rangle$ electron from the sample fills the $| \downarrow \rangle$ hole together with changing the localized spin-state from $| {\psi }_{3}\rangle$ back to $| {\psi }_{1}\rangle$ (figure 14(d)). During this process the substrate electrons are correlated with the spin state. While the magnetic moment for $| {\psi }_{1}\rangle$ is approximately ${m}_{z}\approx +2$ and for $| {\psi }_{3}\rangle$ ${m}_{z}\approx +1$, respectively, at low enough temperature and at sufficiently negative bias and tunneling current, a correlated out-of-equilibrium Kondo-state develops which has the total wavefunction

$$\begin{eqnarray}&&{\Psi }^{\mathrm{total}}={\displaystyle \frac{1}{\sqrt{2}}}\left(\left|{\downarrow }^{\rm{s}}\right.\rangle| +2\rangle -\left|{\uparrow }^{\rm{s}}\right.\rangle| +1\rangle \right).\end{eqnarray} \tag{ 43 }$$

Even though the experimental data only give hints that this correlated state has formed, a similar non-equilibrium Kondo formation has been found in transport measurements on carbon nanotubes [120]. The physical properties of such a non-equilibrium correlated state are very intriguing [121, 122]. Entering this state, the substrate electrons partly screen the magnetic moment of the Fe spin reducing it to $S-1/2$. Such an underscreened Kondo state should show a very particular temperature and energy dependence [123–129]. Interestingly, in a system like the one discussed here, the transition between weak and strong coupling is not only governed by a characteristic temperature T_K [130], but the formation of this correlated state is deliberately tuned by drastically increasing the probability of its creation when applying a spin-polarized current of sufficient energy by the probing tip. This enables us to envision for example pump-probe experiments [14] in which the time-evolution of the formation and the decay of this correlated state is measured in detail.

In this manuscript I have shown that applying perturbation theory to quantum spin systems enables one to describe experimentally measured differential conductance spectra with very high accuracy. This enables one to obtain a profound understanding of the physical processes on play and to separate single- as well as many-electron effects.

The versatility of low-temperature scanning tunneling measurements on single and complexly coupled spin systems led me to believe that we should expect a multitude of exciting new experiments for the future, which will further deepening our fundamental knowledge on quantum systems in general and, in particular, quantum magnetism. Perturbative models, like the one outlined here, might be of help in such systems. For that, the supplemental material to this manuscript include an easy usable software package that allows not only to simulate the differential conductance spectra of arbitrary complex spin systems, but additionally allows one to fit experimental data to the model.

In this manuscript we have restricted ourselves to the adiabatic limit and neglected any time dependence in the parameters. However, it is straightforward to expand this model to capture also time dependent pump-probe measurements. In such a framework, coherent state superpositions could be accounted for by using, for example, a Bloch–Redfield approach in which the spin system is coupled to an open quantum-system [114, 131, 132] and in which interactions up to third order are included to account not only for the decay but also for the creation of coherences.

Additionally, the study of the coupling between the spin-system and other (quasi)-particles, like phonons, photons, or magnons would be very interesting. Furthermore, while the perturbative approach fails when the system enters the strong-coupling Kondo-regime, it would be very intriguing to combine the simple model outlined here, with exact quantum impurity models.

I thank Peter Wahl for invaluable discussions about the original Appelbaum paper. Furthermore, I am very thankful for fruitful discussions and for providing me with experimental data to Christian Ast, Katharina Franke, Andreas Heinrich, Cyrus Hirjibehedin, Alexander Khajetoorians, Nicolas Lorente, Sebastian Loth, Alexander Otte, and Marteen Wegewijs. This work was supported by the SFB 767.

The two arbitrary spin density matrices ${\varrho }^{\rm{t}}$ and ${\varrho }^{\rm{s}}$ are a full description of the ensemble states in the tip and sample electron bath close to the Fermi energy. The eigenvectors $\left|{\varphi }_{i}^{\rm{t}}\right.\rangle$ and $\left|{\varphi }_{i}^{\rm{s}}\right.\rangle$ of ${\varrho }^{\rm{t}}$ and ${\varrho }^{\rm{s}}$ are representative eigenstates of these incoherent ensembles and enables one to calculate the interaction transition intensities between them as:

$$\begin{eqnarray}&&{\xi }_{i\prime f\prime }^{x,y,z,I}=\sqrt{{\lambda }_{{i}^{\prime }}^{\rm{t}}{\lambda }_{{f}^{\prime }}^{\rm{s}}}\left.\langle {\varphi }_{{f}^{\prime }}^{\rm{s}}\right|{\hat{\sigma }}_{x,y,z,I}\left|{\varphi }_{{i}^{\prime }}^{\rm{t}}\right.\rangle,\end{eqnarray} \tag{ A.1 }$$

with ${\lambda }_{i}^{\rm{t}}$ and ${\lambda }_{i}^{\rm{s}}$ as the eigenvalues of the corresponding eigenvectors. Note, that due to the in general incoherent spin ensembles in tip and sample equation (A.1) has to be evaluated independently for all combinations of $i\prime =1,2$ and $f\prime =1,2$.

For the localized spin system we can write similarly the transition intensities between the eigenstates:

$$\begin{eqnarray}&&{\Xi }_{{if}}^{x,y,z}=\left.\langle {\psi }_{f}\right|{\hat{S}}_{x,y,z}\left|{\psi }_{i}\right.\rangle,\end{eqnarray} \tag{ A.2 }$$

so that the absolute square of the matrix element in equation (13) calculates to:

$$\begin{eqnarray}&&{\left|{\mathcal{M}}_{(1){if}}^{\rm{t}\to \rm{s}}\right|}^{2}={\displaystyle \sum }_{i\prime ,f\prime }{\left|{\displaystyle \frac{1}{2}}\left({\xi }_{i\prime f\prime }^{x}{\Xi }_{{if}}^{x}+{\xi }_{i\prime f\prime }^{y}{\Xi }_{{if}}^{y}+{\xi }_{i\prime f\prime }^{z}{\Xi }_{{if}}^{z}\right)+U{\xi }_{i\prime f\prime }^{I}{\delta }_{{if}}\right|}^{2}.\end{eqnarray} \tag{ A.3 }$$

The interference term between ${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{2}$ (equation (23)) can be evaluated in an analogous way leading to:

$$\begin{eqnarray}&&{M}_{{fi}}{\tilde{M}}_{{mf}}{M}_{{im}}={\displaystyle \sum }_{\underset{\{x,y,z,I\}}{j,k,l=}}\;{\displaystyle \sum }_{i\prime ,m\prime ,f\prime }{\bar{\xi }}_{f\prime i\prime }^{l}{\Xi }_{{fi}}^{l}\times {\tilde{\xi }}_{m\prime f\prime }^{k}{\Xi }_{{mf}}^{k}\times {\xi }_{i\prime m\prime }^{j}{\Xi }_{{im}}^{j},\end{eqnarray} \tag{ A.4 }$$

and for the reversed order processes:

$$\begin{eqnarray}&&{M}_{{fi}}{M}_{{mf}}{\tilde{M}}_{{im}}={\displaystyle \sum }_{\underset{\{x,y,z,I\}}{j,k,l=}}\;{\displaystyle \sum }_{i\prime ,m\prime ,f\prime }{\bar{\xi }}_{f\prime i\prime }^{l}{\Xi }_{{fi}}^{l}\times {\xi }_{m\prime f\prime }^{k}{\Xi }_{{mf}}^{k}\times {\tilde{\xi }}_{i\prime m\prime }^{j}{\Xi }_{{im}}^{j}.\end{eqnarray} \tag{ A.5 }$$

The tilde ($\tilde{\xi }$) and bar ($\bar{\xi }={\xi }^{\dagger }$) account for processes which starts and ends in the sample, and for processes in which tunneling electrons traverse the junction in the opposite direction, i.e. which starts in the sample and ends in the tip, respectively. They can be calculated in an analogous way with equation (A.1).