Predicting the conductance of strongly correlated molecules: the Kondo effect in perchlorotriphenylmethyl/Au junctions

Abstract

Stable organic radicals integrated into molecular junctions represent a practical realization of the single-orbital Anderson impurity model. Motivated by recent experiments for perchlorotriphenylmethyl (PTM) molecules contacted to gold electrodes, we develop a method that combines density functional theory (DFT), quantum transport theory, numerical renormalization group (NRG) calculations and renormalized super-perturbation theory (rSPT) to compute both equilibrium and non-equilibrium properties of strongly correlated nanoscale systems at low temperatures effectively from first principles. We determine the possible atomic structures of the interfaces between the molecule and the electrodes, which allow us to estimate the Kondo temperature and the characteristic transport properties, which compare well with experiments. By using the non-equilibrium rSPT results we assess the range of validity of equilibrium DFT + NRG-based transmission calculations for the evaluation of the finite voltage conductance. The results demonstrate that our method can provide qualitative insights into the properties of molecular junctions when the molecule–metal contacts are amorphous or generally ill-defined, and that it can further give a fully quantitative description when the experimental contact structures are well characterized.

---

1. Introduction

Molecular electronics holds great promise for future applications in computing, sensing, clean-energy, and even data-storage technologies.^1–3 However, a general difficulty so far has been the poor characterization of the device structures and their relationship with the measured conductances and functionalities. For this problem, ab initio simulations based on density functional theory (DFT)⁴ have proven very successful in supporting experiments, and they have played a key role in advancing the field during the last decade.^5–9 Yet, standard DFT-based transport schemes for simulations of experimental molecular junctions have several limitations. The most prominent of these is the failure to account for the strong electron correlations leading to the Kondo effect in devices comprising magnetic molecules, and rigorous treatments and extensions overcoming this problem are currently under active development.^10–15

In this article, we establish a suitable combination of DFT and many-body techniques to achieve an unprecedented quantitative description of the equilibrium and non-equilibrium conductance of molecular devices showing Kondo effect. By using gold/perchlorotriphenylmethyl (PTM)/gold junctions as a specific example we relate the Kondo temperature to the electrode–molecule contact geometries, thus matching the range of variability of the experimental results.¹⁶ Furthermore we address the dependence of the conductance at finite temperature and extend the method to finite bias.

Our multi-scale approach combines DFT, non-equilibrium Green's functions (NEGF),¹⁷ numerical renormalization group (NRG) methods,^18–21 and renormalized superperturbation theory (rSPT).^22,23 First the contact geometry and electronic structure of molecular junctions are obtained by DFT + NEGF. Then the DFT Kohn–Sham (KS) states are projected onto an effective Anderson impurity model,^24–33 which is solved exactly to obtain the Kondo temperature and the equilibrium zero-temperature conductance via NRG. Based on these results we finally compute the non-equilibrium rSPT transport coefficients, which encode the behavior of the junctions at low temperature, finite magnetic field, and finite bias voltage.²²

Stable organic radicals contacted to metal electrodes, such as the PTM molecule on Au, form a practical realization of the prototypical single-orbital Anderson impurity model,^16,31,34,35 and are therefore ideally suited to study the fundamental aspects of the interaction of magnetic impurities with metallic surfaces. These aspects include the interplay between the binding geometry and the energy level alignment with respect to the surface Fermi energy, as well as the electron correlations leading to the Kondo effect.

In recent experiments^16,36 PTM-radicals were functionalized with thiophene linkers producing the PTM-bis-thiophene radical (called PTM-BT in the following to distinguish it from the bare PTM; see also Fig. 1 for their atomic structures). These molecules were then integrated into gold mechanically-controlled break-junctions (MCBJs) and gold electromigrated break-junctions (EMBJs) to measure their transport properties. While at room temperature very low conductance values were reported,³⁶ at low temperature a zero-bias conductance resonance was observed in many of the junctions, and its Kondo character verified by temperature- and magnetic field-dependent measurements.¹⁶ The low-temperature results indicate that the PTM radical can preserve the unpaired spin in a solid state three-terminal configuration, and that it is stable under mechanical stretching of the electrodes. One of the remarkable features is the rather high Kondo temperature of about 3 K, which is largely constant upon stretching of the junction. This implies that for the junctions that exhibiting Kondo behavior the contact of the molecule to one of the electrodes is very strong, and is not affected by the elongation of the junction in the MCBJ process. In contrast, the background conductance shows large variations. This can happen upon stretching when the contact to the second electrode varies significantly, or else when one of the two electrodes changes its Au–Au bond conformation significantly.^37,38 Overall the low-temperature experimental results point to a structure with highly asymmetric coupling to the electrodes. In the following we will show that this hypothesis is indeed confirmed by our calculations, thus providing a detailed understanding of the electronic and transport properties of the PTM/gold junctions at the atomic scale.

Fig. 1 Relaxed atomic structures of the bare PTM molecule (a) and of PTM-bis-thiophene (b) (green spheres represent Cl atoms, blue spheres H, large yellow spheres S, and smaller dark yellow spheres represent C).

The paper is organized as follows. We first discuss the equilibrium DFT results for a number of possible junction structures (section 2), and then provide estimates for the Kondo temperature for these geometries (section 3). For a set of geometries we then present the linear response transport properties including the strong electron–electron correlations obtained by DFT + NEGF + NRG (section 4) and finally extend the results to finite temperature and finite bias via rSPT (section 5).

2. DFT calculations

PTM has a propeller-like structure with a central carbon atom coordinated by the three phenyl rings. In the gas phase, it has the typical electronic structure of a radical.^39–41 The energy spectrum has doubly occupied electronic states filled up to the highest occupied molecular orbital (HOMO). Above the HOMO there is a further well-separated, singly occupied molecular orbital (SOMO) with an unpaired electron, giving a total molecular spin quantum number of 1/2. In the PTM, the charge isosurface indicates that the SOMO is mainly confined to the central carbon, while the HOMO and the lowest unoccupied molecular orbital (LUMO) are largely located on the rest of the molecule. This is presented more extensively in the ESI section S2,† while the computational details of our DFT calculations are given in the ESI section S1.† The difference between the ionization potential and electron affinity of the molecule defines the fundamental gap and corresponds to the charging energy U. In the absence of any experimental results, we calculate U via total energy differences⁴² to be about 4 eV. PTM-BT has a very similar electronic structure to that of the bare PTM, although the SOMO is slightly delocalized over the thiophene ligands,¹⁶ and this results in a charging energy smaller by about 0.4 eV. Note that when the molecule is placed between Au electrodes there is a significant renormalization of the energy levels and consequently a reduction of the charging energy, which we discuss in section S2 of the ESI as well as in section 3.2.†

In order to understand the electronic structure of the molecule/Au contact and how this determines the key parameters affecting the Kondo temperature, we consider a number of qualitatively different model structures, which are shown in Fig. 2.

Fig. 2 Junction geometries for bare-PTM on an Au surface (B1–B4), and for PTM-BT between two Au electrodes (T1–T8), investigated in this paper. For each junction we specify the broadening of the singly occupied molecular orbital induced by the coupling to the electrodes (Γ), its position with respect to E_F (ε), and the coupling to the left and right electrodes (Γ_L and Γ_R, respectively; Γ = Γ_L + Γ_R). All units of the specified quantities are meV.

To start, we look at the ideal case of a bare PTM molecule on a flat Au(111) surface, which we denote as configuration (CFG) B1 in Fig. 2. The 3-atom Au tip is placed at a rather large distance, so that the electronic coupling between the molecule and the tip is negligible with respect to that to the substrate. Since in MCBJ and EMBJ experiments the Au stretched surface is expected to be highly corrugated rather than perfectly flat,^37,38 we then model a rough Au surface by removing a number of Au atoms from the perfect Au(111) surface (CFGs B2 to B4). Finally, we consider a number of break-junction setups comprising PTM-BT (CFGs T1 to T8). The detailed contact structure is expected to be different for each individual experimental conductance trace measurement. The model junctions considered here include cases with both symmetric and asymmetric molecule–electrodes coupling. For some structures the PTM central core is located inside the junction's empty gap, whereas for other structures it is physisorbed on one of the electrodes. Furthermore, the thiophene linkers can be connected to the electrodes either non-covalently or covalently via a sulfur–Au adatom direct bond.

A representative DFT projected density of states (PDOS) is shown in Fig. 3 (see ESI section S1† for the computational details). When the molecule is in contact with the Au electrodes, the SOMO DOS can be modeled approximately by a half-filled Lorentzian-like peak close to the Fermi energy, E_F. Note that while we refer to the state as SOMO also when the molecule is on the Au substrate for consistency, its occupation can generally deviate from one in this case. The full width at half maximum (FWHM) of the SOMO peak corresponds to its electronic coupling to the Au substrate, Γ,²⁴ which can be calculated by using the projection scheme recently developed in ref. 24. The results for each model geometry considered are indicated in Fig. 2 along with the DFT SOMO on-site energy, ε, relative to E_F. These values are the parameters required for the evaluation of the Kondo temperature.

Fig. 3 LDA projected density of states (PDOS) for the configuration T1. The peak at the Fermi energy corresponds to the singly occupied molecular orbital, which defines our Anderson impurity, and is located mainly on the central carbon atom of the PTM.

As a matter of notation we label the structures with the propeller-like PTM parallel (perpendicular) to the surface, as “parallel” (“perpendicular”) configurations. For the idealized case of a bare PTM on a flat Au(111) surface, we find that the molecule is physisorbed with an energy difference between the “parallel” configuration (B1) and the “perpendicular” configuration (not shown) of about 335 meV, favoring the “parallel” configuration. The equilibrium position of the central C atom is located at about 5.18 Å from the top Au layer. For this configuration there is a negligible charge-transfer from the surface to the molecule, and the PTM preserves its unpaired electron, with Γ ≈ 7 meV and therefore very small.

On the corrugated surface (CFG B2) the molecule can bind better to the Au, since part of its phenyl rings can move into regions where the Au surface has a dip. In CFG B2 an Au atom is located below the central C atom of the bare-PTM. This atom is then removed in the CFG B3, while it is kept as the only atom from the top-most Au surface in the CFG B4. Comparing the Γ-values for these structures allows us to estimate the effect of Au atoms directly in contact with the central C atom of the PTM. For CFG B2 we find the occupation of the SOMO to be 1.40 electrons, indicating that a partial electron transfer between the gold and the molecule has occurred. In fact, the SOMO DOS peak lies below E_F (ε = −54 meV). The increased charge transfer indicates an increased screening of the transferred electrons by the Au surface atoms, which is due to the molecule moving closer to the Au surface, in particular to the Au atom closest to the core of the PTM molecule. In general an increase in the screening also leads to a reduction of U (see the discussion in the ESI section S2†). We note that the results for the charge transfer obtained for non-spin-polarized calculations are approximately the same as those obtained in spin-polarized calculations in the ESI section S2.† The electronic coupling of 114 meV for CFG B2 is much larger than the one for the PTM on flat Au(111). An analysis of the origin of such a large coupling shows that it is mainly due to the Au atom underneath the central C atom of the PTM. In fact, for CFG B3, where this central Au adatom is removed, the coupling drops to 23 meV, while it remains large for CFG B4, where only this Au adatom is kept of the top Au surface layer.

Finally, we consider the model break-junctions (CFGs T1 to T8). In these cases, we use only the “perpendicular” configuration, since the “parallel” PTM-BT configuration would require very large simulation cells, which are beyond our current computing resources. The results allow us to infer the general trends for the electronic coupling of the radical center to the Au electrodes through the thiophene linkers (see Fig. 2). As can be seen the computed values of Γ vary over almost one order of magnitude, from 26 to 126 meV. We note that in break-junctions Γ is the sum of two contributions, Γ_L and Γ_R, representing the electronic coupling to the left and right lead, which we calculate individually with the method outlined in ref. 24. In general we find that Γ_L or Γ_R are large when there are Au atoms close to the thiophene linkers, such as for Γ_L in CFG T7. A bond between the sulfur atoms and a protruding Au atom also increases the coupling. On the other hand, the coupling is low when such a bond is absent, and when the angle between the thiophene and the Au is larger, such as for Γ_L in CFG T1.

3. Kondo effect

3.1. Formulation of the single impurity Anderson problem

The PTM in contact with the leads is modeled by a single impurity Anderson Model (SIAM), which has the Hamiltonian⁴³

H_SIAM = H_d + H_c + H_hyb, (1)

where H_d describes the electrons of spin σ localized at the impurity site, which are created (annihilated) by the operator d^†_σ (d_σ), with n_dσ = d^†_σd_σ being the corresponding number operator; ε_d is the orbital energy, U the charging energy, and the occupation, where the bracket 〈⋯〉 denotes the thermal expectation value. For the PTM molecule the impurity site is the SOMO. Note that ε_d does not coincide with ε in Fig. 3, since the on-site Coulomb interaction is already partially accounted for in KS-DFT, and this contribution has to be subtracted, so that ε_d = ε − ε_dc.^24,44 Here ε_dc is the so-called double counting correction, whose exact expression is not known except for certain limiting cases, and several approximations have been introduced in the literature.⁴⁵ In general ε_dc depends on U, and in the commonly used “fully localized limit” it is given by ε_dc = U(n^DFT_d−1/2),¹² where n^DFT_d is the DFT occupation of the impurity. A comprehensive discussion of the difficulties arising when combining DFT with such an Anderson impurity model and more generally the dynamical mean field theory is given in ref. 45 and 46. Note that instead of the Anderson impurity model one can also use other methods to treat the highly correlated subsystem, such as for example embedded correlated wavefunction schemes.^47,48 A review of the advantages and limitations of various embedding schemes that link many-body calculations for a subsystem to an environment treated at the DFT level is given in ref. 49.

Since in an experimental setting the occupation can be set by applying a gate voltage, here we treat ε_d as an adjustable parameter, independent of the DFT results, and choose its value to ensure a specified occupation of the impurity orbital. We will also investigate how the results depend on the charging energy U, and will provide estimates of possible values of U for PTM/Au geometries.

In eqn (1)H_c describes the effective bath of electrons with momentum k and spin σ, which are created (annihilated) by the operator c^†_{k, σ} (c_k,σ) and with number operator n_k,σ = c^†_k,σc_k,σ. The effective bath includes the electrons in the Au leads, as well as those in the molecular orbitals, except for the SOMO. We note that the chemical potential in the Hamiltonian eqn (1) is set to zero by shifting both the bath and impurity energies ε_d and ε_k by an additive constant. This does not affect the properties of the system. Furthermore, in zero-temperature calculations we will refer to the chemical potential μ = 0 as the Fermi energy E_F = 0, which is most commonly used in first principles calculations.

Finally, H_hyb accounts for the hybridization between the bath and the impurity, with V_k corresponding to the hybridization matrix element. Accordingly, we can define the hybridization function Δ(E) = ReΔ(E) + iImΔ(E), with

and the coupling strength Γ(E) = −2ImΔ(E). The DFT results for Γ(E_F = 0) for several PTM/Au contacts are presented in the previous section, and the results are shown in Fig. 2.

3.2. Estimation of the Kondo temperature

In order to obtain a first estimate of the Kondo temperature for different junctions presented in Fig. 2, we assume a constant (energy independent) coupling Γ = Γ(E_F = 0). The Anderson model then maps onto the Kondo model while approaching the so-called local moment limit, where |ε_d| ≫ Γ and |ε_d + U| ≫ Γ, (n_d ≈ 1)⁵⁰ (see ESI section S3† for details), and the Kondo temperature is given by the Haldane equation^51,52with −U ≤ ε_d ≤ 0. The results obtained with this expression are shown in Fig. 4, and are compared with the NRG calculations in the next subsection. The experiments in ref. 16 show that the SOMO of the PTM is close to half-filling and that it can be brought to exact half-filling by applying a gate voltage to the system. Here we therefore consider only this half-filled case, and the effects of small deviations from half filling are presented in section 5. We note that if the molecule is partially charged, then in general θ_L increases compared to the charge neutral state,²⁴ so that the values for half filling represent a lower limit for the theoretical results. The values of Γ for each structure are taken from Fig. 2, and at half-filling for a particle-hole symmetric SIAM we have ε_d = −U/2. For large enough U all curves in Fig. 4 decay exponentially with U, and the slope of the exponential decay is inversely proportional to Γ, so that the configurations with the largest Γ have the slowest decay, and therefore the highest θ_L, for a given value of U.

Fig. 4 Kondo temperature, θ_L, as function of the charging energy, U, obtained using eqn (4) with ε_d = −U/2, and NRG solutions for the exact system- and energy-dependent hybridization function. The corresponding electrodes–molecule configurations are illustrated in Fig. 2. The horizontal orange region illustrates the experimental range of Kondo temperatures of about 1–3 K.¹⁶

Experimentally it is found that θ_L is approximately constant upon stretching of the junction, which indicates a highly asymmetric coupling, where the molecule preserves the contact geometry to one electrode, while the contact with the other electrode is elongated. In other words, the molecule is strongly bound to one of the electrodes, which may correspond to the core of the PTM-BT lying flat on a rough Au surface with the thiophene linkers bridging both sides of the junction.

This conclusion is supported by calculations for the CFGs B2, B4 and T4 structures, which have the largest values of Γ, and which all show asymmetric couplings. The calculated θ_L values lie in the experimental range if U is equal to about 1 eV. This charging energy is considerably smaller than the gas phase value of about 4 eV, and we ascribe this reduction of U to the charge screening by the electrons in the Au surface (see ESI section S2†). A value for the change of U due to screening can be calculated using a number of methods,⁵³ for molecules on general corrugated and irregular metal surfaces constrained DFT (cDFT) has been shown to give good results.^54,55 Alternatively, here we estimate it by approximating the metal surface as a plane, and by using a classical image charge model with a molecule between two metal electrodes⁵⁴ to capture this effect. In this way we calculate that a gap reduction of about 3 eV corresponds to ideal planar Au electrodes at a distance of about 2.7 Å from the center of the molecule. This number is similar to the distance of 3.4 Å for the CFG B2 structure. The remaining difference can be due to either an overestimated theoretical gas phase gap, or due to the experimental atomic structures having an even stronger binding between molecule and electrodes than CFG B2.

For the structures with small Γ the value of U that brings θ_L in the experimental range is very small, and goes below the expected possible range. Such junctions are therefore expected to exhibit a θ_L well below the experimentally accessible temperatures. This is consistent with the experimental evidence that only a fraction of the molecular junctions, which we attribute to those with the largest Γ, exhibit a Kondo state at an experimentally accessible temperature. Overall our results confirm that the molecule lies flat on a rough Au surface when it exhibits Kondo behavior, since only such structures allow for small binding distances and strong electronic coupling.

3.3. NRG calculations

In order to confirm the trends for θ_L obtained with the simplified model eqn (4), and to evaluate the conductance in the presence of electronic correlations, we integrate NRG calculations in the method. We consider the SIAM representing the PTM/Au structures with the largest Γ (CFGs B2, B4 and T4), for which in the previous subsection have estimated the Kondo temperatures to lie in the experimental range. For each junction we calculate ImΔ(E) in eqn (2) by using DFT + NEGF with the method presented in ref. 24, and the results are shown in Fig. 5. While the value around E_F is similar for all cases, there are pronounced differences in the energy dependence. The NRG calculations then allow us to verify whether the approximation of a constant Γ used so far is applicable for these realistic atomic structures. The real part is obtained with the Kramers–Kronig relation, eqn (3). Further details about the NRG calculations are presented in the ESI section S4.† The many-body self-energy calculated with NRG is then used to evaluate the zero-bias and zero-temperature transmission in the presence of strong correlations in the next section.

Fig. 5 Negative imaginary part of the hybridization function calculated using the DFT + NEGF method, and used as input for the NRG calculation, for the three configurations with the largest hybridizations (B2, B4, T4) (see Fig. 2). High energy contributions are truncated as outlined in the ESI section S4.†

The Kondo temperature is extracted from the impurity contribution to the magnetic susceptibility χ_s(θ)¹⁹ (see ESI section S5†). In Fig. 6 we present χ_s(θ) for the B4 geometry, where the inset shows the small deviation of the impurity occupation n_d from the half-filled case (n_d = 1). We find that χ_s(θ) always follows the same universal behavior as long as the interaction strength is large enough (U > 0.5 eV). A crossover is observed from the high θ local moment regime, where k_Bθχ_s/(gμ_B)² = 1/4, to the low θ strong correlation limit, where k_Bθχ_s/(gμ_B)² = 0; here g is the Landè-factor and μ_B the Bohr magneton.¹⁹ We find that for U values above U = 0.5 eV the curves can be collapsed onto a single universal function. Note that for U = 0.5 eV one can already recognize the deviation from the universal behavior as a dip in the high θ susceptibility.

Fig. 6 Main graph: The universal function displaying the temperature dependent scaling function (gμ_B)²F(θ/c_Wθ_L) = k_Bθχ_s against log(θ/c_Wθ_L). Inset: The impurity occupation as a function of the rescaled on-site energy ε_d + U/2. The on-site energy obtained within DFT is indicated as the vertical dashed line.

The collapse of the susceptibilities is interpreted as a universality due to the formation of a Kondo-singlet. In the local moment regime the static spin-susceptibility scaled by the Kondo-temperature follows the same universal curve,²⁰ where the scaling function F(x) is defined by⁵¹

and where c_W is the so called Wilson number, which is a model-dependent constant (see ESI section S3†). Here, the Kondo temperature θ_L plays the role of a scale invariant in the renormalization group (RG) language. This means that systems with different initial parameters end up in the same low temperature fixed point after mode elimination(RG-flow towards the same fixed point).⁵¹ This gives rise to the universal behavior in Fig. 6 at low θ. The value for θ_L is obtained in the standard way from the condition that the universal function at θ = c_Wθ_L is F(1) = 0.07.⁴³ In Fig. 4 the values of θ_L calculated in this way are displayed as dashed lines. Importantly, they agree rather well with those obtained using the approximate eqn (4), showing that the approximation of a constant Γ is valid for this system. The NRG results therefore also confirm the conclusion that for the three structures with large Γ the value of θ_L is in the experimental range for U ≈ 1 eV.

4. Electron transmission

To evaluate the transport properties of this system we add the zero-temperature NRG self-energy, Σ(E, θ = 0), to the DFT + NEGF Green's function via the Dyson equation and compute the resulting energy-dependent transmission function, T_t(E, θ = 0), in the presence of many-body correlations not captured at the standard DFT-KS level.^24,29,30 As outlined in the ESI sections S7 and S8,† the linear response zero-temperature conductance, G₀ = G(V = 0, θ = 0) = dI(V, θ = 0)/dV|_{V = 0}, is given by:where e is the electron charge, h the Planck constant and 2e²/h the quantum of conductance. Note that we have also implicitly assumed that there is no external magnetic field, whose effect will be considered in the next section. When Γ_L ≫ Γ_R (Γ_L ≪ Γ_R) this can be extended to finite V as G(V, 0) ≈ (2e²/h)T_t(−eV, 0) (G(V, 0) ≈ (2e²/h)T_t(+eV, 0)). As discussed in the previous sections, we expect the Au/PTM/Au system exhibiting Kondo behavior to have such a highly asymmetric coupling. This condition is indeed fulfilled for CFGs B2 and B4, and to a minor extent also for CFG T4, so that the energy dependence of the transmission approximately corresponds to the voltage dependence of the conductance. Note that for such highly asymmetric coupling the dominant effect of the voltage is a shift of the molecular energies due to its induced local electric field, while for the case of approximately symmetric coupling (Γ_L ≈ Γ_R) the current induced non-equilibrium change of occupation gives an additional important contribution and therefore needs to be taken into account. In the next section we will therefore generalize these relations and provide the non-equilibrium relations for the conductance that are also valid for arbitrary values of Γ_L and Γ_R.

As outlined in ref. 24 and in the ESI section S7,† the total transmission function is the sum of the elastic transmission, T, and of the inelastic impurity transmission, T_R,AI, so that T_t = T + T_R,AI. The elastic transmission has contributions from electrons flowing through the impurity, T_AI, from the background transmission, T_B, and from interference terms, T_I (T = T_AI + T_B + T_I). Notably, at zero-temperature, for a system in the Kondo regime one has T_R,AI(E_F) = 0, because the imaginary part of the impurity many-body self-energy vanishes at E_F in accordance with the Fermi-liquid picture.⁵¹

The calculated low energy transmissions for the B2, B4 and T4 configurations are presented in Fig. 7. Here U is set to 1 eV, since this is the charging energy that provides a Kondo temperature in the experimental range. The results for different values of U are shown in the ESI section S7.† One can clearly identify the Kondo peak around E_F, which has a width of the order of 1 meV, in good agreement with the experiments.¹⁶ The overall dominant contribution comes from T_AI for all cases. While the background transmission and interference terms are negligible in the highly asymmetric setups (CFGs B2 and B4), they are rather large in the break-junction geometry T4. Importantly, while in the B2 and B4 geometry the transmission values are very small, for the T4 break-junction configuration they can reach values up to 0.8, and such variations are indeed found in experiments.¹⁶ In the present case the magnitude of both the background transmission and of the Kondo peak become large for symmetric coupling (Γ_L ≈ Γ_R), while they progressively decrease as the coupling becomes more asymmetric. However, we point out that the background transmission may generally be very large if the overlap between the Au electrodes is large or if the electrodes are very broad. In that case one may still have Γ_L ≫ Γ_R for the molecule itself, but the background current will be much larger than that flowing through the molecule. Therefore, for a comparison between theory and experiments for the Kondo conductance itself ideally one needs to separate out the background conductance. While this is difficult to do in experiment, our simulation scheme allows to perform this separation for each atomic configuration. In Fig. 7 we also plot the incoherent transmission T_R,AI and T_t = T + T_R,AI, which determines the measured conductance. As stated above, T_R,AI vanishes at E_F, while it leads to a further overall enhancement of the transmission spectrum away from it. It therefore does not affect the zero-bias and zero-temperature conductance, but it plays an important role at finite bias and finite temperatures, as discussed below.

Fig. 7 Zero-bias transmission including the zero-temperature NRG self-energy for the B2, B4 and T4 structures, and for U = 1.0 eV. Here T is the total coherent transmission, T_t is the total transmission including incoherent effects, T_AI is the coherent transmission component of the AI itself, T_B is the coherent background transmission, T_I is the interference term, and T_R,AI is the incoherent transmission. The total transmission is then T_t = T + T_R,AI, with T = T_AI + T_B + T_I. Note the different scales of the transmission-axis for B2, B4 and T4.

Although the results shown so far are obtained for zero temperature, we can obtain an estimate of the temperature dependence of the full width at half maximum (FWHM, W) of the Kondo peak in the DOS by performing a low energy expansion of the SIAM DOS. For the system investigated here we consider the half-filled particle-hole symmetric case, and moreover, since Γ ≪ U, we are in the so-called strong correlation regime.⁵⁶ As shown in the ESI section S6,† in such a regime the dependence of the FWHM on temperature for a SIAM with energy-independent hybridization Δ = Γ/2 is approximately given by:

Here Δ˜ is the renormalized quasi-particle spectral width, Δ˜ = zΔ, and is the so called wave-function renormalization factor.^22,56 Note that here we use the zero temperature limit of Σ(E, θ), since we perform the perturbation expansion around θ = 0, but in general z can also be evaluated at finite temperature by using the finite-temperature Σ(E, θ) in its definition above. Furthermore, in the particle-hole symmetric regime Δ˜ is related to the Kondo temperature as⁵⁶

Note that the relation in eqn (7) is different from the widely used form given in ref. 57, since in that reference the energy dependence of the real part of the many-body self-energy is neglected. In the ESI section S6† we show that W(θ, Δ˜) from eqn (7) reproduces rather well the NRG results up to temperatures of about 2θ_L.

In experiments, θ_L can be obtained by fitting eqn (7) and (8) to the measured temperature dependent data for the FWHM of the conductance peak. Note that this is somewhat larger than the FWHM of the DOS due to the additional temperature induced broadening of the Fermi distribution of the electrons (see ESI section S8†).

From our NRG calculations we can extract effective values of Δ˜ for the three configurations B2, B4 and T4, which take into account the energy-dependent hybridization at an average level (see ESI section S6†). The resulting values, together with the corresponding Kondo temperatures, are shown in Table 1. Note that the values for θ_L calculated in this way are in good agreement with the values calculated directly from the NRG susceptibility (Fig. 4). In Fig. 8(a) we present the resulting temperature dependent FWHM for all three systems calculated using eqn (7) and the parameters in Table 1.

Fig. 8 (a) Full width at half maximum of the Kondo peak in the DOS, W, calculated with the model eqn (7); (b) normalized impurity conductance as function of temperature using the function in eqn (10), for the B2, B4 and T4 configurations. The results are compared to experimental data from ref. 16, denoted as “Exp₁” and “Exp₂”.

Table 1 Renormalized quasi-particle spectral width, Δ˜, and corresponding Kondo temperature θ_L calculated with eqn (8), as well as wave-function renormalization factor, z, for the configurations B2, B4, T4 (note that Δ˜ and z given here are denoted as Δ˜_∑ and z_∑ in the ESI section S6 and in Table S1). For U = 1 eV we also give the value of U/πΔ, with the values of Δ = Γ/2 taken from Fig. 2

	B2	B4	T4
θ L (K)	1.91	1.96	4.09
Δ˜ (meV)	0.210	0.215	0.449
z	0.00315	0.00326	0.00598
U/πΔ	5.63	5.44	5.05

---

Finally, within the approximation considered in this section we also estimate the temperature dependence of the normalized conductance of the Anderson impurity at zero-bias. If one neglects the interference terms (T_I ≈ 0), then one can write G(V, θ) ≈ G_AI(V, θ) + G_B(V, θ), where G_B is the background conductance originating from T_B, and G_AI is the conductance due to T_AI + T_R,AI. The temperature dependence of G_B is usually small, and for small V also the voltage dependence can be neglected, so that we set G_B to be a constant background conductance. Within the approximations used in this section, the temperature dependence of G_AI is derived in the ESI section S8 (eqn (S29) of the ESI†) to be

where G₀ = G_AI(0, 0). If accurate experimental data are available at low θ, then the mapping of the measured temperature dependent conductance profile to this equation allows to determine the experimental θ_L. However, in many experiments including also those for Au/PTM/Au junctions in ref. 16, the low temperature conductance data is too noisy, so that θ_L is estimated from the high temperature data. Since no analytic expression is available for the whole temperature range, in ref. 58 a functional form is introduced in order to fit calculated NRG results in ref. 59. The proposed fitting curve is:where [small theta, Greek, tilde]_K = (2^1/s − 1)^−1/2θ_K, and θ_K and s are phenomenological parameters. The value of θ_K sets the temperature at which the conductance is reduced by a factor 2 (G_AI(0, θ_K) = G_AI(0, 0)/2). The second order expansion of this relation leads to . As outlined in the ESI section S8,† for the particle-hole symmetric SIAM one can approximate θ_K ≈ θ_L. Furthermore, the condition that the second order expansion needs to be equal to the form given in eqn (9) then sets the value of s to be s ≈ 0.20.

A comparison of the temperature dependent conductance obtained using eqn (10) for the B2, B4 and T4 structures with the experimental data in ref. 16 is plotted in Fig. 8(b). The experimental normalized conductance agrees rather well with the calculated curves, in particular with the one for T4, which has the highest Kondo temperature of all the calculated structures. We denote as “Exp₁” and “Exp₂” the data for the two sets of experiments presented in Fig. 3c and d of ref. 16, respectively. When extracting the experimental Anderson Impurity conductance one has to first subtract the background conductance, G_B. Our calculations show that the background conductance depends significantly on the detailed atomic structure, as shown by the values of T_B in Fig. 7. However, in experiments only the total conductance is accessible. One can approximate the background conductance by the conductance at zero bias for a very large applied magnetic field, which can be extracted from Fig. 3g–h of ref. 16. In this way we extract the ratio of background conductance to the total conductance at zero bias and zero temperature to be about 0.29 for Exp₁, and 0.34 for Exp₂.

While the results presented in this section show good agreement with the experiments in ref. 16, the limitation is that the equations are all based on the assumption of a particle-hole symmetric system, which is not generally the case. Indeed, in ref. 16 it is also shown that by applying a gate voltage the occupation of the SIAM can be systematically changed. At particle-hole symmetry the system is characterized by a single energy scale, k_Bθ_L, and eqn (7)–(10) reflect this property. Away from particle-hole symmetry, however, this no longer holds and corrections to these formulas enter. Furthermore, the condition that Γ_L is very different from Γ_R does not apply for a general system. In the next section we will therefore extend the method to the general non-equilibrium case, and also to the case away from particle-hole symmetry within a perturbative approach.

5. Non-equilibrium relations: renormalized super-perturbation theory

In this section we account for finite-temperature (θ > 0) and general finite-bias (V ≠ 0) effects by using the renormalized super perturbation theory (rSPT) described in ref. 22, 23 and 60. The rSPT corresponds to a perturbative method organized around the particle-hole symmetric strong coupling fixed point considered in the previous sections. While for the PTM/Au system considered here we always have U ≫ Δ, the rSPT relations are in principle valid for arbitrary values of U, and account for deviations from the particle-hole symmetry at a perturbative level. It is based on the insight that at the strong-coupling fixed point the equations have the form of an Anderson model, albeit with renormalized parameters.⁵⁶ These parameters are the renormalized hybridization, Δ˜, which has been introduced in the previous section (Δ˜ = zΔ), the renormalized energy level, [small epsilon, Greek, tilde]_d, which is given by [small epsilon, Greek, tilde]_d = (ε_d + U/2)/Δ, and the renormalized interaction energy, Ũ, defined in the ESI section S9.† We introduce the rescaled renormalized interaction ũ = Ũ/πΔ˜, which lies in the range from 0 for small U to 1 for very large U (see Fig. S9†). In this section we present results as function of [small epsilon, Greek, tilde]_d, which determines the deviation from the particle-hole symmetric case, and which can be tuned experimentally by applying a gate voltage.¹⁶

The Kondo temperature θ_L near the strong coupling fixed point is obtained as ,⁶⁰ with the θ = 0 limit of the static spin susceptibility⁵⁶

and where Ã_AI(E = 0, θ = 0) = z⁻¹A_AI(E = 0, θ = 0) denotes the equilibrium quasi-particle renormalized spectral density at the Fermi energy. Note that for the particle-hole symmetric reference system this definition of θ_L is equivalent to the one presented in section 3 (see also ESI section S3†). Up to second order in ũ[small epsilon, Greek, tilde]_d we have

Ã_AI(0, 0) ≈ [πΔ˜(1 + (1 − ũ)²[small epsilon, Greek, tilde]_d²)]⁻¹. (12)

Inserting this into eqn (11) yields the Kondo temperature

which is a generalization to finite [small epsilon, Greek, tilde]_d and to arbitrary U of the result for the symmetric SIAM in the strong coupling limit given in eqn (8).

A central issue is the relation between renormalized and bare parameters, which is encoded in the wave-function renormalization factor z. The renormalization factor z = Δ˜/Δ can be obtained from NRG for a general energy-dependent hybridization function, and from Bethe ansatz for the case of a constant energy-independent hybridization function.⁵⁶ A comparison between z calculated for a constant hybridization function Δ(E) = Δ(E_F) = Γ/2 using Bethe ansatz and the NRG is shown as function of the interaction energy in Fig. 9, and demonstrates that they agree well. To address the question of the effect of an energy-dependent hybridization function on z, we also calculate z using NRG and with the full energy-dependent Δ(E) for the B2, B4, and T4 structures (Fig. 5). Importantly, we find that they also agree rather well with the results for constant hybridization, showing that the low energy SIAM is largely dominated by the hybridization function around the Fermi energy. Based on these results we therefore calculate z and ũ for the rSPT expansion using the Bethe ansatz for the particle-hole symmetric SIAM with constant hybridization Δ(E_F) = Γ/2 (see Fig. S9 in the ESI section S9†), where the values of Γ for each configuration are given in Fig. 2.

Fig. 9 (Main graph) Comparison between the wave-function renormalization factor, z, calculated using NRG, for a constant hybridization function (black diamonds), and for the three configurations B2 (turquoise filled disks), B4 (pink open triangles), and T4 (white open rectangles). The Bethe ansatz solution for the same constant hybridization used in the NRG calculation is shown as the red solid line. Inset: Wave-function renormalization factor from the main graph plotted on a semi-logarithmic scale.

In order to generalize the relations for the conductance, we first evaluate the equilibrium conductance, G₀ = G_AI(V = 0, θ = 0, B = 0), defined in eqn (6) and (S18) of the ESI,† away from the particle-hole symmetry. Here we have explicitly noted that we consider the reference case with zero magnetic field (B). This results to

Then the extension of rSPT to current-carrying steady states allows us to evaluate the non-linear low-voltage conductance for finite temperatures and also magnetic fields, which has the form:^22,23,60

This result can be obtained by expanding G_AI(θ, V, B) up to second order in eV/Δ˜, k_Bθ/Δ˜, and gμ_BB/Δ˜. The relations for the expansion coefficients are presented in the ESI section S9,† and extend the second order coefficients in U given in ref. 22 to arbitrarily large values of U. Note that the equilibrium transmission calculations, presented in the previous section and in the ESI section S8,† allow to extract the values of c_θ = π² and also c_V = 3/2 in the strong coupling limit (ũ = 1) and at particle-hole symmetry (ε˜_d = 0), and for highly asymmetric coupling to the electrodes (eqn (S29) in the ESI†). Using the general rSPT relations given in the ESI section S9† one can see that as long as ũ = 1 and ε˜_d = 0 these values are valid for arbitrary Γ_L and Γ_R, so that they are independent of the level of asymmetry in the electronic coupling to the electrodes. Note that an important advantage of the rSPT approach is that it is not restricted to these limiting cases, and it is valid for arbitrary values of the parameters, which is a consequence of the fact that it is a truly non-equilibrium method.

The rSPT expansion coefficients calculated for the B2, B4, and T4 structures are displayed in Fig. 10 as a function of the local level energy [small epsilon, Greek, tilde]_d = (ε_d + U/2)/Δ. As noted above, in an experiment this can be modified by applying a gate voltage. We use the Bethe ansatz ε˜_d = 0 for the values of U/πΔ give in Table 1, which then result to ũ = 0.99999418 for the B2 structure, ũ = 0.99999088 for B4, and ũ = 0.99997705 for T4. These values are all very close to 1, and indeed replacing them with 1 leads to essentially the same results, confirming that the Au-PTM system is in the strong coupling limit. The coefficients therefore differ only due to the changes in Γ_L/Γ_R, for which we use the DFT values given in Fig. 2. Since c_θ and c_B are linear-response properties and do not depend on Γ_L/Γ_R, they are identical for all configurations. Consequently, c_θ and c_B can also be calculated via NRG. A comparison for these two quantities between rSPT and NRG is given in ref. 23, where a rather good agreement is found up to moderate values of [small epsilon, Greek, tilde]_d.

Fig. 10 The dependence of the conductance coefficients in eqn (15) on the deviation from particle-hole symmetry, determined by [small epsilon, Greek, tilde]_d = (ε_d + U/2)/Δ(E_F), for the configurations B2, B4 and T4. The mathematical relations for the coefficients are given in the ESI section S9,† and the parameter ζ = 3(Γ_L/Γ_R)/(1 + Γ_L/Γ_R)² in those equations, which determines the asymmetry of the electronic coupling to the left and right electrodes, follows from the values of Γ_R and Γ_L in Fig. 2 as ζ_B2 = 2.6 × 10^–4, ζ_B4 = 2.6 × 10^–4, and ζ_T4 = 0.568. The dimensionless Coulomb repulsion U/(πΔ(E_F)) for the effective Anderson model applicable to each system is presented in Table 1.

The effect of the contact asymmetry, as captured by the ratio Γ_L/Γ_R, affects the value of the finite voltage coefficients, as clearly seen in the lower part of Fig. 10. At particle hole symmetry (ε˜_d = 0) the influence of the contact asymmetry vanishes, except for c_θV. A more detailed analysis of this effect is presented in Fig. 11, where we show c_V as function of Γ_L/Γ_R for different value of [small epsilon, Greek, tilde]_d and U. It can be seen that the overall variations of c_V are rather large, and only as the system goes into the strongly interacting regime (large U) the effect of contact asymmetry becomes small, and it completely vanishes for very large U and ε˜_d = 0, where it reaches the limiting value of 3/2 discussed above. Note that around Γ_L/Γ_R = 1 (symmetric coupling) c_V varies quadratically for small variations of Γ_L/Γ_R around 1 (see also ESI section S9†).

Fig. 11 The dependence of the transport coefficient c_V, obtained from the rSPT, is displayed as a function of the asymmetry in the contacts Γ_L/Γ_R, for different values of the dimensionless Coulomb repulsion U/(πΔ) and local level energy (ε_d + U/2)/Δ.

The rSPT provides a consistent description of the low-temperature, low-field, low-bias transport properties of the Anderson model. When the parameters are calculated from DFT + NEGF, and combined with NRG and/or Bethe ansatz, the method allows for an effectively first principles calculation of all the transport parameters. If the atomic structure is well defined, as is the case in STM experiment of molecules or other adsorbates on flat surfaces,²⁴ the approach is predictive on a quantitative level. When the structure is not known, as is the case for the PTM/Au system considered here, the approach allows to estimate ranges of possible electronic coupling coefficients, interactions energies and deviations from the particle-hole asymmetry. In this case the results give a qualitative guidance to experiments as to which atomic structures are expected to lead to Kondo physics in a measurement.

6. Conclusions

The theoretical modeling of Kondo physics in nanoscale devices is usually limited to fitting the parameters of a SIAM to conductivity measurements. Due to this adjustment of the parameters to the experiment such an approach is therefore not predictive, and the question whether it captures the right physics for a given experiment is therefore open. Moreover, it does not provide any information on the relationship between the device structure and its conductance as well as its electronic interactions. In order to overcome this limitation and provide a predictive model here we present a scheme that obtains the required parameters of the SIAM from DFT calculations for realistic atomic structures. Importantly, conductance measurements are inherently a non-equilibrium process, and our novel scheme combining DFT, NEGF, NRG and rSPT is designed to capture such effects. We derive the equations that relate the equilibrium density of states to the non-equilibrium conductance versus voltage curves, which is necessary to interpret experimental conductance measurements in terms of the electronic and atomic structure of the system. With this approach it is therefore possible to calculate the electronic and non-equilibrium transport properties of strongly correlated molecular junctions in a systematic and predictive way effectively from first principles.

We employ the method for the description of the recently measured Au/PTM/Au break-junctions. The main limitation of break-junction experiments is that the statistical nature of the measurements does not allow a direct understanding of the atomic structures responsible for the conductance and its variations. First-principles calculations are therefore essential to gain a full atomistic insight on the system properties. While state of the art DFT + NEGF can only be applied to weakly correlated systems, the method presented here is proven to overcome this limitation. In fact, for the Au/PTM/Au break-junction we show how the molecule–electrode contacts affect the energy level alignment, charge transfer, hybridization and, ultimately, the Kondo temperature and conductance. Importantly, we show that while the Kondo temperature depends only on the total hybridization of the molecules with the electrodes, the experimental conductance depends also on the relative coupling to left and right electrodes, since those determine the current flow. Our projection scheme allows us to obtain these required individual electronic couplings from DFT, and with these we are able to evaluate the low bias conductance versus voltage curves by means of the rSPT. For PTM molecules weakly coupled to the electrodes, as is the case for a molecule on an idealized perfectly flat Au surface, we predict the Kondo temperature to lie below the experimentally accessible limit. In contrast, for asymmetric junctions with molecules on a corrugated Au surface, where the central carbon atom has a good electronic contact with the Au, the calculated Kondo temperature is in good agreement with experiments. These results are consistent with the experimental finding, where only a limited number of junctions exhibit Kondo features in the conductance at the accessible low temperatures.

Finally, we note that for experimental setups, where the atomic structure is well characterized, such as for certain adsorbates or defects on flat metal surfaces, the method will enable quantitative comparisons with low-noise experiments. By eliminating free parameters it can therefore lead to a systematic understanding of the non-equilibrium Kondo physics of molecular systems. The inclusion of the rSPT allows to predict systematic changes in non-linear transport at low voltage, temperature and magnetic field, which cannot be addressed directly from state of the art calculations of the transmission coefficient alone. Such changes can be induced experimentally, for example by varying the scanning tip height, which modifies the asymmetry in the electronic coupling to the electrodes, and these can then be calculated effectively from first principles with the approach presented here. Our method therefore paves the way toward the rational design of Kondo systems, and the possibility of performing systematic comparisons with unprecedented accuracy between theory and experiments.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

A. D. and I. R. would like to thank E. Burzurí and A. V. Rudnev for useful discussions about transport experiments in Au/PTM/Au junctions, and C. Weber for insightful discussions on the Kondo physics. A. D. and I. R. acknowledge the financial support from the EU project ACMOL (FET Young Explorers, No. 618082). A. D. received additional support from EU Marie Sklodowska-Curie project SPINMAN (No. SEP-210189940) and from the Ministerio de Economía y Competitividad de España (No. FPDI-2013-16641). I. R. acknowledges additional financial support from the EU H2020 programme PETMEM project (Grant No. 688282). I. R. thanks the Cambridge CSD3 HPC centre for providing part of the computing resources. W. H. A., L. C. and D. V. acknowledge the financial support from the Deutsche Forschungsgemeinschaft through TRR80/F6, TRR80/G7 and the FOR1346/P3. E. Muñoz acknowledges financial support by Fondecyt (Chile) No. 1141146. M. M. Radonjić acknowledges the support from Ministry of Education, Science, and Technological Development of the Republic of Serbia under project ON171017. S. Kirchner acknowledges support by the National Key R&D Program of the MOST of China, grant No. 2016YFA0300202, the National Science Foundation of China, grant No. 11774307 and No. 11474250, and the U.S. Army RDECOM – Atlantic Grant No. W911NF-17-1-0108.

References

1. S. V. Aradhya and L. Venkataraman, Nat. Nanotechnol., 2013, 8, 399.
2. L. Bogani and W. Wernsdorfer, Nat. Mater., 2008, 7, 179.
3. D. Xiang, X. Wang, C. Jia, T. Lee and X. Guo, Chem. Rev., 2016, 116, 4318.
4. W. Kohn, Rev. Mod. Phys., 1999, 71, 1253.
5. J. Taylor, H. Guo and J. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 63, 245407.
6. M. Brandbyge, J.-L. Mozos, P. Ordejón, J. Taylor and K. Stokbro, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 165401.
7. A. R. Rocha, V. M. García-Suárez, S. Bailey, C. Lambert, J. Ferrer and S. Sanvito, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 085414.
8. J. J. Palacios, A. J. Pérez-Jiménez, E. Louis, E. SanFabián and J. A. Vergés, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 66, 035322.
9. A. Pecchia and A. di Carlo, Rep. Prog. Phys., 2004, 67, 1497.
10. G. Stefanucci and S. Kurth, Phys. Rev. Lett., 2011, 107, 216401.
11. P. Tröster, P. Schmitteckert and F. Evers, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 115409.
12. J. P. Bergfield, Z.-F. Liu, K. Burke and C. A. Stafford, Phys. Rev. Lett., 2012, 108, 066801.
13. G. Stefanucci and S. Kurth, Nano Lett., 2015, 15, 8020–8025.
14. S. Kurth and G. Stefanucci, Phys. Rev. B: Condens. Matter Mater. Phys., 2016, 94, 241103.
15. D. Jacob and S. Kurth, Nano Lett., 2018, 18, 2086–2090.
16. R. Frisenda, R. Gaudenzi, C. Franco, M. Mas-Torrent, C. Rovira, J. Veciana, I. Alcon, S. T. Bromley, E. Burzurí and H. S. Van der Zant, Nano Lett., 2015, 15, 3109–3114.
17. S. Datta, Electronic Transport in Mesoscopic Systems, Cambridge University Press, Cambridge, UK, 1995.
18. K. G. Wilson, Rev. Mod. Phys., 1975, 47, 773–840.
19. H. R. Krishna-murthy, J. W. Wilkins and K. G. Wilson, Phys. Rev. B: Condens. Matter Mater. Phys., 1980, 21, 1003–1043.
20. H. R. Krishna-murthy, J. W. Wilkins and K. G. Wilson, Phys. Rev. B: Condens. Matter Mater. Phys., 1980, 21, 1044–1083.
21. R. Bulla, T. A. Costi and T. Pruschke, Rev. Mod. Phys., 2008, 80, 395–450.
22. E. Muñoz, C. J. Bolech and S. Kirchner, Phys. Rev. Lett., 2013, 110, 016601.
23. E. Muñoz, F. Zamani, L. Merker, T. A. Costi and S. Kirchner, J. Phys.: Conf. Ser., 2017, 807, 092001.
24. A. Droghetti and I. Rungger, Phys. Rev. B: Condens. Matter Mater. Phys., 2017, 95, 085131.
25. D. Jacob, K. Haule and G. Kotliar, Phys. Rev. Lett., 2009, 103, 016803.
26. D. Jacob, K. Haule and G. Kotliar, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 195115.
27. M. Karolak, D. Jacob and A. I. Lichtenstein, Phys. Rev. Lett., 2011, 107, 146604.
28. D. Jacob, M. Soriano and J. J. Palacios, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 134417.
29. D. Jacob, J. Phys.: Condens. Matter, 2015, 27, 245606.
30. L. Chioncel, C. Morari, A. Östlin, W. H. Appelt, A. Droghetti, M. M. Radonjić, I. Rungger, L. Vitos, U. Eckern and A. V. Postnikov, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 054431.
31. R. Requist, S. Modesti, P. P. Baruselli, A. Smogunov, M. Fabrizio and E. Tosatti, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 69.
32. P. Lucignano, R. Mazzarello, A. Smogunov, M. Fabrizio and E. Tosatti, Nat. Mater., 2009, 8, 563.
33. P. P. Baruselli, A. Smogunov, M. Fabrizio and E. Tosatti, Phys. Rev. Lett., 2012, 108, 206807.
34. J. Liu, H. Isshiki, K. Katoh, T. Morita, B. K. Breedlove, M. Yamashita and T. Komeda, J. Am. Chem. Soc., 2013, 135, 651.
35. Y.-h. Zhang, S. Kahle, T. Herden, C. Stroh, M. Mayor, U. Schlickum, M. Ternes, P. Wahl and K. Kern, Nat. Commun., 2013, 4, 2110.
36. F. Bejarano, I. J. Olavarria-Contreras, A. Droghetti, I. Rungger, A. Rudnev, D. Gutiérrez, M. Mas-Torrent, J. Veciana, H. S. J. van der Zant, C. Rovira, E. Burzurí and N. Crivillers, J. Am. Chem. Soc., 2018, 140, 1691–1696.
37. W. R. French, C. R. Iacovella, I. Rungger, A. M. Souza, S. Sanvito and P. T. Cummings, Nanoscale, 2013, 5, 3654.
38. W. R. French, C. R. Iacovella, I. Rungger, A. M. Souza, S. Sanvito and P. T. Cummings, J. Phys. Chem. Lett., 2013, 4, 887.
39. N. Crivillers, C. Munuera, M. Mas-Torrent, C. Simao, S. T. Bromley, C. Ocal, C. Rovira and J. Veciana, Adv. Mater., 2009, 21, 1177.
40. N. Crivillers, M. Paradinas, M. Mas-Torrent, S. T. Bromley, C. Rovira, C. Ocal and J. Veciana, Chem. Commun., 2011, 47, 4664.
41. G. Seber, A. V. Rudnev, A. Droghetti, I. Rungger, J. Veciana, M. Mas-Torrent, C. Rovira and N. Crivillers, Chem. – Eur. J., 2017, 23, 1415.
42. A. Droghetti, I. Rungger, C. Das Pemmaraju and S. Sanvito, Phys. Rev. B: Condens. Matter Mater. Phys., 2016, 93, 195208.
43. P. W. Anderson, Phys. Rev., 1961, 124, 41–53.
44. G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet and C. A. Marianetti, Rev. Mod. Phys., 2006, 78, 865–951.
45. A. G. Petukhov, I. I. Mazin, L. Chioncel and A. I. Lichtenstein, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 67, 153106.
46. M. I. Katsnelson, V. Y. Irkhin, L. Chioncel, A. I. Lichtenstein and R. A. de Groot, Rev. Mod. Phys., 2008, 80(315), 315–378.
47. P. Huang and E. A. Carter, Nano Lett., 2006, 6, 1146.
48. F. Libisch, C. Huang and E. A. Carter, Acc. Chem. Res., 2014, 47, 2768.
49. Q. Sun and G. K.-L. Chan, Acc. Chem. Res., 2705, 49, 2705.
50. J. Schrieffer and P. Wolff, Phys. Rev., 1966, 149, 491.
51. A. Hewson, The Kondo Problem to Heavy Fermions, Cambridge Univ. Press, Cambridge, UK, 1993.
52. F. D. M. Haldane, J. Phys. C: Solid State Phys., 1978, 11, 5015.
53. B. Himmetoglu, A. Floris, S. de Gironcoli and M. Cococcioni, Int. J. Quantum Chem., 2014, 114, 14.
54. A. M. Souza, I. Rungger, C. D. Pemmaraju, U. Schwingenschlögl and S. Sanvito, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 165112.
55. A. M. Souza, I. Rungger, R. B. Pontes, A. R. Rocha, A. J. Roque da Silva, U. Schwingenschlögl and S. Sanvito, Nanoscale, 2014, 6, 14495.
56. A. C. Hewson, Phys. Rev. Lett., 1993, 70, 4007.
57. K. Nagaoka, T. Jamneala, M. Grobis and M. F. Crommie, Phys. Rev. Lett., 2002, 88, 077205.
58. D. Goldhaber-Gordon, J. Gores, M. A. Kastner, H. Shtrikman, D. Mahalu and U. Meirav, Phys. Rev. Lett., 1998, 81, 5225.
59. T. A. Costi, A. C. Hewson and V. Zlatić, J. Phys.: Condens. Matter, 1994, 6, 2519.
60. L. Merker, S. Kirchner, E. Muñoz and T. A. Costi, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 165132.

---

Footnotes

† Electronic supplementary information (ESI) available: DETAILS. See DOI: 10.1039/C8NR03991G

‡ These authors contributed equally to this work.

---