Spin Excitations and Correlations in Scanning Tunneling Spectroscopy

In recent years inelastic spin-flip spectroscopy using a lowtemperature scanning tunneling microscope has been a very successful tool for studying not only individual spins but also complex coupled systems. When these systems interact with the electrons of the supporting substrate correlated manyparticle states can emerge, making them ideal prototypical quantum systems. The spin systems, which can be constructed by arranging individual atoms on appropriate surfaces or embedded in synthesized molecular structures, can reveal very rich spectral features. Up to now the spectral complexity has only been partly described. This manuscript shows that perturbation theory enables one to describe the tunneling transport, reproducing the differential conductance with surprisingly high accuracy. Well established scattering models, which include Kondo-like spin-spin and potential interactions, are expanded to enable calculation of arbitrary complex spin systems in reasonable time scale and the extraction of important physical properties. The emergence of correlations between spins and, in particular, between the localized spins and the supporting bath electrons are discussed and related to experimentally tunable parameters. These results might stimulate new experiments by providing experimentalists with an easily applicable modeling tool.


Introduction
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 quantummechanical 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 2 O 3 on a NiAl surface [6]. Since then many experiments have focused on transition metal atoms on a thin layer of Cu 2 N on Cu (100) [1,[7][8][9][10][11][12][13][14][15][16][17][18][19][20]. In these experiments, patches of Cu 2 N are formed by sputtering a clean Cu (100) surface with N + ions and subsequent annealing. This leads to a monolayer of Cu 2 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 I V d d . All spectra show characteristic features that can be classified into two groups: (1) step-like increases in the differential Apart from Cu 2 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 3 and a high magnetic anisotropy of 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 2 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 2 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 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 2 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][64][65][66][67][68][69][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][75][76][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][80][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][87][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]: 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  can have both signs, either positive or negative. However, this result can be different for higher spins.
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 2 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.

The model
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: The Hamiltonian of the tip Ĥ t and sample Ĥ s , respectively, can be described in the framework of second quantization using creation and annihilation operators â † and â: with k t ϵ σ and k s ϵ σ 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 a a ( )ˆk k , with · 〈 〉as the time averaged expectation value. In general, the electronic states φ | 〉 in these electrodes might be spin-polarized in an arbitrary direction which we account for by the corresponding spin density matrices I n Here, n⃗ describes the direction and n 1 1 − ⩽ | ⃗ | ⩽ the amplitude of the polarization in the chosen coordinate system, (ˆ,ˆ,ˆ) x y z σ σ σ σ = are the standard Pauli matrices, and Iˆis the (2 2) × identity matrix. With this convention the spin polarization is identical to the relative imbalance between majority and minority spin densities , where ↑ and ↓ 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 Ĥ a that includes Zeeman and magnetic anisotropy energy and-in the case of more than one spin-the Heisenberg coupling terms ij  ⃗ 1 and the non-collinear Dzyaloshinskii-Moriya coupling ij  ⃗ between individual spins: In this equation 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 2 . The externally applied magnetic field is B ⃗ and the total spin operators S S S Ŝ (ˆ,ˆ,ˆ) Diagonalizing the Hamiltonian of the localized spin system (equation (5)) leads to discrete energy eigenvalues n ϵ and eigenstates n ψ | 〉 . For a single spin there are S 2 1 + eigenstates, while for complex spin structures that consist of several spins the number of eigenstates increases quickly as S (2 1) i i ∏ + . 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 φ | 〉 and the discrete spin states ψ | 〉, i.e. , t,s φ ψ | 〉 t,s φ ψ =| 〉| 〉. 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 1 We write ij  ⃗ as vector to enable also anisotropic Heisenberg or Ising-like couplings. 2 Note, that instead of using the rather phenomenological D and E values the model can be easily adapted to more physical operators that connect to the spin-orbit couplings [15] or to the symmetries of the system [61]. 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)): with T ta i as the coupling constants between tip and the ith adsorbate spin, and U J 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: 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][96][97].

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 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 if f 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: 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 | contains all experimental parameters that will determine the strength of the tunneling current. In particular, T 0 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 1 0 2 ≪ 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 tipsample 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 (zerocurrent 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, 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- 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]: , so that equation (9) becomes: The main goal of this paper lies in deriving the transition matrix elements 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): This matrix is energy independent and connects the initial states in the electron baths i φ | 〉 and spin system i ψ | 〉 with their final states f φ | 〉and f ψ | 〉. 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 ( if δ is the delta distribution). Computing the absolute square of if (1)  results in three terms that contribute to the tunneling conductance: denotes the real part of the matrix A. The first term consist of the operators 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: where, for completeness, we have additionally included the probability amplitudes for interacting with the unity operator Iˆ. Since ŝ 2σ = , 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]: 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 η = along the z-axis, the transition matrix elements are tunneldirection dependent because either the initial or the final electron bath has now a spin imbalance in the density of states [99]: 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 g B B μ | | > | |. These steps are symmetrical around zero voltage and thermally smeared out by about k T 5 B [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 ( 0 η > ) the conductance step at positive bias decreases, while it increases by the same amount at negative bias ( figure 4(b)).
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 1 η = and U S 1 2 = − (figures 4(c)-(d)). This magneto-resistive elastic tunneling [99] arises from the third term in equation (14) and scales with the expectation value Sˆz 〈 〉 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].

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: Here, the tilde on the matrix 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 m ψ | 〉 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 ψ | 〉 and an integration over the continuous states of the intermediate electron states φ | 〉 in the substrate. For the integration, we consider states in an energy range of 0 ω ± around E F as possible scatterer leading to the following characteristic function for electron-like processes [97,100]: and for hole-like processes: 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 0 ϵ = leads to a logarithmic singularity that is broadened by the temperature. Thus, as long as the energy m ϵ of the intermediate state and the temperature k B T are small compared to the integration bandwidth 0 ω , the equations (20) and (21) (with a change of sign) can be rewritten as: Figure 5. Here, ∂ is the derivative of the temperature broadened step function (equation (11)) and a small 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 B T 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 . Due to the restriction to states in an energy interval 0 ω ± the function remains analytical even at ϵ → ±∞ where F approaches zero. Note, that the precise value of the cut-off energy 0 ω is not critical and mainly changes only the background offset. If not otherwise noted we use 20 0 ω = meV throughout the paper.
When calculating the conductance we now have to consider both processes depicted in figures 3(b) and (c) and thus have to in equation (12) leading to: The evaluation of the matrix elements up to third order yields two new terms due to the interference between the processes described by 1  and 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 ( 0) ≠ the term (23b) will, in addition, produce a bias-asymmetry in the differential conductance even without spinpolarized 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 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: , , Here, equation (24) accounts for the direct diagrams (normal order) and (25) for the exchange diagrams (reversed order), respectively. In these equations 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.
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 Sˆ0 z 〈 〉 ≠ (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 | ↓ 〉 electron originating from the tip is scattered into a | ↑ 〉 state in the sample exciting the spin system from 1 2 ψ ψ | 〉 → | 〉. Second, a | ↑ 〉 electron in the sample is scattered into | ↓ 〉 bringing the spin system back to 1 ψ | 〉. Third, the sample | ↓ 〉 is scattered back into the tip as | ↓ 〉 electron. The last step of this process can either take place via theˆz σ or, in the case of potential scattering, via theˆI σ term, respectively. For an electron tunneling in the reverse direction we first have a | ↓ 〉 from the sample being scattered into an electron | ↑ 〉 state of the tip simultaneously exciting the spin system from 1 2 ψ ψ | 〉 → | 〉. Then, the second process flips a | ↑ 〉 sample electron into the | ↓ 〉 hole from the first process. Finally, a | ↑ 〉 tip electron traverses the junction and fills the | ↑ 〉 hole in the sample. Comparing both tunneling directions we see Figure 6. Interaction diagrams of order two (a) and three (b) for an electron tunneling from tip to sample into a two level S 1 2 = spin system. The large (small) spheres depict the state of localized spin (interaction electron) and the color their spin directions. Schematic spectra show their contributions to the conductance at positive bias. The numbers label the processes with the state order of the localized spin-system. An appended 'R' label processes in which the scattering into the intermediate state is performed before the tunneling electron interacts (exchange diagrams). Note that the time order of the processes influences crucially the conductance spectra as schematically displayed in the small graphs (vertical line is E F ; the * means multiplication).
that it is the third process which differs in the initial and final state, i.e. the matrix elements are eitherˆz σ 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:  where we assumed J 0 s ρ < and U 0 > . Note, that the preceding sign change at the tunneling direction from sample to tip (s t → ) is due to the rearrangement of the interaction order together with the switching from holelike to electron-like scattering. The contribution to the conductance from the processes in which only Kondolike 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.

Some single spin examples
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][96][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][84][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.

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 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 0 < ) the energetically lowest eigenstates have weights only for m 1 z = ± . 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:  (110) surface [65], for Fe atoms adsorbed on InSb(110) [53], or for CoH complexes adsorbed on the h-BN/Rh (111) surface [51]. 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 Sˆ+, Sˆ−, and Sˆz and are thus forbidden. Nevertheless, additional potential scattering enables transitions involving the operators Sˆ+, Sˆ−, and Iˆor Sˆz, Sˆz, and 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 2 ϵ = ± , but due to energy conservation it contributes to the spectrum only at eV 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 2 ϵ and 3 ϵ ( figure 8(c)). Note, the overall general behavior would not alter when changing from easy-axis to easy-plane anisotropy (D 0 > ) as long as all degeneracies are lifted.

Single Mn and Fe atoms on Cu 2 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 2 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). 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 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 2 |± 〉as ground and first excited state and weights in 1 |± 〉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 m 2 Δ = ± . Experimental I V d d 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 2 ϵ and 3 ϵ or 4 ϵ . In contrast, the processes (132) and (142) scale with J 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 3 . 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 0.087 s ρ = − , 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 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 2 N, which have a spin of S 5 2 = and only a small easy-axis anisotropy of D 40 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 5 2 z = ± (figure 9(c)). At zero field a typical spectrum shows only one step, which belongs to the transition between the 5 2 ± and the 3 2 ± states that have superimposed asymmetric peak structures ( figure 9(d)). The fit to the model yields J 0.029 s ρ = − 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 0.0091 s ρ = − . Interestingly, we find for both atoms a potential scattering value of U S 1 2 ≈ , 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 1 | | < mV is evident while the differential conductance at V 1.5 > mV stays, for all tips, constant.
4. The limit of the perturbative approach: the Kondo system Co on Cu 2 N When a half-integer spin with S 1 2 > has easy-plane anisotropy D 0 > , its ground state at zero field is a doublet with its main weights in m 1 2 z = ± . 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 2 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 2 N at temperatures down to T = 550 mK [9]. Apart from D 0 > the system also has a small in-plane anisotropy (E 0 ≠ ) which creates an easy axis (x) along the nitrogen row and a hard axis (z) along the vacancy rows ( figure 10(a)).
Interestingly, in this system the coupling to the substrate J s ρ changes with the position of the Co atom on larger Cu 2 N patches, concomitant with a change in the anisotropy energies which separates the 1 2 |± 〉 states from the 3 2 |± 〉 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 0.1 s ρ ≈ − ), 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 D 2 | | = | |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 2 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 0.25 s ρ ≈ − , 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 D 2| |. 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 2 N system have weights in states that are separated by m 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 i ψ | 〉and f ψ | 〉of the spin system that have the energies i ϵ and f ϵ is Γ ≔ , i.e. the transition probability per time, is analogous to the calculation of the current (equation (9)) given by: Here we have set the bias to eV = 0 and replaced the coupling constant T 0 between adsorbate and tip with J s ρ, the coupling between adsorbate and substrate. When we now consider only processes up to second order, the matrix if  is independent of the energy and the integral can be evaluated to Thus, the scattering rate between the two degenerate states (i.e. 0 12 ϵ = ) is directly proportional to the temperature ( figure 11(a)): B 12 (2)s s s 2 12 2 Γ π ρ = →  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, k T B decoh Γ ∝ [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: if m fi mf im Due to the growing intensity of F ( 0) ϵ = at reduced temperatures (figure 5), for temperatures T 0 → the scattering (3) Γ decreases significantly more slowly than (2) Γ (see figure 11(b)) so that their ratio steadily increases: In contrast to that, the decoherent scattering rate with the bath 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 (3) Γ and higher order scattering terms become the dominant processes and perturbation theory breaks down [34,36,106,107]: 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.
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 2 N system. At temperatures T 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 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 k T B 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 3 2 z = ± 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.

Spin dynamics and rate equations
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 I V d d 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 i ψ | 〉and f ψ | 〉due to the tunneling current flowing between tip and sample we can again use equation (9) and the relation eV , which leads, in first order Born approximation, to: The integral can be solved using equation (29) resulting in: (36) can be further simplified to an equation that is linear in the energy difference: 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 0 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 ( ) t, 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 1  and 2  of equation (23) and leads to an additional contribution to the scattering rate of: 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)  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 i τ for the eigenstate i ψ | 〉 is inversely proportional to the sum of all scattering processes which leaves the state: 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 2 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 5 2, 3 2, , 5 2 z = − − … + . At low tip sample couplings the I V d d 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 For the highest T 0 the average magnetization becomes even negative at negative bias showing clearly the inversion of the state populations. Interestingly, the strong bias asymmetry in the I V d d 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 2.7 S S μ = was found [13]. This is much higher than the spin-substrate scattering determined solely by the spectroscopically estimated J 0.02 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 d 3 x y 2 2 − orbitals than experimentally observed [116], while the approach discussed here neglect any orbital symmetry.

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 2 N using spin-averaged and spinpolarized 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)).
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 0 t η ≈ to 0.4 t η ≈ 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 I V d d spectrum. However, note that the small shift toward higher energies of the step at V 4 | | ≈ mV, which has its origin in transitions between the states 1 ψ | 〉 and 3 ψ | 〉, 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 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 I V d d at V 4 < − mV is maintained. Furthermore, we notice that the calculated magnetization m z 〈 〉 decreased abruptly below this tunneling bias voltage, much faster than for the spin-averaging tip, due to the increased population of the states 2 ψ | 〉and 3 ψ | 〉 ( 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 I V d d data show a significantly stronger peak at V 5 ≈ − mV than one would expect based on the simulation, resembling of Even though the experimental data only give hints that this correlated state has formed, a similar nonequilibrium 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][124][125][126][127][128][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.

Summary
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 manyelectron 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 pumpprobe 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.