Multi-orbital non-crossing approximation from maximally localized Wannier functions: the Kondo signature of copper phthalocyanine on Ag(100)

Density functional calculations of gas phase copper phthalocyanine (CuPc) capture the relevant electronic and geometric properties of the molecule [36]. Briefly, CuPc is a D_4h molecule, figure 1, where the d electrons of the central copper atom are split by the ligand field of the surrounding atoms. The ligands capture ∼2 electrons of the Cu atom rendering the molecule in a magnetic d⁹ configuration [36]. The d manifold is split such that the d_x²−y² orbital is singly occupied (SOMO, the singly occupied molecular orbital). The following empty orbital (LUMO, lowest unoccupied molecular orbital) has π character and double degeneracy, because it constitutes the e_g representation [36] of the point group.

Zoom In Zoom Out Reset image size

We have used the Siesta code [37], relaxing the molecule and first two surface layers to forces below 0.04 eV  Å⁻¹using the geometry of [31]. The super-cell contains five layers of 7 × 7 Ag atoms. Norm-conserving pseudo-potentials [38] and strictly localized DZP basis set [39] are used, including an improved noble-metal surface description [40]. The Kohn–Sham Bloch functions were calculated on a Monkhorst–Pack grid 2 × 2 × 1 and transformed to the basis of maximally localized Wannier functions [41–43] (MLWF). Our method for obtaining MLWF from Siesta is described in [32].

The main reasons for working with a MLWF basis set are: (i) the reduction of the computational problem by projecting the M-dimensional Hilbert space of Kohn–Sham states, where M is the full dimension of the CuPc plus Ag(100) electronic problem, onto a smaller N-dimensional space which faithfully represents [42] the energy bands around the Fermi level taking part in the Kondo effect; (ii) the orthogonality of the MLWF basis set that permits us to use standard multi-orbital NCA [16, 22]; and (iii) the extreme localization of MLWF plus their orthogonality gives us a natural way to partition the problem into impurity and substrate subspaces [32].

Maximally localized Wannier functions have been used in the study of strongly correlated matter (for instance see [44–47]). In some cases, the spread minimization is skipped and the Kohn–Sham electronic structure is projected directly onto trial orbitals. This strategy has been used in the study of transition metal oxides [48]. In our system, we found that the projection onto trial orbitals gives unreliable results. The disentanglement method [42] can be contrasted with a direct selection of certain bands of interest. This is often done in the study of transition metal oxides where bands bear strong orbital character [48, 49]. In the adsorbate problem, the band structure is far more complex, which impedes direct band selection. We conclude that the MLWF are an optimal choice for the problem of a large molecule on a metallic substrate.

For the substrate's MLWF, we choose to describe the s–p bands only. In order to achieve this, the s–p bands of the slab are generated by interstitial Wannier functions, in the same way as the surface state of Cu(111) was achieved in [32]. This is a reasonable selection, as long as Kondo physics is considered, because the d bands start at ∼−3 eV below the Fermi energy. This is away from the relevant region in energy with an energy scale several orders of magnitude larger than the typical Kondo scale. Additionally, d states are rather localized, presenting small coupling with molecular states. For these two reasons d states can be safely omitted in the description of the substrate electronic structure taking place in the Kondo effect. The choice of the interstitial centers for the MLWF is delicate because CuPc on Ag(100) displays a C₄ symmetry [31] and it is crucial to ensure that MLWF do not reduce the symmetry of the problem. An unexpected problem of the obtention of MLWF for an extended surface is that convergence to MLWF is more difficult as the lateral dimensions of the super-cell are increased. The surface had to be repeatedly tested to achieve a good description of its electronic bands within 2 eV of the Fermi energy.

The molecule's MLWF are obtained jointly with the substrate ones by using the disentanglement method [42]. In order to achieve the clear partitioning between molecule and substrate, the initial set of trial functions consists of 5 d orbitals of the central copper atom (see figure 1), 16 s states for C–H bonds, 32 p_z orbitals for every carbon atom, 24 orbitals for the nitrogen atoms and 40 s orbitals for C–C and C–N bonds. This Wannier function description was tested in the two separate systems: the clean Ag(100) surface and the gas phase copper phthalocyanine molecule.

2.2.1. One-particle Hamiltonian.

In [32], we showed how to obtain an Anderson-like Hamiltonian from a Kohn–Sham Hamiltonian in a MLWF basis set. There, the evaluation of the main quantities of an atomic magnetic impurity in a non-magnetic host was described. A special emphasis was put on obtaining the intra-atomic Coulomb energy. In the present section, we extend that study to the case of molecules. The Hamiltonian terms have to be expressed in a suitable way to be able to apply the NCA scheme to account for Kondo physics on molecular orbitals.

Thanks to the localization and orthogonality of the MLWF basis, the Wannier-projected Kohn–Sham Hamiltonian can be organized into four blocks

where H_mol contains matrix elements between molecular Wannier functions, the second block following the diagonal, H_slab, refers to Wannier functions of the silver slab, while the off-diagonal blocks, V', involve couplings between molecule and slab.

The eigenstates of H_mol are the molecular orbitals now obtained for the MLWF transformed Hamiltonian. As we shall see below, CuPc is a magnetic molecule because one electron is kept in copper's d_x²−y² orbital. Furthermore, upon adsorption CuPc captures one more electron in the first π-orbital, the doubly degenerated LUMO. Hence, this open-shell structure, when hybridized with the substrate, gives rise to the electron fluctuations of the Kondo effect. These states are singled out of the problem and we define a special subspace for them, _imp, the 'impurity' subspace. The rest of the states are lumped together into the substrate's subspace, _subs.

Hence, we first apply to (1) a unitary transformation that diagonalizes H_mol and leaves H_slab untouched. In this way, we obtain the molecular orbitals of the molecular part of the Hamiltonian, so that the first block of (1) becomes diagonal. Explicitly,

The unitary matrix T defines the transformation from MLWF of the molecule to molecular orbitals. The second step is to choose the molecular orbitals that play a role in the Kondo effect. These orbitals are placed in the first rows and columns of (2) and define the subspace _imp.

After this rearrangement, the matrix of the Wannier-projected Kohn–Sham Hamiltonian reads

where the first block is diagonal and contains selected eigenenergies of H_mol, the block H_subs contains Hamiltonian matrix elements in _subs and the block V contains couplings between states in _imp and _subs. The molecular orbitals that do not directly intervene in the Kondo effect are included in the substrate. Hence, the full mean-field structure of the DFT calculation is preserved within the inner window of the MLWF transformation which is 2 eV around the Fermi energy in the present calculation.

In the present case, the obtention of an Anderson-like Hamiltonian is facilitated by the small values of the intrinsic U term in LDA. As shown in [32], one can extract the intrinsic U term when writing the full Kohn–Sham Hamiltonian in a Wannier basis set. In that reference, the Co atom impurity had values of U around 1 eV. For the present method to work, one should subtract the intrinsic U term such as is done in the LDA + U method [50, 51] or in recent parameterizations of model Hamiltonians [52]. Anisimov and co-workers show that this intrinsic U term is related to the Hund's rule exchange and, thus, is much smaller than the actual Coulomb U values. This is particularly true in the case of molecules. In the case of CuPc, we have evaluated the intrinsic U-term to be ∼30 meV, and hence negligible in front of the molecular level values and hybridization function.

2.2.2. Ab initio calculation of a hybridization function.

The hybridization function has a fundamental role in the impurity physics [53]. It is of uttermost importance to make accurate evaluations of this function for numerical applications. In order to achieve this, we diagonalize H_subs, (3), obtaining Bloch states |kn〉 and Bloch energies _kn, n is the band index and k is from the discretized Brillouin zone with 50 × 50 × 1 k-points. Convergency on k-points was checked by a four-fold increase of the k-point sampling. The couplings V of (3) are transformed to the |kn〉 basis accordingly; we label them V_kn,m, where m indices states in _imp. After these manipulations, we can calculate

from the one-particle Hamiltonian (3). It is a matrix in _imp.

We emphasize that the Hamiltonian (3) is obtained by projecting the Kohn–Sham electronic structure onto a set of N MLWF. The smaller N-dimensional subspace spanned by MLWF is designed in order to represent the selected energy bands around the Fermi level, with the same level of accuracy as the original bands of the M-dimensional space of Kohn–Sham states. The total number of Wannier functions N is the sum of the MLWF of the molecule and of the substrate. The Wannier description has been tested in separate CuPc and Ag(100) systems. M usually means the highest band output from the ab initio calculation. Interestingly, we found that it is very important to verify the convergence in M. Inclusion of bands with energy of even tens of eV above the Fermi level is vital and preconditions the correct projection onto the Wannier space. As a consequence, the important value of Γ_mm(ω) at the Fermi level turns incorrect if M is too small. In the present case, convergence was attained for N = 705 and M = 3000.

Formally, (3) becomes the one-particle part of the Anderson Hamiltonian on bringing H_subs to a diagonal form. Although the procedure we have just introduced is straightforward and essentially algorithmic, it is not correct in principle, because the on-site energies of molecular orbitals in H_imp contain a certain part of the Coulomb energy that is included in the mean-field-like LDA framework [32]. Furthermore, the inadequacy of Kohn–Sham energies to describe true one-electron energies yields an important uncertainty when evaluating the Kondo temperature. At the same time, the hybridization is strongly dependent on the distance of the molecule to the surface. Small variations of these two quantities (one-electron energies and molecule–surface hybridization) exponentially propagate in the Kondo temperature. Present ab initio methods cannot then predict the Kondo temperature of a molecular system.

In this work, the one-particle Hamiltonian H_imp is replaced by a Hamiltonian with many-body interactions. This subsection discusses the physical grounds on which is designed.

The experimental analysis on the Kondo features in CuPc/Ag(100) [13] indicates that the two-fold degenerate e_g LUMO has the main role in Kondo fluctuations. However, DFT calculations [31] show that the spin in the SOMO orbital is not quenched by charge transfer from the molecule. Hence, the observed Kondo resonance is interpreted [13] as due to a spin S = 1 of the molecule, formed by an electron in LUMO and one electron in SOMO. The LUMO electron is captured from the substrate upon adsorption.

We adopt this interpretation and model the molecule as an impurity with a two-electron S = 1 ground state. Furthermore, the lifetime of the SOMO will be assumed an order of magnitude longer than the LUMO's, as will be confirmed by ab initio calculation in section 3.1.2.

The electron in the SOMO will be represented by a local magnetic moment which interacts via exchange with the hybridized LUMO. The multi-orbital structure of the problem appears in the two-fold degeneracy of the LUMO. The configurations that we will study are then the ones coming from the filling of the LUMO, figure 2.

Zoom In Zoom Out Reset image size

We choose the relevant configurations by estimating their energetic accessibility in the charge fluctuation process. Then, the relevant parameters determining the electronic configuration are: the LUMO on-site energy _L and the charging energy U. The lowest-energy molecular configuration for the adsorbed molecule is singly occupied, as suggested by DFT calculations [31]. We can estimate the free molecule U by evaluating the energies of the neutral (E₀), singly (E_I) and doubly (E_II) negatively charged molecule. Assuming a simple impurity Hamiltonian, then the first affinity is given by

and the second affinity by

From these equations and the total energy calculations for the free molecular species we obtain that U = 2.08 eV. However, U is substantially screened on the metallic surface. A constrained DFT calculation for π-systems on silver surfaces leads to a reduction of U to values between 0.5 and 1.0 eV [54]. In the present case, the CuPc molecule lies at ∼3 Å from the surface, reducing the metal screening. Hence, since _L = −0.35 eV (see section 3.2), which leads to |_L| |< U and the doubly occupied configurations e_g ², e_g ¹¹, figures 2(c) and (d), are higher in energy than the neutral one, figure 2(a). Thus, we will consider e_g ¹–e_g ⁰ fluctuations.

Hence, the impurity Hamiltonian contains the eigenvalues of the two LUMO as given by the diagonalization of the MLWF Hamiltonian (2), plus the exchange interaction I with the spin of the SOMO:

We introduce Hubbard operators which automatically restrict the LUMO occupancy to zero or one. The sums are over spin and orbital degrees of freedom, the latter indexed by a. The second term in the Hamiltonian is the direct exchange interaction involving the spin S₂ of the SOMO and the spin operator of the LUMO expressed through Pauli matrices τ as

If I = 0 the Schrieffer–Wolff projection of the model onto the e_g ¹ manifold leads to a SU(4) Coqblin–Schrieffer model [55]. The evolution of spectral features with I leads to a quantum phase transition of the singlet–triplet impurity problem, that has recently received much attention experimentally [56] and theoretically [29, 30]. Here, we will study three cases: I < 0, I = 0 and I>0.

The above impurity Hamiltonian (7) commutes both with the LUMO occupancy operator, n_a = ∑_σ|aσ〉〈aσ|, and the total spin operator S = S₁ + S₂ showing that it properly describes the molecular properties. Due to the hybridization with the substrate, the molecule will exchange electrons. Here, we assume just fluctuations between the two above configurations, figures 2(a) and (b) plus the fluctuations of the SOMO spin, as included in the definition of .

The full Anderson-like Hamiltonian reads

The hybridization part is written with Hubbard operators which change the state of LUMO from empty to occupied or vice versa. Substrate electrons are described by fermionic operators c_kσ and energies _kσ. For the sake of brevity, all substrate electronic degrees of freedom are encapsulated in the k symbol. The band index will be introduced when necessary.

Let . Following Bickers [18], we write down a resolvent operator of the impurity

as a result of averaging over eigenstates |Ω〉 of the non-interacting substrate with energies E_Ω. The substrate partition function is denoted by Z^subs. Since does not mix different occupancies, it can be written as a block matrix

where refers to one-electron occupancy (i.e. LUMO empty) and acts on two-electron occupancies of the impurity.

These quantities are calculated via self-energies defined by

Following Kuramoto [16], we introduce basis states labeled by α for the configurations with one electron and β for the two-electron configurations. In the non-crossing approximation the self-energies for the fixed occupations are given by all diagrams without crossings of substrate electron lines,

The hybridization vertex V_kσ(α|β) comes when emitting an electron from LUMO to the band state kσ which is accompanied by impurity transition from β to the state with label α. Similarly, V_kσ(β|α) is brought about when annihilating a one-electron state α of the impurity and creating a two-electron state by absorbing a substrate electron. Explicitly,

It is convenient to re-express the fixed-occupation self-energies (also called bosonic and fermionic self-energies in slave-boson approaches [17, 22, 30]) in terms of a hybridization function that is directly related to the level broadening of the impurity configurations

Please, notice that we have not included π or 2π factors as sometimes is done in the literature. In order to evaluate lifetimes a factor 2π will be added to the diagonal terms of the hybridization function (10).

We can now represent the fixed-occupation self-energies in the form

This permits us to efficiently perform all computations using fast Fourier transforms.

Let us write β = (a,S,S^z), where a = 1,2 indicates the two orbitals of LUMO, S the total spin and S^z one of its components. Similarly, α denotes the projection of the spin of SOMO on the z-axis, the only degree of freedom of the one-electron configurations.

In the present case, the substrate is non-magnetic, _kσ and V_kσ,a do not depend on the electron spin, σ, and the expression for the hybridization factorizes into spin and orbital parts (see appendix). This considerably simplifies the equations (11). Introducing

and defining the orbital part of the hybridization function by (see (A.1))

we can write

The resolvent for the two-electron configurations now depends on the total spin and is a matrix in the orbital space. The equations (12) and (14) constitute a self-consistent system.

The main quantity of interest, the Green's function of LUMO, will be calculated from the general time-ordered correlation function of Hubbard operators, defined through the thermal average

Starting from the latter expression, we average over the singly occupied configurations αα' and take the Fourier transform [16] to obtain the Green's function of LUMO for the given spin multiplet

The quantities A_2,aa' and A₁ are the spectral functions of the resolvents and Z_i is the impurity partition function. We follow Kuramoto [16] and use defect propagators in the numerical implementation of (16), see appendix.

The spectral function of LUMO is calculated by tracing over orbital and spin degrees of freedom,

3.1.1. Molecular orbitals.

Thanks to the partitioning made possible by the MLWF basis, we can extract the molecular Hamiltonian H_mol (1) computed from the LDA calculation of CuPc on Ag(100) [31, 60]. Table 1 presents eigenenergies of H_mol the molecular orbitals close to the substrate's Fermi level. The orbitals are then closely related to the gas phase CuPc ones and hence, we denote them by SOMO, LUMO and highest occupied molecular orbital (HOMO) labels. For comparison, the results for the LDA Siesta calculation of the gas phase molecule of the same molecular geometry are also given. The MLWF qualitatively reproduce the gas phase molecular levels. We emphasize that the (MLWF) refer to molecular orbitals screened by the metal but without direct hybridization and the (Siesta) are for a true gas phase molecule without any screening or coupling with an external metal.

Table 1.  Eigenenergies of the molecular orbitals for the impurity Hamiltonian evaluated with MLWF using (2), compared with the gas phase LDA calculation using Siesta, for the same geometry of the CuPc molecule. All energies in the third column were shifted, so that the energies of SOMO (↓) coincide. Energies of the second LUMO state are in parenthesis. The comparison is only indicative that the MLWF capture the molecular properties since the two calculations are performed for different systems: the MLWF refer to molecular orbitals screened by the metal but without direct hybridization and the Siesta column is for a true gas phase molecule without any screening or coupling with an external metal.

	(MLWF) (eV)	(Siesta) (eV)
HOMO (↑)	−1.110	−0.888
HOMO (↓)	−1.112	−0.890
SOMO (↑)	−0.549	−0.714
SOMO (↓)	0.040	0.040
LUMO (↑)	0.185 (0.190)	0.496 (0.498)
LUMO (↓)	0.211 (0.216)	0.520 (0.522)

Figure 3 shows the PDOS onto these four molecular orbitals. The calculation is performed for the LDA k-point grid (2 × 2 × 1) which is clearly insufficient to account for the continuum character of the substrate electronic states. Hence, the peaks have been slightly broadened using a Gaussian broadening of 50 meV. Despite the qualitative character of this figure, much information can be gleaned from it. The difference in widths between SOMO on the one hand and LUMO and HOMO on the other hand is substantial. The SOMO peak basically presents the numerical width, and is a featureless peak, while the LUMO shows oscillations and spreads over several hundreds of meV. Then, the latter orbitals hybridize more strongly with the substrate. From this figure we can conclude that while the LUMO width is in the range of a few hundreds of meV, the SOMO is significantly less. This is in agreement with the move involved calculations of the hybridization function that we present in section 3.1.2.

Zoom In Zoom Out Reset image size

The PDOS are in excellent agreement with the PDOS on molecular orbitals from DFT results [60]. Both the eigenenergies and the PDOS give support to our method of transforming the DFT Hamiltonian to a MLWF and selecting the impurity Hamiltonian using (3).

It is difficult to conclude on the actual occupancies from the PDOS. Indeed, the PDOS numerical broadening thwarts any precise calculation of orbital occupancies, and the definition itself of occupancy of a molecular orbital in a chemisorbed system is somewhat arbitrary. As we show in the next section, the actual peak width of SOMO is negligible; the SOMO then remains singly occupied as in the gas phase. The LUMO peaks in figure 3 cross the Fermi level, hence these orbitals capture charge from the substrate. The LUMO are then the only orbitals that participate in charge transfer from the surface.

3.1.2. Hybridization function.

The calculation of (4) proceeds along the lines described in section 2.2.2. The main ingredients of this calculation are the Bloch eigenvalues and the molecule–substrate couplings. The Bloch eigenstates span the substrate continuum, and this task is much simplified thanks to the reduced computational effort that the Wannier basis set introduces. Indeed, this allows us to use a very dense k-point sampling mesh and hence, accurate results. The couplings are obtained after rotation into the molecular orbital basis set (3) of the couplings obtained after the Wannierization process. This procedure is numerically very reliable because it assures perfect agreement with the LDA results within an energy window of a few eV about the Fermi energy. Hence, the shortcomings of our procedure to compute the hybridization function are the ones inherent to the LDA in the present case. The molecule–substrate couplings depend exponentially on the molecule–surface distance. This means that the determination of the system's geometry is very important to obtain accurate hybridization values. Otherwise, we expect LDA values of the hybridization function to be accurate because the hybridization function evaluates the transparency of the barrier between the molecule and the substrate and good numerical results have been previously obtained  [61].

The off-diagonal parts of Γ_mm'(ω) are very small (<1 meV). This is a very strong result that shows that the substrate does not mix the different molecular orbitals among themselves. Hence, every molecular orbital defines an electronic channel of the system.

The diagonal elements for SOMO and LUMO orbitals are presented in figure 4, along with the density of states of the substrate. The three curves are very similar in shape. By comparing the hybridization functions with the substrate DOS, we notice that the hybridization function grows more slowly than the DOS as the electron energy increases. This is an effect due to the couplings, V, between the molecular orbitals and the substrate electronic structure. However, at higher energy, it is the DOS that controls the behavior of the hybridization function with energy. Figure 4 shows that the SOMO width has to be multiplied by 20 to be comparable to the LUMO one. The SOMO orbital has then a very small mixing with the substrate as compared to the LUMO. Finally, in the U = 0 picture, the FWHM for the LUMO is 2πΓ_aa(ω = _L) = 440 meV and is comparable to the FWHM of the PDOS peak, figure 3, obtained with an insufficient k-point sampling.

Zoom In Zoom Out Reset image size

We have also evaluated the hybridization function of the LUMO with d orbitals by preparing a Wannier basis set with d orbitals. In the region of interest here (some 2 eV about the Fermi energy), this hybridization is strictly zero due to the lack of d states at these energies. However, when resonant with Ag d-band (3 eV below the Fermi energy), the hybridization function becomes larger than the corresponding values for the sp Wannier functions. At 4 eV below the Fermi energy, the d-electron hybridization has a maximum of 0.032 eV, while the sp one is at 0.035 eV. The consequence of this is that the spectral function at −4 eV will not be just a simple Lorentzian tail. However, the LUMO orbital is some hundreds of meV away from the Fermi energy, and the effect of the d-band contribution to the overall shape of the LUMO spectral function is negligible both for Kondo physics and for the one-electron spectral shape.

Finally, we comment on the problem of broken spin symmetry in DFT. Spin-polarized LDA implies two subsystems: minority and majority spin. This in turn says that we have two sets of Wannier functions, two distinct substrates, hybridization functions, etc. The main effect of breaking the spin invariance in our DFT calculation is that the SOMO occupancy is n_SOMO≈1. As a secondary effect, the substrate and molecule become slightly spin-polarized as well. However, this effect is perturbational [14] and is not an intrinsic property of the bare substrate and impurity of the Anderson model. In what follows we drop the minority spin data of the hybridization function and restore the spin symmetry.

Our LDA calculations yield a hybridization function for the SOMO, Γ_S, ten times smaller than the one for the LUMO levels. From these values we obtain Γ_S≈4.5 meV. We use the standard expression for the Kondo temperature [62]

and take the ab initio value for the on-site energy from the table 1. The bandwidth is 2D, where a rectangular DOS is typically assumed. Here we have taken D as 10 eV because the DOS of figure 3 integrates to 703 states, while the calculated DOS at the Fermi energy is 32.28 eV ⁻¹. Hence 2D = 703/32.28 = 21.8 eV, and D turns out roughly 10 eV. We get T_K,S of the range 10⁻²⁶ eV . Thus, the energy scale that would correspond to a Kondo effect on the SOMO orbital (without exchange coupling to LUMO) is unobservable.

These arguments show that the Kondo coupling in this system is indeed given by virtual charge fluctuations of the two-fold degenerate LUMO. Hence, we are dealing with the impurity problem described in section 2.3, for which we can calculate the spectral function according to the section 2.4.

However, the on-site energies of LUMO in (7), as given by their LDA values in the table 1, lie very close to the Fermi level, which would correspond to a fluctuating valence regime. That has not been observed in experiment [13]. This is related to the fact that the LDA energies of LUMO do not reflect the considerable Coulomb repulsion in their spin splitting. We fix these deficiencies by rescaling the _L/Γ_aa(ω) ratio in order to achieve a SU(4) Kondo temperature, T_K,L = Dexp(−|_L|/4Γ_aa), of ∼20 K as observed in the measurements [13]. The results presented in this section are calculated with _L = −0.35 eV and Γ_aa(_L) rescaled to 0.01 eV.

For temperatures below T_K,L the model Hamiltonian (7) and (8) exhibits a rich variety of physics as the value of I changes. Four regimes can be identified: (i) I>0, I>T_K,L, (ii) I < 0, |I|>T_K,L, (iii) I>0, I < T_K,L and (iv) I < 0, |I| |< T_K,L.

In the first case, taking the limit I→∞, the Kondo physics is that of the underscreened Kondo effect. Hence, the model will show singular Fermi-liquid characteristics [63–65] in the low-energy domain. The physics in (ii) is dominated by the molecular spin zero states. In the limit I→−∞ we obtain a spin-less molecule with orbital pseudo-spin, which is over-screened by two substrate channels. The cross-over temperature is, however, exponentially smaller than T_K,L, due to the reduced degeneracy. In the regime (iii) the SOMO is effectively decoupled from the LUMO and we recover the case of two degenerate LUMO on the same footing as the electron spin: we have an SU(4) system as described above. Finally, in (iv), the negative I allows for the Kondo screening of the SOMO by a two-stage Kondo effect [66] at very low temperatures. Here we present NCA spectral functions in the I = 0 case and in the intermediate regimes Γ_aa>|I|>T_K,L which are dominated by inelastic spin transitions.

Figure 5 shows the spectral function of LUMO at T = 7 K when the exchange interaction between SOMO and LUMO is turned off, I = 0. The ground state of is four-fold degenerate in this case. Since each of the two orbitals of LUMO defines an independent scattering channel in the substrate (i.e. Γ_aa'(ω) is diagonal), the LUMO is subject to a SU(4) Kondo effect. The spectral function shows a broad charge-excitation peak of the FWHM given by 4×2πΓ_LL(_L) and a narrow Kondo peak.

Zoom In Zoom Out Reset image size

When the SOMO and LUMO spins are subject to a ferromagnetic interaction (I>0) the ground state of the molecule without couplings to the substrate is an orbitally degenerate S = 1. The excited state is a S = 0 orbital doublet. Now, two satellites develop at ±I from the Kondo peak, figure 6. These satellites are inelastic replicas of the Kondo peak since the spin-excitation energy is exactly I. Decomposition into spin channels yields that the Stokes ( + I) satellite is in the S = 0 channel, while the anti-Stokes peak as well as the Kondo peak are in the S = 1 channel. This result can be understood by recalling the meaning of the spectral function. The spectral function at T = 0 yields the probability density to inject one electron in the system when ω>0 or to inject a hole when ω < 0. Hence, the positive energy satellite corresponds to an excited Kondo effect triggered by the injection of an electron to the S = 0 state, which is an excited state of . Similar results have been found by Roura Bas and Aligia [29, 30].

Zoom In Zoom Out Reset image size

When the interaction is anti-ferromagnetic (I < 0), the ground state of is an orbital doublet and is followed by six excited states of spin one. A noteworthy feature of the spectral function are the steps typical for an inelastic spin-flip transition (see for example [57, 58] and references therein). The inelastic steps enhanced by Kondo effect have already been explored in the case of singlet–triplet transitions in nanotubes [67]. The spin channel analysis yields that the Stokes peak is now S = 1 which again can be rationalized by noticing that it corresponds to injecting an electron in an excited S = 1 state and the transition is enhanced by an excited Kondo effect. By the same token, the anti-Stokes peak is S = 0 since a hole is created in this channel.

The spectral function, figure 7, shows a small peak on the Fermi level. This peak cannot correspond to a spin-flip Kondo effect. The orbital Kondo resonance is ruled out by its exponentially suppressed energy scale ∝exp(−|_L|/2Γ_aa) as compared to T_K,L. NCA is known to produce spurious peaks at the Fermi level [15, 23, 18]. The temperature scale at which they appear is given by [62] T_S = T_K,L/5[T_K,L/Γ_aa]^5/3 = 0.18 K (strictly valid for I = 0 only). All calculations shown here are performed at 7 K, well above the pathology temperature T_S. We conjecture that the zero energy structure is related to the problems of NCA to reproduce spin-split Kondo peaks due to the artificial flow of the marginal potential scattering term as shown in [22, 68].

Zoom In Zoom Out Reset image size

Experiments with the STM measure the conductance as a function of bias. As a first approximation, the conductance is basically proportional to the substrate spectral function, justifying our using of the spectral function as a first approach to the experiment. However, a more detailed description of the electron transmission process from the STM tip [69] shows that the transmission channels can interfere, giving rise to a complex Fano structure instead of the above peaked spectral function. In [8] the different Fano structures are classified according to the prevalence of one transmission channel above the others. Hence, peaks are expected if the channel involving the adsorbate is preferential in the conductance process, while dips appear in the conductance spectra if the direct contribution of the surface overwhelms the transport process. In the present case [13], peaks are found, suggesting that transport in such a large molecule is molecule-dominated, however recent reports on porphyrins [8] show a dip in the conductance despite the similarity between both systems. Hence, a detailed conductance calculation is necessary to be able to compare with the experimental data.

The couplings between α and β (see (9)) are given by the expression

with V_kσ,a from the Hamiltonian (8) and the bracket is a Clebsch–Gordan coefficient. When the substrate is non-magnetic, _kσ and V_kσ,a do not depend on σ and the expression for hybridization intensity (10) factorizes

into orbital

and spin parts

The subsequent identities have proven significant

The first one is due to completeness (A.2), the second and third can be proven using Wigner 3jm symbols [74].

With these identities it is easy to show that the following property holds: let and be some functions diagonal in the total spin representation, so will the self-energies calculated from NCA (11a) and (11b),

Since we solve the equations of NCA iteratively, starting with resolvents having zero self-energies, we conclude that spin off-diagonal terms of vanish.

Equation (16) has to be expressed in terms of defect propagators. For convenience we introduce boldface notation for matrices in the orbital space of LUMO, i.e. for the resolvent of the two-electron configuration, its spectral density, hybridization function (A.1) and the Green's function of LUMO. The starting expression (16) reads

We introduce the operator , whose effect on an arbitrary matrix function X(ω) is given by

Applying on both sides of NCA equations (14) yields a self-consistent system

for the defect propagators [23], defined by

When these equations are solved, the Green's function (A.6) can be expressed as