Including the probe tip in theoretical models of inelastic scanning tunneling spectroscopy: CO on Cu(100)

G. Teobaldi,1 M. Peñalba,2 A. Arnau,3 N. Lorente,4 and W. A. Hofer1 1Surface Science Research Centre, The University of Liverpool, Liverpoool L69 3BX, United Kingdom 2Departamento de Fisica Aplicada, Escuela Universitaria Politecnica UPV/EHU, 20018 San Sebastian, Spain 3Departamento de Fisica de Materiales y Centro Mixto CSIC-UPV/EHU, 20018 San Sebastian, Spain 4Centro de Investigaciones en Nanociencia y Nanotecnología, CSIC-ICN, Campus UAB, Bellaterra, Barcelona 08913, Spain Received 31 August 2007; revised manuscript received 17 October 2007; published 6 December 2007


I. INTRODUCTION
Since the pioneering work by Jaklevic and Lambe, 1 where the changes in conductance across a metal-metal oxide-metal tunneling junction could be referred to the vibrational modes of molecules in the tunneling barrier, the field of inelastic electron tunneling spectroscopy ͑IETS͒ has considerably matured. 2 More recently, the new modes of operation in scanning tunneling microscopy ͑STM͒ have permitted chemical identification of molecules at metallic surfaces by their vibrational signature, manipulation of individual atoms, bond cleavage, and induced chemical reactions. [3][4][5][6] However, the inherent complexity of IETS makes data interpretation difficult in many cases and, therefore, represents a challenge for theory. The main problem is that even though several thorough theoretical methods have been developed, [7][8][9][10][11] in practice, very few approaches are compatible with state-ofthe-art ab initio methods for many-electron systems.
In this paper, we describe the practical implementation of a Bardeen approach for IETS, compatible with firstprinciples electronic structure simulations of the whole system. The method permits us to include the tip structure explicitly and is, therefore, particularly well suited to describe experimental observations in which the tip structure plays a decisive role. The model system in our study is a CO molecule on Cu͑100͒, which has been widely researched in experiments. 2,12 In Sec. II, we describe the essential details of our first-principles method to treat the electronic structure of tip and sample together with the IETS method itself. To improve readability of the paper, the details of the theoretical formulation are presented in the Appendix. The results and their discussion are reported in Sec. III. We first describe the chemisorption of CO on Cu͑100͒ and the frustrated rotation and C-O stretch vibrational modes. Subsequently, we present STM topography simulations. Finally, we perform the IETS analysis of the two frustrated rotations and the C-O stretch mode. The paper concludes with a short summary.

A. Electronic structure and vibrational modes
CO chemisorption on Cu͑100͒ has been modeled with the density functional theory projected augmented wave ͑DFT-PAW͒ approach as implemented in the VASP program. 13 Exchange and correlation energies were described within the PW91-generalized gradient approximation ͑GGA͒. 14 The Cu͑100͒ surface was described by a periodic supercell consisting of a four-layer slab using the optimized bulk lattice parameter for Cu ͑fcc͒ of 3.65 Å. In order to simulate a molecular coverage low enough to account for single molecule properties, a single molecule of CO has been optimized on top of a Cu͑100͒-͑4 ϫ 4͒ unit cell ͓see Fig. 1͑a͔͒. 15,16 The molecule and the two topmost Cu͑100͒ layers have been fully relaxed until the forces on individual ions were less than 0.01 eV/ Å. The surface Brillouin zone was sampled in all simulations with a grid of 4 ϫ 4 ϫ 1 special k points. The energy cutoff was 400 eV. Vibrational energies and displacements have been simulated in a two-step approach: first, the dynamical matrix was computed on the basis of small finite atomic displacements along the three Cartesian axes, and then the dynamical matrix was diagonalized yielding the ensuing vibrational eigendisplacements and frequencies ͑eigen-values͒. All frequencies have been obtained on the basis of a fixed atomic displacement of 0.02 Å. In analogy with Refs. 15, 21, and 30, the root mean square displacement of each mode has been used for the calculation of the deformation potential as explained in these references. The tip structure has been modeled by a ͑46 atoms͒ periodic 13.53ϫ 9.57 Å 2 W͑110͒ atomically sharp bipyramid ͑see Fig. 1͑b͒ and Refs. 17 and 18 for further details͒. The vacuum range in the calculations was 15 Å. The ensuing inelastic efficiencies were found to be converged within 0.15% with respect to a 20 Å vacuum range, and the maximum level of accuracy for the fast Fourier transform algorithm currently implemented in VASP. 13 STM simulations both from the local density of states of the isolated molecule-surface system 19 and including the STM tip 20 were performed with BSKAN. 17,18

B. Inelastic tunneling spectroscopy simulations
Assuming low temperature ͑harmonic, zero-phonon ap-proximation͒, weak electron-vibration coupling and weak tip-sample interactions, the change in conductance across a vibrational threshold ͑eV = ប, where is the mode frequency, e the electron charge, and V the STM bias voltage͒ can be evaluated from the first order change in the sample wave functions ͑␦ n ͒ due to the perturbing deformation potential ͑U͒ produced by the vibration: 16 Here, m ͑ n ͒ are the unperturbed sample wave functions and ⑀ m ͑⑀ n ͒ the corresponding eigenvalues. The evaluation of the Bardeen matrix elements at the Fermi level ͑⑀ F ͒ requires also the unperturbed tip wave functions and corresponding eigenenergies ͑⑀ ͒. In the quasistatic approximation ͑ប Ӷ ⑀ F ͒, the change in conductance is given by ͑we use atomic units in the following ប = m = e =1͒ ⌬ͩ dI and the surface ͑S͒ integration is carried out in the vacuum region. 20 The strength of the inelastic contribution is given by the modulation of the tunneling amplitude due to the vibration. In addition to the evaluation of Eq. ͑1͒, and substitution in Eq. ͑2͒, one has to account for another contribution to the change in conductance due to the electron-vibration mediated exchange interaction that leads to interference in the elastic channel. In practice, one evaluates two terms, P and I , 21 to calculate the total inelastic efficiency, tot = P − I , defined as the ratio of the total change in conductance ͑due to the crossing of the vibrational threshold͒ to the unperturbed elastic conductance. In the Appendix, we derive explicit expressions for P and I in terms of the overlap integrals.
The proposed method to simulate inelastic tunneling spectroscopy thus requires ͑i͒ the electronic structure of the tip and sample from first-principles DFT calculations, ͑ii͒ the vibrational modes of the sample ͑molecule on a substrate͒, ͑iii͒ the change in the sample wave function due to the electron-vibration perturbation for each mode, described in the Appendix, and ͑iv͒ the inelastic efficiency maps from the elastic conductance and the changes due to the crossing of the vibrational threshold implemented in BSKAN. 17,18

A. Chemisorption of CO on Cu(100) and vibrational modes
In order to model CO chemisorption on Cu͑100͒, we have focused on the experimental on-top adsorption geometry. The adsorption geometry calculated with the PAW method is in good agreement with previously published pseudopotential data. The deviations between the two approaches are less than 1%. 15 The bond length of C-O is 1.20 Å, and the distance to the surface, C-Cu, is 1.84 Å resulting in a global deviation of roughly 5% from experimental low-energy electron diffraction data. 22 For the adsorption energy, defined as the difference between the total energy of the full optimized system minus the sum of the energies for noninteracting relaxed molecule and surface, we find −0.96 eV.
The vibrations of the adsorbed molecules have been determined on the basis of a reduced dynamical matrix displacing CO and the topmost metallic layer from the equilibrium position. Table I reports the results for selected modes of CO / Cu͑100͒. We labeled the C-O stretch mode by S and the frustrated rotations by R. They correspond to motion of the C and O atoms nearly parallel to the surface but in opposite directions. The two frustrated rotations are quasidegenerate for CO / Cu͑100͒. This aspect can be readily understood in terms of surface symmetry. In fact, for a fcc ͑100͒ surface, the geometrical corrugation is the same for both crystal directions, ͓011͔ and ͓011͔. With deviations from the experimental STM-IETS data 25 of 11 meV and 6͑3͒ meV for the S ͓256 meV ͑Ref. 25͔͒ and R ͓36.3/ 33.2 meV ͑Ref. 25͔͒ modes, respectively, we note that the agreement between gradient corrected vibrational energies and the experimental values is in line with standard DFT accuracy for vibrational calculations ͑error bar=5-10%͒.
The weak chemisorption of CO on Cu surfaces agrees well with the classical Blyholder model of chemisorption: 23 Upon adsorption, charge is transferred from the occupied 5 CO orbital to the metal which, in turn, back donates charge density to the molecule into the previously empty 2 * orbital. In view of modeling vibrational tunneling efficiencies, our main interest is the hybridization with the metal of the pristine twofold degenerate 2 * orbital upon adsorption on Cu͑100͒. The calculated partial densities of state ͑PDOSs͒ for CO 2 * orbitals are shown in Fig. 2. The 2 * orbitals broadens into a partially filled resonancelike feature upon adsorption on the metallic surface. The present results for CO / Cu͑100͒ show a 2 * peak centered at 1.80 eV with a full width at half maximum ͑FWHM͒ of 2.5 eV in reasonable agreement with previous simulations ͑2.9 FWHM, main peak at 2.2 eV͒. 15 Compared to experimental inverse photoemission spectroscopy ͑IPS͒ data, suggesting a 2 * peak centered at 3.6 eV above the Fermi energy for CO-Cu͑100͒-2 ϫ 2, 24 the peak positions are shifted by more than 1 eV. The discrepancy between experiment and theory is most likely due to the known deficiencies of the local density approximation: the energy difference between occupied and unoccupied states of molecules is generally underestimated in LDA and GGA-based electronic structure simulations. An additional contribution could be the injection of electrons into the system in IPS, which should give rise to negative ion states.

B. Topographical images of CO on Cu(100)
STM topographies for CO / Cu͑100͒ were first simulated using the Tersoff-Hamann approximation, where the lowbias tunneling conductance is assumed proportional to the sample's local density of states at the Fermi level. 19 CO is imaged as a protrusion with a maximum corrugation of 1.0-1.2 Å for distances in the range of 5 -9 Å ͓see Fig.  3͑a͔͒. Calculating the constant current contours including an STM tip, we find that the appearance of CO / Cu͑100͒ changes substantially with the distance from the surface. A constant current topography is shown in Fig. 3͑b͒. Close to the surface ͑ϳ6 Å͒, CO is still simulated as a protrusion of 0.7 Å with an FWHM value of about 5 Å. However, at a distance of 9 Å, it appears as a depression of 0.3 Å ͑FWHM of 8 Å͒. The change of appearance from protrusion to depression in the constant current contours is due to a k-parallel selection: While for small distances all parts of the surface Brillouin zone contribute in nearly equal measure, at high distances, the larger decay length of states at the ⌫ point leads to a gradual damping of states near the Brillouin-zone edges, responsible for the appearance as a protrusion in the low distance regime. This can be verified by simulating local density contours with a restricted sampling of k points. In our view, this high dependency on the origin of the electron states may account for the widely varying appearances of CO on Cu͑100͒ found in the literature. [25][26][27] However, for the purpose of inelastic tunneling spectroscopy, this appearance, which finds its expression in a variable shape of the molecule depending on the tip electronic structure and gap resistance, does not alter the obtained results. In this case, the important magnitude is the relative change of the conductance due to the opening of an inelastic pathway: this relative change is independent of the exact appearance of the molecule on a metal surface.

C. Inelastic electron tunneling spectroscopy signal of frustrated rotations
We have applied the theory of vibrational inelastic tunneling described in Sec. II to the twofold degenerate R mode of CO on Cu͑100͒ at the Tersoff-Hamann level. The nonadiabatic coupling for this vibrational mode is found to be very efficient. In Table II, we report the maximum calculated values for the total inelastic efficiencies tot , as well as their contributions P and I for different distances from the surface. To illustrate the influence of the energy range on the simulated results, we have calculated the efficiencies with two different values for the thermal broadening, 0.40 and 0.25 eV, respectively. The adiabatic coupling and, therefore, the tunneling efficiency are very weakly dependent on the absolute value of the local density of states ͑LDOS͒: The change of P is of less than 30% over 2 orders of magnitude change for the LDOS. The I contribution contributes only marginally to the total inelastic efficiency tot , as found previously. 15 This suggests that at the vibrational threshold, the nonadiabatic coupling is mainly diagonal ͓see the Appendix, first term of Eq. ͑A7͒, explicitly derived in Eq. ͑A9͔͒ and that the extent of crossband mixing leading to the elastic  The corresponding profile along ͓001͔ is reported in Fig. 4͑c͒ for several LDOS values. In agreement with experiments, 25 we find that the global signal is symmetric with respect to the adsorption site. With a FWHM spanning 1.5-3.8 Å within the considered range of the LDOS, it is narrower than the constant density contour. 15 As a final remark, we note that the calculated efficiency for the R modes ͑17.   that including the tip does not alter the trend of almost constant calculated efficiencies in a wide range of gap resistances. Interestingly, while the absolute change in I with respect to TH is nearly negligible, this is not the case for the P contributions, which are also significantly lower ͑by a factor of 2͒ compared to the TH efficiencies. These results suggest that a real STM tip will not reflect every possible nonadiabatic coupling, but only the part which is suited to the specific experimental tip geometry and lateral resolution. Other tip structures with different terminations can change these values substantially, resulting in effective selection rules for active modes in IETS. 29 In addition, we investigated tip effects with respect to the lateral extent of the inelastic efficiencies tot , P , and I . The P contribution remains symmetrical around the adsorbed CO molecule, while the double-featured I contributions remain perpendicular to each other and aligned along the ͓011͔ and ͓011͔ directions, i.e., the vibrational R 1 and R 2 displacement directions. We have summed up the total R 1 and R 2 contributions to compare with the TH results above and experimental data. The results are shown in Fig. 4͑b͒. Similar to the previous analysis, the global R efficiency is symmetrical around the CO adsorption site, with a FWHM of about 4 Å. The net final efficiency for the R mode is more localized than the topographical feature associated with CO, as shown in Fig.   4͑d͒. The absolute value of the tot efficiency is now in good agreement with experimental value: the maximum calculated value ranges from 11.3% to 11.8% for the simulated currents, while the corresponding experimental value is ϳ8%. These values indicate that including an explicit description of the STM tip in the theoretical model removes the discrepancy of about 50% for the simulated efficiencies. The overestimation is due to neglecting the implicit limitations of the adiabatic coupling. This error is remedied if an explicit electronic tip structure is included in the description. However, a note of caution is in place. The tip, which determines the tunneling matrix elements and acts as a kind of filter concerning the included coupling, is different for every apex geometry and for every chemical composition of the STM tip. The reduction of the efficiencies obtained at the TH level thus implicitly depends on the tip structure. This also implies that the local extension of the available inelastically coupled modes will have an effect on the variation of the TH values with respect to the Bardeen values. This introduces a second dimension into the filtering problem: if the local extension of the available modes is high, the efficiency due to the filtering will be substantially reduced. If, however, the extension is low to start with, the values will not change as dramatically. It will be seen in our analysis of the stretch mode that this also has an impact on the results obtained with a realistic tip structure. However, at present, this behavior has only been established for single modes ͑S and R 1 , R 2 ͒ on a single system ͓CO on Cu͑100͔͒. It will need considerable efforts simulating the IETS signal of other molecules to verify a general trend.

D. Inelastic electron tunneling spectroscopy signal of C-O stretching
Also, for this mode, the inelastic tunneling efficiencies tot , P , and I have been calculated at both the TH and the Bardeen level. The results are summarized in Table III. In contrast to the rotational modes, the largest P values are much closer to the I values, resulting in a strongly reduced global total inelastic efficiency tot . As for the R case, also for the S mode, the computed efficiencies depend only weakly on the specific value of the considered LDOS and/or tunneling current contour. However, in contrast to the R results, the TH efficiencies are lower than the corresponding Bardeen efficiencies by a factor of about 1.2-1.8. This is due, as already mentioned, to the higher localization of the stretch mode compared to the rotational modes. The tip in this case has practically no effect on the obtained efficiencies. The results are in good agreement with the experimentally detected value of ϳ1.5%. 25 The slight underestimation of the calculated data with respect to the experiment is probably due the fact that in the present case, we have not introduced the long-range dipolar term, which has been previously recognized as an important term for the S mode on Cu͑100͒. 15 Due to the symmetry of the mode displacements and the influence of the 2 * CO orbitals, the appearance of P and I is similar: both contributions are symmetrically displaced along the CO intermolecular axis. Over the simulated LDOS values, the FWHM of the main protrusion changes from 5.6 to 8.8 Å. Thus, at TH level, the calculated lateral extent of the total efficiency is slightly overestimated and evaluated to be wider than the corresponding topography feature. This is also an aspect of our Bardeen results ͓see Fig. 5͑d͔͒. Over the range of simulated tunneling currents, the tot signal changes from a flat protrusion of 5.7 Å FWHM with a 3.2 Å top plateau to a protrusion of 8.8 Å FWHM with a small dip of 0.05 Å on its top. This change of appearance with a maximum at the position of CO reflects the geometry of the 5 orbital: due to its short decay length, it is only visible at short distances, while at long distances, the 2 * orbital contributes the majority of the inelastic signal.

IV. SUMMARY
We have developed a method that permits us to include tip effects in IETS simulations at the Bardeen level of approximation using first-principles electronic structure calculations of the tip and sample wave functions. Its application to CO / Cu͑100͒ using a sharp pyramidal tip structure with W͑110͒ stacking confirms a quantitative improvement with measured data for the inelastic efficiency values, as compared to calculations done at the TH level. The implementation of the method is simple because it only requires us to evaluate the first order change in the sample wave functions due to the perturbing deformation potential, which can be expressed in terms of overlap integrals. Other tip structures with different terminations are foreseen to explore the possible selection rules for active vibrational modes in IETS.

APPENDIX
By assuming a linear dependence of the electron-vibration coupling in the nuclear coordinates ͑Q͒, 30 it is possible to approximate the wave functions for the displaced configurations as where a general scaling factor f ͑f ͓0,1͔͒ has been introduced for the normal mode ␦Q i that displaces the system out of the optimized geometry Q 0 . The deformation potential is obtained by effecting a displacement of the nuclear coordinates along the normal mode, of modulus the root mean square displacement. However, this value can be small for high-energy modes, and the deformation potential will have values within the error bar of the numerical calculation. To avoid this problem, it is more convenient to perform calculations for larger displacements, scaling them back to the root mean square displacement. Since this theory is based on the linear scaling of the deformation potential with the nuclear displacements, we can linearly scale the potential with the displacement by means of the constant factor f. The adopted labeling is the following: m, n indices run over bands, k over k points, and finally is used to specify the given spin state ͑ = 2 for spin polarized calculation, 1 oth-erwise͒. It is straightforward to show that upon summation of Eq. ͑A1͒, the unperturbed wave function can be recovered: FIG. 5. ͑Color online͒ Comparison between ͑a͒ TH and ͑b͒ Bardeen calculated total efficiencies maps for CO stretching on Cu͑100͒. ͑c͒ TH and ͑d͒ Bardeen efficiency profiles along ͓001͔ are also displayed for several contour values.
and that the change to the wave function itself can be explicitly evaluated by subtracting the wave functions for the displaced configurations appearing in Eq. ͑A1͒: At the right end of both Eqs. ͑A2͒ and ͑A3͒, we have introduced a shorter notation ͑that we are going to use in the following͒ where the atomic configurations, displaced by the ith normal mode, are referred to as "Ϯ" instead of "±f␦Q i ." Extra care must be exerted when dealing with wave functions coming from diagonalization of two different Hamiltonians as those appearing in Eqs. ͑A2͒ and ͑A3͒. In fact, in order to properly account for the different phases characterizing + and − , it is necessary to renormalize the wave function of one of the two displaced geometries with respect to the overlap between the expanded and the contracted geometries. 15 In doing so, Eq. ͑A2͒ is extended to read and analogously, Eq. ͑A3͒ yields Of course, it is also possible to renormalize + with respect to − and consistently working with the computed phase of − , obtaining in the end equivalent results, i.e., ͗ nk Taking the limits of Eq. ͑1͒ for ␦ → 0, the first order perturbation to the nth band at the kth k point for the spin of the sample can be written as where, in analogy to Ref. 15, the first and second term, on the right hand side of Eq. ͑A6͒ are referred to as principal ͑P͒ and secondary ͑I͒ contributions to the inelastic perturbation, respectively, Upon application of the Hellman-Feynman theorem and following the same approach as in Ref. 8, it is possible to define the coupling matrix element as which provides the possibility of evaluating the perturbation to the wave function ͑and eventually the change in conduc-tance͒ on the basis of overlap integrals between the displaced configurations. In fact, by substitution of Eq. ͑A8͒ into Eq. ͑A6͒ and using Eqs. ͑A4͒ and ͑A5͒, after some simple algebra, it is possible to explicitly rearrange the two terms appearing in Eq. ͑A7͒ to read where the ␦ functions have been approximated by a broadened Gaussian function of width for computational purposes, ⑀ nks = 1 2 ͑⑀ nks + + ⑀ nks − ͒, and the interband mixing coefficient C m,n is defined as follows: Once the two terms ͑principal P and secondary I͒ that constitute the inelastic perturbation are known, also the corresponding efficiency, i.e., the change in differential conductance relative to the unperturbed system, can be split into two terms on the basis of Eq. ͑A7͒: ine = P + I .

͑A12͒
The numerical handiness of Eq. ͑A7͒ emerges in that it allows us to evaluate the perturbation and eventually the efficiency due to exchange of two electrons mediated by the electron-vibration coupling ͑ ela ͒ on the basis of the secondary contribution to the inelastic perturbations, Eq. ͑A10͒: 15 ela = − 2 I .

͑A13͒
Being the total efficiency the sum of the inelastic and elastic terms, it straightforwardly follows that tot = ine + ela = P − I .

͑A14͒
When working in a Tersoff-Hamann framework, defined ͑r 0 , ⑀ F ͒ as the LDOS at the position r 0 for the unperturbed system:

͑A15͒
The total efficiency can therefore be calculated as

͑A19͒
where f͑⑀ n ͒ is the Fermi distribution function for the electronic state of energy ⑀ n and M ,n refers to the tunneling matrix element between tip ͑ ͒ and sample ͑ n ͒: In analogy with Eq. ͑A16͒, it is then possible to calculate the total efficiency by using Bardeen