# Orbital dependent electron tunneling within the atom superposition approach: Theory and application to W(110)

###### Abstract

We introduce an orbital dependent electron tunneling model and implement it within the atom superposition approach for simulating scanning tunneling microscopy (STM) and spectroscopy (STS). Applying our method, we analyze the convergence and the orbital contributions to the tunneling current and the corrugation of constant current STM images above the W(110) surface. In accordance with a previous study [Heinze et al., Phys. Rev. B 58, 16432 (1998)], we find atomic contrast reversal depending on the bias voltage. Additionally, we analyze this effect depending on the tip-sample distance using different tip models, and find two qualitatively different behaviors based on the tip orbital composition. As an explanation, we highlight the role of the real space shape of the orbitals involved in the tunneling. STM images calculated by our model agree well with Tersoff-Hamann and Bardeen results. The computational efficiency of our model is remarkable as the k-point samplings of the surface and tip Brillouin zones do not affect the computation time, in contrast to the Bardeen method.

###### pacs:

68.37.Ef, 71.15.-m, 73.63.-b## I Introduction

The experimental use of scanning tunneling microscopy (STM) and spectroscopy (STS) has recently gained a boost. The explanation of experimentally observed effects is not straightforward without a proper theoretical support hofer03rmp (); hofer03pssci (). One direction of recent theoretical research is focused on extracting surface local electronic properties from experimental STS data ukraintsev96 (); koslowski07 (); passoni09 (); ziegler09 (); koslowski09 (), which is the convolution of tip and sample electronic structures. Another research direction is concerned with the simulation of STM and STS by using different models mostly based on electronic structure data obtained from first principles.

Much work has been devoted to analyze the electron tunneling properties depending on the scanning tip. In Ref. hofer05sts () a theoretical method has been presented to separate the tip and sample contributions to STS. Ness and Gautier studied different metal tips and their interaction with metal surfaces in a tight-binding framework ness95jpcm1 (); ness95jpcm2 (); ness95prb (). Chen and Sacks investigated the effect of the tip orbitals on the corrugation of constant current STM images theoretically chen90 (); sacks00 (). While Chen pointed out that corrugation enhancement is expected for tip orbitals localized along the surface normal () direction ( and ), Sacks argued that tip states (, , , ) are responsible for this effect. Another work of Chen explained the corrugation inversion found on low Miller index metal surfaces by tip states chen92 (). Atomic contrast reversal has also been found above Xe atomic adsorbates mingo96 () and oxygen overlayers calleja04 () on metal surfaces. It was established that the character of the contrast depends on the tip-sample distance and on the tip geometry and electronic structure. The quality of the tip also plays a crucial role in the inelastic tunneling spectroscopy (IETS). Studying the CO molecule on Cu surfaces it has been found that the IETS intensity is close to the experiment including the full electronic structure of the tip teobaldi07 (), and the tip position and orbital symmetry can change the selection rules for the active vibrational modes in IETS garcia-lekue11 ().

The role of the electron orbitals has been considered in different electron transport models. Sirvent et al. developed a tight-binding model based on the Keldysh formalism for calculating the conductance in atomic point contacts and analyzed the effect of the orbitals sirvent96 (). The same methodology has been applied to STM junctions by Mingo et al. mingo96 (). Cerdá et al. presented an STM simulation method based on the Landauer-Büttiker formula buttiker85 () and the surface Green function matching technique cerda97 (). In these methods the decomposition of the conductance/current with respect to electron orbitals has been provided.

In the present work we consider a simple model for orbital dependent tunneling within the atom superposition approach for simulating STM palotas11stm () and STS palotas12sts (). The main idea of the paper is the introduction of a geometrical factor to account for a modified transmission due to electron orbital orientational overlap effects within a three-dimensional (3D) Wentzel-Kramers-Brillouin (WKB)-based theory. The reliability of this new method is demonstrated by the analysis of the tip-sample distance and bias voltage dependent corrugation reversal effect observed on the W(110) surface, where we find excellent agreement with a previous work heinze98 (). The computational efficiency of our method, which is justified in section III.5, enables us to study this effect in much more detail. We particularly focus on tip effects, and consider ideal tip models with different orbital symmetries and a more realistic W(110) tip. We find two qualitatively different corrugation inversion behaviors based on the tip orbital composition. Our results indicate that anticorrugation on the W(110) surface can not only be observed at negative bias voltages but also at positive bias at reasonably short tip-sample distances.

The paper is organized as follows: The theoretical model of the orbital dependent tunneling within the atom superposition approach is presented in section II. Applying this method we investigate the convergence and the orbital contributions of the tunneling current as well as the corrugation reversal of the W(110) surface depending on the applied bias voltage and the tip-sample distance in section III. We also report a comparison of STM images calculated by our model to Tersoff-Hamann tersoff83 (); tersoff85 () and Bardeen bardeen61 () results in section III.5. Summary of our findings is found in section IV.

## Ii Orbital dependent tunneling model within 3D WKB theory

Recently, Palotás et al. developed a 3D atom superposition approach for simulating spin-polarized STM (SP-STM) palotas11stm () and spin-polarized STS (SP-STS) palotas12sts () based on previous theories passoni09 (); wortmann01 (); yang02 (); smith04 (); heinze06 (). In the model it is assumed that electrons tunnel through one tip apex atom, and contributions from individual transitions between this apex atom and each of the surface atoms are summed up assuming the one-dimensional (1D) WKB approximation for electron tunneling processes in all these transitions. The key input is the projected electron density of states (PDOS) of the tip apex atom and of the sample surface atoms obtained from ab initio electronic structure calculations. In the present paper first we review the tunneling current and the differential conductance expressions based on the independent orbital approximation for the vacuum decay of electron states, and then we extend this model to include a simple orbital dependent tunneling transmission. We consider the non-spin-polarized part of the tunneling only, however, this theory can be applied to SP-STM and SP-STS in the future.

Assuming elastic tunneling at K temperature, the tunneling current at the tip position and at bias voltage is given by

(1) |

The integrand is the so-called virtual differential conductance,

Here is the elementary charge, is the Planck constant, and and are the Fermi energies of the tip and the sample surface, respectively. ensures that the is correctly measured in the units of . has been chosen to 1 eV but its actual value has to be determined comparing simulation results to experiments or to calculations using more sophisticated tunneling models. Note that the exact choice of this scaling factor has absolutely no qualitative influence on the reported results, and the comparison of current values to Bardeen results confirms our choice, see section III.5. The sum over corresponds to the atomic superposition and has to be carried out, in principle, over all surface atoms. However, convergence tests showed that including a relatively small number of atoms in the sum provides converged values palotas11sts (). is the energy and bias dependent tunneling transmission function, which also depends on the distance between the tip apex and the surface atom labeled by with position vector . The tip and sample electronic structures are included into this model via projected DOS (PDOS) onto the atoms, i.e., and denote projected charge DOS onto the tip apex and the th surface atom, respectively. They can be obtained from any suitable electronic structure method.

Taking the derivative of Eq.(1) with respect to the bias voltage the differential conductance is obtained. It can be written at the tip position and at bias voltage as the sum of three terms,

(3) |

Here has been defined in Eq.(II), and and are the background and tip-derivative terms,

respectively palotas12sts (). The background term, which contains the bias-derivative of the transmission function, is usually taken into account in recent STS theories passoni09 (); passoni07 (); donati11 (), while the tip-derivative term containing the energy derivative of the tip DOS is rarely considered in the recent literature.

In the spirit of the independent orbital approximation the transmission probability for electrons tunneling between states of atom on the surface and the tip apex is of the simple form,

(6) |

This corresponds to a spherical exponential decay of the electron wave functions irrespective of their orbital symmetry tersoff83 (); tersoff85 (); heinze06 (). Assuming the same effective rectangular potential barrier between the tip apex and each surface atom, the vacuum decay can be written as

(7) |

where the electron’s mass is , is the reduced Planck constant, and and are the average electron work function of the sample surface and the local electron work function of the tip apex, respectively. The method of determining the electron work functions from the calculated local electrostatic potential is reported in Ref. palotas11stm ().

Next, we extend this tunneling model by taking advantage of the orbital decomposition of the electronic structure data and the real space shape of the electron orbitals. The PDOS of the sample surface atoms and the tip apex can be decomposed according to orbital symmetry, i.e., real spherical harmonics, , as

(8) | |||

(9) |

Similar decomposition of the Green functions has been employed in the linear combination of atomic orbitals (LCAO) scheme by Refs. mingo96 (); sirvent96 (). Assuming such an orbital decomposition, the virtual differential conductance can be generalized as

where, additionally to the atomic superposition (sum over ) we sum up each transition with an orbital dependent tunneling transmission: gives the electron tunneling probability at energy from the tip apex orbital to the orbital of the th surface atom at positive bias voltage (), and from the orbital of the th surface atom to the tip apex orbital at negative bias (). can be defined in different ways based on physical arguments. We consider the following form,

(11) |

for each surface atom tip apex 1D electron transition. Here, the energy and bias dependent part corresponds to the spherical exponential decay considered in Eq.(6), and is independent of the orbital symmetry. This is multiplied by an orbital dependent expression , which depends on the spatial arrangement of the sample atoms relative to the tip apex and all the orbital shapes involved in the tunneling. The angular dependence on and comes into play in the following way: Let us consider one transition between surface atom and the tip apex along the line connecting the two atoms. A particular geometry is shown in Figure 1. For brevity, we omitted the notation of the surface atom. For both atoms a local coordinate system can be set up, and the angular dependence of the atomic orbital wave functions are defined in the corresponding coordinate system, as summarized in Table 1. defines a vector pointing from the surface atom toward the tip apex, and it can be represented by the coordinates in the spherical coordinate system centered on the th surface atom. Taking the opposite connecting vector from the tip apex toward the surface atom, its coordinates are in the spherical coordinate system centered on the tip apex. According to Figure 1, , , and can be calculated as

(12) | |||

(13) | |||

(14) |

using the global coordinates and . Considering above, accounts for the modification of the perfect spherical exponential decay along the connecting line through the angular dependence of the atomic orbitals as

(15) |

where are the real spherical harmonics summarized in Table 1. They were chosen in such a way that . This factor takes the effect of the directional tunneling between real space orbitals into account. The physical motivation is the angular dependence of the electron density of the orbitals in the given tunneling direction, which modifies the tunneling transmission. The maximum is obtained if the angular distributions of the electron density according to the given orbital symmetries () around both the sample surface atom and the tip apex have maxima along the line of the tunneling direction. This is always the case for - type of tunneling irrespective of the relative position of the tip apex and sample surface atoms, i.e., we observe perfect spherical exponential decay between tip and sample orbitals. In some particular geometries can be obtained even for other type of orbitals, e.g., if the tip apex is precisely above surface atom , i.e., , then for all of the following combinations: . On the other hand, if the tip apex is above surface atom then orbitals with nodal planes orthogonal to the surface have zero contribution to the tunneling from this particular surface atom, i.e., a reduced effective tunneling transmission is obtained sacks00 (). Note that the independent orbital approximation corresponds to for all transitions, i.e., the same tunneling transmission is assumed between all orbitals. Within our orbital dependent tunneling approach ideal tip models with particular orbital symmetries can be considered, i.e., orbital symmetry corresponds to the choice of and . More realistic tips can be obtained by explicitly calculating the orbital decomposition of the tip apex PDOS in model tip geometries, e.g., in the present work a model W(110) tip is used.

Our theory is thus an extension of the atom superposition STM/STS approach considering tunneling between directional orbitals. The tunneling current and the differential conductance can be calculated at the tip position and at bias voltage as the sum of all transitions between real space orbitals,

(16) | |||

(17) |

respectively, with

This decomposition enables the analysis of the orbital contributions to the total tunneling current and to the differential conductance. In relation to Chen’s derivative rule chen90 () we can state that while it is formulated inspired by the Tersoff-Hamann model, and calculates the tunneling transmission as the absolute value square of the tunneling matrix element that is proportional to the sample wave function derivative with respect to the real space coordinate corresponding to the given tip orbital symmetry (), our transmission function also depends on the sample surface atoms’ orbital symmetry . Moreover, the electronic structure of the tip apex is included in our theory via the PDOS.

Our method does not account for multiple scattering palotas05 () and interference effects mingo96 (); cerda97 (), which could be important for certain systems. Therefore, it is expected that our model works well on simple metals, possibly on molecular adsorbates on surfaces as well but not on materials with strong band structure or Fermi surface effects. Other limitation is the uniform vacuum decay of the electron states for different orbital symmetries in Eq.(11). This could be improved in the future. Still, the model in its present form provides comparable results to more sophisticated tunneling models (Tersoff-Hamann, Bardeen) as will be presented in section III.5.

Note that the presented method can also be applied to magnetic systems taking into account the orbital-decomposed magnetization PDOS of the tip and sample palotas11stm (); palotas12sts () together with the orbital dependent tunneling transmission in Eq.(11). As it was pointed out by Ferriani et al. ferriani10tip () the spin polarization in the vacuum can have an opposite sign than within the tip apex atom, and this sign change is also accompanied by different dominating orbital characters. Thus, the consideration of an orbital dependent tunneling transmission might be a better model for describing electron transport through a magnetic tunnel junction. We return to the related spin-polarized STM/STS model in the future.

## Iii Results and discussion

In order to demonstrate the reliability of our orbital dependent tunneling model we consider a W(110) surface. This surface has technological importance as it is widely used as substrate for thin film growth, see e.g., Refs. heinze98 (); bode07 (). As it was pointed out by Heinze et al. heinze98 () the determination of the position of surface atomic sites is not straightforward as atomic resolution is lost at negative bias voltages, and a bias-dependent corrugation reversal has been predicted. This means that normal and anticorrugated constant current STM images can be obtained in certain bias voltage ranges, and the W atoms do not always appear as protrusions in the images. It was shown that a competition between states from different parts of the surface Brillouin zone is responsible for this effect heinze98 (); heinze99 (). We reinvestigate this corrugation reversal effect as it provides a challenge for our orbital dependent tunneling model. We find excellent agreement with the results of Ref. heinze98 (), where an -wave tip has been used, and we study this effect in more detail. We particularly focus on tip effects, and consider ideal tip models with different orbital symmetries, and a more realistic W(110) tip. We find two qualitatively different corrugation inversion behaviors based on the tip orbital composition. Explaining our findings we highlight the role of real space orbital orientational overlaps between the surface and the tip rather than considering electron states in the reciprocal space, thus a different kind of understanding is provided. Finally, by comparing STM images to results of more sophisticated tunneling models we find good agreement.

### iii.1 Computational details

We performed geometry relaxation and electronic structure calculations based on the density functional theory (DFT) within the generalized gradient approximation (GGA) implemented in the Vienna Ab-initio Simulation Package (VASP) VASP2 (); VASP3 (); hafner08 (). A plane wave basis set for electronic wave function expansion together with the projector augmented wave (PAW) method kresse99 () has been applied, and the exchange-correlation functional is parametrized according to Perdew and Wang (PW91) pw91 (). The electronic structures of the sample surface and the tip have been calculated separately.

We model the W(110) surface by a slab of nine layers, where the two topmost W layers have been fully relaxed. After relaxation the W-W interlayer distance between the two topmost layers is reduced by 3.3%, while the underneath W-W interlayer distance increased by 1.1% compared to bulk W. A separating vacuum region of 18 Å width in the surface normal () direction has been set up between neighboring supercell slabs. The average electron work function above the surface is calculated to be eV. We used a Monkhorst-Pack (MP) monkhorst () k-point grid for obtaining the orbital-decomposed projected electron DOS onto the surface W atom, . The same k-set has been used for calculating the sample electron wave functions for the Tersoff-Hamann and Bardeen simulations. The unit cell of the W(110) surface (shaded area) and the rectangular scan area for the tunneling current simulation are shown in Figure 2. In our calculations we used the experimental lattice constant pm. Moreover, the surface top (T) and hollow (H) positions are explicitly shown.

We considered different tip models. The orbital-independent ideal tip is characterized by and , so that . This ideal electronically flat tip represents the limiting case of the independent orbital approximation used in previous atom superposition tunneling models palotas11stm (); palotas12sts (); heinze06 (); palotas11sts (). In order to study the effect of the orbital dependent tunneling other tip models are needed. First, we consider ideal tip models having a particular orbital symmetry . In this case is calculated following Eq.(15), and for the energy dependence of the tip PDOS, and are assumed. More realistic tips can also be employed by calculating the orbital decomposition of the tip apex PDOS in model tip geometries, and using Eq.(15) for the orbital dependent transmission factor. We used a blunt W(110) tip. Motivated by a previous simulation teobaldi12 (), it has been modeled by a slab consisting of three atomic layers having one W apex atom on both surfaces, i.e., with a double vacuum boundary. In this system the apex atoms have been relaxed on both sides. The adatom-topmost layer vertical distance decreased by 19.3% compared to bulk W. The interaction between apex atoms in neighboring supercells in the lateral direction is minimized by choosing a surface cell, and a 17.9 Å wide separating vacuum region in the direction. The local electron work function above the tip apex was assumed to be eV. Moreover, an MP k-point grid has been chosen for calculating the orbital-decomposed projected DOS onto the apex atom, . The same k-point sampling has been used for obtaining the tip electron wave functions for the Bardeen calculation.

STM images were simulated employing our model, and the Tersoff-Hamann tersoff83 (); tersoff85 () and Bardeen bardeen61 () methods implemented in the BSKAN code hofer03pssci (); palotas05 (). Scattering up to first order palotas05 () did not affect the quality of the images. Using our model the tunneling current has been calculated in a box above the rectangular scan area shown in Figure 2 containing 99000 () grid points with a Å lateral and Å vertical resolution. The electron local density of states (LDOS) was calculated above the same scan area in a box of grid points using the Tersoff-Hamann method with the same spatial resolution as above. For the calculation of the tunneling current employing the Bardeen method a box of grid points above the half of the rectangular scan area has been chosen in order to speed up the simulation. In this case the lateral resolution remains Å, and the vertical resolution is Å. The constant current contours are extracted following the method described in Ref. palotas11stm (). All of the STM images will be presented above the full rectangular scan area.

### iii.2 Convergence properties

Previously, the convergence of the part of the differential conductance has been investigated with respect to the number of surface atoms involved in the summation of the orbital-independent atomic superposition formula palotas11sts (). Due to the spherical exponential decay assumed for the electron wave functions a rapid convergence was found. We report a similar convergence test for the orbital-dependent tunneling approach comparing different tip models. In order to take into account a wide energy range around the Fermi level we calculated the tunneling current at -2.5 V and +2.5 V bias voltages at Å above a surface W atom, and averaged these current values. We considered ideal tips of the orbital-independent model, and with , , and symmetry, as well as the W(110) tip. In order to obtain comparable results we normalized the averaged current for each tip calculation. The convergences of the normalized averaged current with respect to the lateral distance on the surface, , characteristic for the number of atoms involved in the atomic superposition, are shown in Figure 3. By calculating the current, contributions from surface atoms within a radius of measured from the W atom below the tip apex are summed up (sum over ). We find that the orbital-independent, the -type, and the W(110) tips behave quite similarly concerning the current convergence, while for the - and -type tips a faster convergence is found. This rapid convergence can be explained by the more localized character of the latter tip orbitals in the direction normal to the sample surface (). On the other hand, the orbital-independent tip with is a good approximation for the -type tip (with index ), where the spherically decaying transmission function part is still dominant, i.e., because . In case of the W(110) tip, electronic states of all considered symmetries have a contribution, and their relative importance is not only determined by the transmission function via the orbital shapes but also by the product of the symmetry-decomposed electron PDOS of the surface and the tip. In general, the orbitals localized in different than the direction can show a slower current convergence than the orbitals. However, the partial PDOS of such states is relatively low, and interestingly we obtain a similar current convergence in the studied energy range as for the -type tip. Choosing different bias voltages for the W(110) tip, thus different electron states involved in the tunneling, we found current convergences dissimilar to the -type tip behavior. The convergence can be slower or faster than obtained for the -type tip depending on the partial PDOS of each directional orbital in the given energy range.

Based on the convergence tests, atom contributions within at least Å distance from the surface-projected tip position shall be considered. In case of calculating STM images, Å has to be measured from the edge of the scan area in all directions in order to avoid distortion of the image, thus involving 67 surface atoms in the atomic superposition. For brevity, in the following we use the same surface atoms to calculate single point tunneling properties as well.

### iii.3 Orbital Contributions

Let us analyze the relative importance of all transitions in determining the total tunneling current at different tip positions. From this analysis we obtain a qualitative picture about the role of the different atomic orbitals in the construction of the tunneling current. The current contributions can be calculated according to Eq.(II). These can be represented by a current histogram that gives the percentual contributions of the individual transitions to the total current. Figure 4 shows such histograms using the W(110) tip at = -0.1 V bias voltage Å above two different tip positions: part a) corresponds to the tip apex above the surface top position, and part b) to the tip apex above the surface hollow position, T and H in Figure 2, respectively. We obtain a matrix from the considered orbitals, which are denoted by numbers 1 to 9 following the indices reported in Table 1. We find that most contributions are due to the (1), (3), (6), (7), and (8) orbitals and their combinations. The largest contribution to the current is given by the (7-7) transition, 31 and 20 per cent above the top and hollow positions, respectively. Concomitantly, above the hollow position, the relative importance of both tip and sample (6) and (8) orbitals is increased as it is expected from the geometrical setup, i.e., the (6-6), (6-7 and 7-6), (8-8), and (7-8 and 8-7) contributions correspond to larger orientational overlap of the mentioned tip and sample orbitals if the tip is above the hollow position rather than above the top position as suggested by the geometry in Figure 2 and Eq.(15). Thus, our simple orbital dependent tunneling model captures the effect of the localized orbitals and goes beyond the spherical Tersoff-Hamann model. Note that if a larger bias voltage is considered, i.e., the electronic states are somewhat averaged, then the independent orbital approach might turn out to be a good approximation heinze06 ().

### iii.4 Atomic contrast reversal

The role of the localized orbitals can best be demonstrated by reinvestigating the corrugation inversion phenomenon found, e.g., on (100) mingo96 (), (110) heinze98 (), and (111) ondracek12 () metal surfaces. Chen explained this effect as a consequence of tip states chen92 (). According to Heinze et al. heinze98 () under certain circumstances the apparent height of W atoms at the surface top position () can be larger or smaller than the apparent height of the surface hollow position () at constant current () condition. (For the surface top (T) and hollow (H) positions, see Figure 2.) Thus, the sign change of is indicative for the corrugation inversion. Obviously, corresponds to a normal STM image, where the W atoms appear as protrusions, and to an anticorrugated image. Since the tunneling current is monotonically decreasing with the increasing tip-sample distance, we can obtain information about the occurrence of the corrugation inversion simply by calculating the current difference between tip positions above the top and hollow sites of the W(110) surface. The current difference at tip-sample distance and at bias voltage is defined as

(20) |

This quantity can be calculated for specific tips, and we call the contour as the corrugation inversion map. This gives the combinations where the corrugation inversion occurs. The sign of provides the corrugation character of an STM image in the given regime. Due to the monotonically decreasing character of the tunneling current, corresponds to , i.e., normal corrugation, and similarly corresponds to and anticorrugation.

First, we calculated using the independent orbital approximation and Eq.(6) for the tunneling transmission, and found that is always positive. This means that the spherical exponential decay itself can not account for the observed corrugation inversion effect, and the W atoms always appear as protrusions in STM images calculated with this model. However, considering the orbital dependent tunneling transmission in Eq.(11) we find evidence for the corrugation inversion effect, thus highlighting the role of the real space shape of electron orbitals involved in the tunneling. Figure 5 shows contours calculated with different tip models in the Å Å tip-sample distance and [-2 V,+2 V] bias voltage range. Before turning to the analysis of the results obtained with previously not considered tip models let us compare our results with those of Heinze et al. heinze98 (), where an -wave tip model has been used. They found corrugation reversal at around -0.4 V at Å tip-sample distance, and above that voltage normal while below anticorrugated STM images were obtained. Our model with an -tip provides the same type of corrugation reversal at -0.21 V at the same distance as can be seen in part a) of Figure 5 (curve with filled square symbol). These bias values are in reasonable agreement particularly concerning their negative sign. At this range atomic resolution is difficult to achieve experimentally, which is an indication for being close to the corrugation inversion regime heinze98 (). On the other hand a linear dependence of the corrugation reversal voltage and the tip-sample distance has been reported by Heinze et al. Å V to Å V. Our model qualitatively reproduces this linear dependence in the same bias range though the quantitative values are somewhat different.

Calculating the corrugation inversion maps with more tip models, we find two distinct behaviors depending on the tip orbital composition. Parts a) and b) of Figure 5 show these. While the tip models in part a) can show corrugation inversion in the whole studied bias range, this effect does not occur at positive bias voltages for tips in part b). Moreover, anticorrugation () is observed in the large tip-sample distance region ( Å) in both parts. This is in accordance with the prediction of Ref. heinze99 () based on the analysis of the competing electron states in the surface Brillouin zone of an Fe(001) surface. In the Å range, however, the corrugation character in the two parts of Figure 5 is remarkably different. In part a), normal corrugation is found close to the surface, which reverts only once with increasing tip-sample distance for the tip models with a single orbital symmetry in the full studied bias range. The W(110) tip behaves similarly below +1.7 V, while above there is a double reversal of the corrugation character as the tip-sample distance increases. This indicates that anticorrugation can be expected at short tip-sample distances (3.5 Å-5 Å) at around +2 V. On the other hand, the tip models in part b) always show anticorrugation at positive bias voltages, and below -0.05 V they provide corrugation characters starting from anticorrugation, then normal corrugation, and again anticorrugation with increasing tip-sample distance. These different behaviors can be attributed to the tip orbital characters. It is interesting to notice that none of the considered tip orbitals in part b) are localized in the -direction, and they have nodal planes either in the plane ( and ) or in the plane ( and ) or in the and planes (). On the other hand, in part a) there are tips which are localized in the -direction ( and ) or having nodal planes in both the and planes () as well as the spherical tip and the W(110) tip that contains all type of orbitals with energy dependent partial PDOS functions. The particular tip nodal planes restrict the collection of surface atom contributions to specific regions on the sample surface. Furthermore, by changing the tip-sample distance, the orientational overlaps between the tip and sample orbitals change, and according to our model some localized orbitals gain more importance in the tunneling contribution, see also Figure 4. Since we calculate the current difference between tip positions above the surface top and hollow sites, the complex tip-sample and bias voltage dependent effect of the real space orbitals on the tunneling can be visualized via the corrugation inversion maps.

Concerning tips with and orbital symmetry, Heinze et al. heinze98 () calculated a corrugation enhancement factor of 2 and 6.25, respectively, based on Chen’s derivative rule chen90 (). Moreover, they argued that the corrugation inversion map should be practically identical to the one obtained by using the -tip model, and the corrugation values just have to be scaled up by these factors. On the contrary based on our orbital dependent tunneling model we find that the and tips provide qualitatively different corrugation inversion maps, i.e., although their bias dependent shape is similar to the one of the -tip, their tip-sample distance is systematically pushed to larger values, see part a) of Figure 5. This is due to the more localized character of these tip orbitals in the -direction.

Corrugation inversion with the tip occurs at the largest tip-sample distance. A possible explanation can be based on its and nodal planes. While above the top position only the underlying W atom, above the hollow position all four nearest neighbor W atoms give zero contribution to the current, thus is expected to be higher than at small tip-sample distances. To overcome this effect the tip has to be moved farther from the surface since then the relative importance of the nearest neighbor contributions decays rapidly compared to other parts of the surface.

Apart from above findings we obtain corrugation inversion also in the positive bias range at enlarged tip-sample distances for the , , , and W(110) tips considered in part a) of Figure 5. This is most probably due to the surface electronic structure. Note that this effect is even more difficult to capture in experiments as the corrugation values themselves decay rapidly with increasing tip-sample distance.

### iii.5 STM images - Comparison to other tunneling models

In order to demonstrate the corrugation inversion more apparently, constant current STM images can be simulated. As it is clear from Figure 5, any type of crossing of the contour results in the occurrence of the corrugation reversal. In experiments two ways can be considered to record STM images in the normal and anticorrugated regimes: 1) keep the tip-sample distance constant, and change the bias voltage ; or 2) keep the bias voltage constant, and change the tip-sample distance. Respectively, these modes correspond to a horizontal and a vertical crossing of the contour in the plane in Figure 5. Heinze et al. followed the first method in their simulations heinze98 (). However, as the second option seems to be experimentally more feasible and needs less calculations as well, we simulated STM images at a fixed bias voltage of -0.25 V.

In Figure 6 STM images are compared using our model assuming an -type tip [first row a)-c)] to those calculated by the Tersoff-Hamann method [second row d)-f)]. We find that the images are in good qualitative agreement for the a)-d), b)-e), and c)-f) pairs, respectively. In parts a) and d), at a tip-sample distance of about 3.80 Å, the apparent height of the W atom is larger than the one of the hollow position, i.e. . This resembles normal corrugation. Moving the tip farther from the surface, we obtain the corrugation inversion and striped images. These are shown in parts b) and e) of Figure 6. We find that our method reproduces the corrugation inversion effect at almost the same tip-sample distance (4.15 Å) as the Tersoff-Hamann model predicts (4.21 Å). Increasing the tip-sample distance further we enter the anticorrugation regime, and the apparent height of the W atom is smaller than the one of the hollow position, i.e., . Such images are shown in parts c) and f). Note that all of the simulated STM images in Figure 6 are in good qualitative agreement with Ref. heinze98 (). The corrugation of the individual current contours has also been calculated: a) pm, b) pm, c) pm, d) pm, e) pm, f) pm. We find that our model gives approximately one order of magnitude less corrugation than the Tersoff-Hamann method. Note, however, that the small corrugation amplitudes using our method are in good agreement with Ref. heinze98 (), where they report pm close to the contrast reversal.

As we have seen, the corrugation inversion effect already occurs considering the electronic structure of the sample only. However, Figure 5 indicates that different tips can modify its tip-sample distance and bias voltage dependence quite dramatically. In Figure 7 STM images are compared using our model [first row a)-c)] to those calculated by the Bardeen method [second row d)-f)] explicitly taking the electronic structure of the W(110) tip in both cases into account. We find that the images are in good qualitative agreement for the b)-e) and c)-f) pairs. In parts a) and d), at a tip-sample distance of about 4.50 Å, the agreement is weaker, however, the normal corrugation is more pronounced in our model: The corrugation amplitude of part a) pm is much larger than that of part d) pm. Moreover, as the current values of 6.3 nA (our model) and 4.4 nA (Bardeen) are comparable to each other at the given tip-sample separation, the choice of eV in Eq.(II) is confirmed. Note that employing our model, a better qualitative agreement to the image of part d) has been found at a larger tip-sample separation, i.e., closer to the corrugation inversion. This inversion is demonstrated in parts b) and e) of Figure 7. Again, we obtain striped images. Note, however, that the stripes with larger apparent height correspond to the atomic rows, in contrast to what has been found in parts b) and e) of Figure 6, where the atomic and hollow sites appeared as depressions. This difference is definitely due to the effect of the W tip, which was not considered in Figure 6. On the other hand, we find good agreement concerning the tip-sample distance of the corrugation inversion: 5.80 Å in our model, and 5.55 Å calculated by the Bardeen method. Parts c) and f) of Figure 7 correspond to anticorrugated images. In this tip-sample distance regime the extremely small corrugation amplitudes are in good agreement between our model and the Bardeen method: pm in parts b), c), f), and pm in part e).

Finally, we compared computation times between our model and the Bardeen method, and found the following:

1) our orbital dependent model, grid points, 1 CPU, time=229 s;

2) Bardeen method in BSKAN code, grid points, 4 CPUs, time=173866 s.

Normalizing to the same real space grid points we obtain that our method is 2425 times faster using 1 CPU than using 4 CPUs
for the Bardeen calculation. As the 4 CPUs calculations are roughly 3.5 times faster than the 1 CPU ones in our computer cluster,
a remarkable 1 CPU equivalent time boost of about 8500 is obtained for our method compared to the Bardeen for the given
surface-tip combination. While the k-point samplings of the surface and tip Brillouin zones affect the computation time of the
Bardeen method due to the enhanced number of transitions as the number of k-points increases, the computation time
of our model is insensitive to the number of k-points as the PDOS of the tip apex and surface atoms are used. The energy dependent
PDOS functions have the same data structure, no matter of the number of the constituting electron states obtained by the
k-summation palotas11stm (). This is a great computational advantage of our model. Of course, the quality of the results
depends on the k-point samplings. Moreover, please note the further potential that our method can be parallelized in the future.

Thus, employing our new computationally efficient orbital dependent tunneling model we could reproduce and reinvestigate the corrugation inversion effect observed on the W(110) surface. Although this effect is driven by the surface electronic structure, we showed that different tips can drastically modify its tip-sample distance and bias voltage dependence.

## Iv Conclusions

We developed an orbital dependent electron tunneling model and implemented it within the atom superposition approach based on 3D WKB theory, for simulating STM and STS. Applying our method we analyzed the convergence and the orbital contributions to the tunneling current above the W(110) surface. We found that the contribution is the largest, and depending on the tip position other states can gain importance as well. We also studied the corrugation inversion effect. Using the independent orbital approximation no corrugation reversal has been obtained at all. Employing the orbital dependent model we found corrugation reversals depending on the bias voltage in accordance with the work of Heinze et al. heinze98 (), and also on the tip-sample distance. Explaining this effect we highlighted the role of the real space shape of the orbitals involved in the tunneling. Moreover, we calculated corrugation inversion maps considering different tip models, and found two qualitatively different behaviors based on the tip orbital composition. Our results indicate that using a W tip anticorrugation can not only be observed at negative bias voltages but also at positive bias at reasonably short tip-sample distances. Simulation of STM images made the corrugation inversion effect more apparent. A good agreement has been found by comparing STM images calculated by our model to Tersoff-Hamann and Bardeen results. The computational efficiency of our model is remarkable as the k-point samplings of the surface and tip Brillouin zones do not affect the computation time, in contrast to the Bardeen method. Extending this orbital dependent tunneling model to magnetic junctions is expected to enable the study of the interplay of real space orbital and spin polarization effects in SP-STM and SP-STS experiments in the future.

## V Acknowledgments

The authors thank W. A. Hofer, G. Teobaldi and C. Panosetti for useful discussions. Financial support of the Magyary Foundation, EEA and Norway Grants, the Hungarian Scientific Research Fund (OTKA PD83353, K77771), the Bolyai Research Grant of the Hungarian Academy of Sciences, and the New Széchenyi Plan of Hungary (Project ID: TÁMOP-4.2.2.B-10/1–2010-0009) is gratefully acknowledged. Furthermore, partial usage of the computing facilities of the Wigner Research Centre for Physics, and the BME HPC Cluster is kindly acknowledged.

## References

- (1)
- (2) W. A. Hofer, A. S. Foster, and A. L. Shluger, Rev. Mod. Phys. 75, 1287 (2003).
- (3) W. A. Hofer, Prog. Surf. Sci. 71, 147 (2003).
- (4) V. A. Ukraintsev, Phys. Rev. B 53, 11176 (1996).
- (5) B. Koslowski, C. Dietrich, A. Tschetschetkin, and P. Ziemann, Phys. Rev. B 75, 035421 (2007).
- (6) M. Passoni, F. Donati, A. Li Bassi, C. S. Casari, and C. E. Bottani, Phys. Rev. B 79, 045404 (2009).
- (7) M. Ziegler, N. Néel, A. Sperl, J. Kröger, and R. Berndt, Phys. Rev. B 80, 125402 (2009).
- (8) B. Koslowski, H. Pfeifer, and P. Ziemann, Phys. Rev. B 80, 165419 (2009).
- (9) W. A. Hofer and A. Garcia-Lekue, Phys. Rev. B 71, 085401 (2005).
- (10) H. Ness and F. Gautier J. Phys. Condens. Matter 7, 6625 (1995).
- (11) H. Ness and F. Gautier J. Phys. Condens. Matter 7, 6641 (1995).
- (12) H. Ness and F. Gautier Phys. Rev. B 52, 7352 (1995).
- (13) C. J. Chen, Phys. Rev. B 42, 8841 (1990).
- (14) W. Sacks, Phys. Rev. B 61, 7656 (2000).
- (15) C. J. Chen, Phys. Rev. Lett. 69, 1656 (1992).
- (16) N. Mingo, L. Jurczyszyn, F. J. Garcia-Vidal, R. Saiz-Pardo, P. L. de Andres, F. Flores, S. Y. Wu, and W. More, Phys. Rev. B 54, 2225 (1996).
- (17) F. Calleja, A. Arnau, J. J. Hinarejos, A. L. Vázquez de Parga, W. A. Hofer, P. M. Echenique, and R. Miranda, Phys. Rev. Lett. 92, 206101 (2004).
- (18) G. Teobaldi, M. Peñalba, A. Arnau, N. Lorente, and W. A. Hofer, Phys. Rev. B 76, 235407 (2007).
- (19) A. Garcia-Lekue, D. Sanchez-Portal, A. Arnau, and T. Frederiksen, Phys. Rev. B 83, 155417 (2011).
- (20) C. Sirvent, J. G. Rodrigo, S. Vieira, L. Jurczyszyn, N. Mingo, and F. Flores, Phys. Rev. B 53, 16086 (1996).
- (21) M. Büttiker, Y. Imry, R. Landauer, and S. Pinhas, Phys. Rev. B 31, 6207 (1985).
- (22) J. Cerdá, M. A. Van Hove, P. Sautet, and M. Salmeron, Phys. Rev. B 56, 15885 (1997).
- (23) K. Palotás, W. A. Hofer, and L. Szunyogh, Phys. Rev. B 84, 174428 (2011).
- (24) K. Palotás, W. A. Hofer, and L. Szunyogh, Phys. Rev. B 85, 205427 (2012).
- (25) S. Heinze, S. Blügel, R. Pascal, M. Bode, and R. Wiesendanger, Phys. Rev. B 58, 16432 (1998).
- (26) J. Tersoff and D. R. Hamann, Phys. Rev. Lett. 50, 1998 (1983).
- (27) J. Tersoff and D. R. Hamann, Phys. Rev. B 31, 805 (1985).
- (28) J. Bardeen, Phys. Rev. Lett. 6, 57 (1961).
- (29) D. Wortmann, S. Heinze, P. Kurz, G. Bihlmayer, and S. Blügel, Phys. Rev. Lett. 86, 4132 (2001).
- (30) H. Yang, A. R. Smith, M. Prikhodko, and W. R. L. Lambrecht, Phys. Rev. Lett. 89, 226101 (2002).
- (31) A. R. Smith, R. Yang, H. Yang, W. R. L. Lambrecht, A. Dick, and J. Neugebauer, Surf. Sci. 561, 154 (2004).
- (32) S. Heinze, Appl. Phys. A 85, 407 (2006).
- (33) K. Palotás, W. A. Hofer, and L. Szunyogh, Phys. Rev. B 83, 214410 (2011).
- (34) M. Passoni and C. E. Bottani, Phys. Rev. B 76, 115404 (2007).
- (35) F. Donati, S. Piccoli, C. E. Bottani, and M. Passoni, New J. Phys. 13, 053058 (2011).
- (36) K. Palotás and W. A. Hofer, J. Phys. Condens. Matter 17, 2705 (2005).
- (37) P. Ferriani, C. Lazo, and S. Heinze, Phys. Rev. B 82, 054411 (2010).
- (38) M. Bode, M. Heide, K. von Bergmann, P. Ferriani, S. Heinze, G. Bihlmayer, A. Kubetzka, O. Pietzsch, S. Blügel, and R. Wiesendanger, Nature 447, 190 (2007).
- (39) S. Heinze, X. Nie, S. Blügel, and M. Weinert, Chem. Phys. Lett. 315, 167 (1999).
- (40) G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996a).
- (41) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b).
- (42) J. Hafner, J. Comput. Chem. 29, 2044 (2008).
- (43) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
- (44) J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
- (45) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
- (46) G. Teobaldi, E. Inami, J. Kanasaki, K. Tanimura, and A. L. Shluger, Phys. Rev. B 85, 085433 (2012).
- (47) M. Ondráček, C. González, and P. Jelínek, J. Phys. Condens. Matter 24, 084003 (2012).

orbital | index | definition | |
---|---|---|---|

1 | |||

2 | |||

3 | |||

4 | |||

5 | |||

6 | |||

7 | |||

8 | |||

9 |