First-principles modelling of scanning tunneling microscopy using non-equilibrium Green’s functions

Abstract

The investigation of electron transport processes in nano-scale architectures plays a crucial role in the development of surface chemistry and nano-technology. Experimentally, an important driving force within this research area has been the concurrent refinements of scanning tunneling microscopy (STM) techniques. The theoretical treatment of the STM operation has traditionally been based on the Bardeen and Tersoff–Hamann methods which take as input the single-particle wave functions and eigenvalues obtained from finite cluster or slabs models of the surface-tip interface. Here, we present a novel STM simulation scheme based on non-equilibrium Green’s functions (NEGF) and Wannier functions which is both accurate and very efficient. The main novelty of the scheme compared to the Bardeen and Tersoff–Hamann approaches is that the coupling to the infinite (macroscopic) electrodes is taken into account. As an illustrating example we apply the NEGF-STM method to the Si(001)-(2×1):H surface with sub-surface P doping and discuss the results in comparison to the Bardeen and Tersoff–Hamann methods.

Graphical abstract

Keywords

Keywords STM simulation / non-equilibrium Green’s function / Wannier function

1 Introduction

The past decades have witnessed substantial advances and considerable technical progress in the field of STM. Its applications, consequently, have extended far beyond the original function of mapping surface topography, now encompassing a wide range of engineering and analytical areas [1–8]. As the STM tip provides direct access to individual atoms or molecules in real space, one of its most innovative applications is to construct artificial nano-structures by means of “bottom up” fabrication on metal or semiconductor surfaces. The electronic properties of these nano-scale architectures can then be investigated by STM measurements. The major impetus of these studies arises from the ambitious goal of replacing conventional semiconductor devices with future molecular electronics. As the device features today have shrunk to the size of 10–100 nanometers, classical models like the Ohmic conductance relation,

\begin{array}{}\text{(1)}& G=\sigma A/L\end{array}

where G is the conductance, σ the conductivity of a conductor of length L and cross section A start to break down. Consequently, the working principles of desired device characteristics must be based on completely new ideas. This represents a profound scientific challenge. One of the key problems in this field is to obtain a fundamental understanding of the electronic structure and conducting properties of semiconductor based nano-structures. Due to the complexity of the studied systems in terms of both electronic structure analysis and electron transport measurement, the collaboration between experiment and theory becomes necessary to arrive at successful interpretations of experimental observations. An accurate STM simulation provides not only the insight on the origins of observed image contrast, but also sheds light on the energy landscape governing the working principles of device systems.

In this paper, we first discuss to some extent the commonly used Tersoff–Hamann (TH) and Bardeen methods [9–13]. Their characteristic working principles and limitations are presented in relation to actual applications of the models. In practice these simulation methods are implemented using wave functions obtained from finite cluster or slab calculations. This introduces a need for artificial broadening of the discrete energy levels – a technique that can affect the result of the STM simulation. During the past few years, we have therefore developed a more complete, first-principles STM simulation scheme based on non-equilibrium Green’s functions (NEGF), which provides a unified description of the electron transport through the whole STM circuit, consisting of substrate, interface, vacuum, interface, and tip.

The problem with the NEGF transport method as implemented in a number of existing codes (see Ref. [14], and references therein) in the context of STM operation is the need for positioning the STM tip at an arbitrary offset with respect to the surface plane. This together with the fact that tip and sample materials are in general different, makes it highly problematic to represent the combined substrate-tip interface within the same supercell. Therefore, our method deviates from the usual L–C–R setup applied in standard NEGF transport calculations, in that it models the tip and substrate as two separate systems in thermal equilibrium and treats the coupling between these perturbatively. The Green’s function of these semi-infinite systems is obtained once and for all at the beginning of the simulation. We use Density Functional Theory (DFT) to model the electronic structure of tip and substrate, and represent the Hamiltonian, self-energies, and Green’s functions in terms of atom-centered Wannier functions. As the tip is scanned over the surface only a small number of coupling matrix elements has to be calculated which can be easily done on-the-fly.

2 Electron tunneling theory: The traditional approach

In general, a complete single electron transport process in an STM operation can be divided into three elementary steps:

1)Propagation of the electron from the bulk of the sample material to the surface–vacuum interface.

2)Tunneling over the vacuum region between tip and surface.

3)Transport of electron from tip–vacuum interface into the bulk of tip material.

In the simplest case, the electron transport in the first and third step can be regarded as due to free electron modes. The value of tunneling current is then proportional to the probability of finding the surface electron in the tip apex region. If the electron–electron interaction can be ignored, and the potential barrier in the vacuum region has a constant value U, the electron can be described by a single particle Schrödinger equation:

\begin{array}{}\text{(2)}& -{\displaystyle \frac{{\mathrm{\hslash }}^2}{2m}}{\displaystyle \frac{{\mathrm{d}}^2}{\mathrm{d}{z}^2}}\psi (z)+U\psi (z)=E\psi (z)\end{array}

Here, z indicates the distance of electron from the sample surface. In the classically forbidden region (), it is known that the surface states decays exponentially with a constant exponent:

\begin{array}{}\text{(3)}& \varphi (z)=\varphi (0){\mathrm{e}}^{-\kappa z}\end{array}

\kappa =\sqrt{2m(U-E)}/\mathrm{\hslash } is the decay constant of the electron state, and for states at the Fermi level it is related to the square root of the work function. In the Tersoff–Hamann model, the unknown tip electronic structure is approximated by an atomic s orbital which has the same work function as the sample. As a consequence, there is a non-vanishing probability |{\psi }_{s_n}(0){\mathrm{e}}^{-2\kappa z}| for a surface electron from the nth state {\psi }_{s_n} with energy E_{s_n} to tunnel into the corresponding tip state. Assuming that the bias voltage is very small so that the electronic properties of tip and sample are not affected, the tunneling current can be evaluated from the Tersoff–Hamann formula [9, 10]:

where E_{\mathrm{F}} denotes the Fermi level of the sample surface states. In practice, the obtained tunneling current map from the TH model is simply the Local Density of States (LDOS) of the sample at the Fermi level. In particular, the electronic structure of the tip is completely disregarded. Due to its simplicity, the Tersoff–Hamann method has been incorporated into most existing DFT codes. Figure 1(b) shows an STM simulation of a hydrogen terminated Si(001)-(2×1) surface with the Tersoff-Hamann method. Despite the fact that the TH method cannot explain the origin of atomic resolution in STM measurements, it often provides a fairly good qualitative prediction of the experimental observations [15].

In reality, however, the work functions of tip and sam-ple are different, and the decay of the electrons states is not constant. More importantly, a large number of experiments show that the tip effects can play a vital role in the STM imaging process. A typical example is shown in Fig. 2 [16]. In this case the adsorbed perylene molecule on an Ag(110) surface is imaged with a clean tungsten tip and a tungsten tip with an attached molecule. In the first case, perylene molecules appear as four bright lobes, seen in Fig. 2(a). However, if this STM tip had been functionalized by a molecule, the molecular eigenstates have a vanishing overlap with the states of the molecule attached to the STM tip. As a consequence, only the substrate surface states are imaged, while the molecules, due to charge transfer between surface and adsorbed molecules, appear as depressions.

If we assume that the tunneling current will be determined by electron transport in the vacuum region, and make no assumption about the tip orbitals contributing to the tunneling process, then tip effects can be included by the Bardeen approach [12]. In this case, tip and sample surface are assumed to be sufficiently far apart, so that the presence of the tip can be regarded as a small perturbation of the electronic properties of the surface. When the bias voltage is low, the tunneling current can be regarded as the net effect of independent scattering events that transfer electrons across the tunneling barrier. Thus, instead of treating the combined surface-tip system as one quantum mechanical system, the Bardeen approach takes the electronic structures of the two subsystems separately and then calculates the electron transfer between the tip and surface using time dependent perturbation theory. The method is based on four assumptions: (i) the surface–tip interactions are very weak, (ii) the electron–electron interactions can be neglected, (iii) occupation probabilities of the surface and the STM tip are independent from each other and do not change during tunneling and (iv) the eigenstates of the two electrodes are approximately orthogonal. These assumptions allow us to state that the electrons are governed by the following single-electron Hamiltonian:

\begin{array}{}\text{(5)}& H=H_L+H_R+H_T\end{array}

H_T is called the “transfer Hamiltonian”. It describes the elastic transport of electrons from one lead to the other. H_L and H_R are the left and right lead Hamiltonian respectively, they are defined as:

\begin{array}{}\text{(6)}& H_L\psi =-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2m}}{\nabla }^2\psi (\mathit{r})+V_L(\mathit{r})\psi (\mathit{r})\end{array}

and

\begin{array}{}\text{(7)}& H_R\psi =-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2m}}{\nabla }^2\psi (\mathit{r})+V_R(\mathit{r})\psi (\mathit{r})\end{array}

Here, V_L(\mathit{r})=V(\mathit{r}) for \mathit{r} in the left lead region and the vacuum. Correspondingly, V_R(\mathit{r})=V(\mathit{r}) for \mathit{r} in the right lead region and the vacuum. Based on first-order perturbation theory, the probability of elastic tunneling from state {\psi }_{\mu } of the left lead to the right lead is

\begin{array}{}\text{(8)}& P={\displaystyle \frac{2\pi }{\mathrm{\hslash }}}{\displaystyle \sum _{\nu }}{|M_{\mu \nu }|}^2\delta (E_{L,\mu }-E_{R,\nu })\end{array}

where the matrix element M_{\mu \nu } is expressed by

\begin{array}{}\text{(9)}& M_{\mu \nu }=⟨{\chi }_{\nu }|H-H_L|{\psi }_{\mu }⟩\end{array}

At the right lead, H=H_R, we therefore can write

\begin{array}{}\text{(10)}& \begin{array}{cc}{M}_{\mu \nu }\hfill & ={\displaystyle {\int }_R}[{\chi }_{\nu }^{\ast }(H-H_L){\psi }_{\mu }-{\psi }_{\mu }(H-H_R){\chi }_{\nu }^{\ast }]{\mathrm{d}}^3\mathit{r}\hfill \\ \hfill & ={\displaystyle {\int }_R}({\chi }_{\nu }^{\ast }H{\psi }_{\mu }-{\chi }_{\nu }^{\ast }{\psi }_{\mu }{E}_{L\mu }-{\psi }_{\mu }H{\chi }_{\nu }^{\ast }+{\psi }_{\mu }{\chi }_{\nu }^{\ast }{E}_{R\nu })\mathrm{d}{}^3\mathit{r}\hfill \end{array}\end{array}

As we are interested only in elastic transport (E_{R\nu }=E_{L\mu }), the second and last term in Eq. (10) cancel. Keeping in mind that H=H_R at the right lead region, we arrive at

\begin{array}{}\text{(11)}& M_{\mu \nu }={\displaystyle {\int }_R}({\chi }_{\nu }^{\ast }{H}_R{\psi }_{\mu }-{\psi }_{\mu }{H}_R{\chi }_{\nu }^{\ast }){\mathrm{d}}^3\mathit{r}\end{array}

If we take the right lead as the STM tip and substitute H_R with -\frac{{\mathrm{\hslash }}^2}{2m}{\nabla }^2+V(\mathit{r}) in the equation above, the tunneling matrix M_{\mu \nu } can be written as:

\begin{array}{}\text{(12)}& M_{\mu \nu }={\displaystyle \frac{{\mathrm{\hslash }}^2}{2m}}{\displaystyle {\int }_{{\Omega }_T}}({\psi }_{\mu }\nabla {\chi }_{\nu }^{\ast }-{\chi }_{\nu }^{\ast }\nabla {\psi }_{\mu })d\mathit{S}\end{array}

Here, {\Omega }_T represents the enclosed surface surrounding the STM tip. In practice, however, it can be replaced by the separation surface between the two electrodes. Taking occupation probabilities into account, the tunneling current can be written as [12]:

\begin{array}{}\text{(13)}& I={\displaystyle \frac{4\pi e}{\mathrm{\hslash }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}[f(E_{\mathrm{F}}-eV+ϵ)-f(E_{\mathrm{F}}+ϵ)]\cdot {\rho }_S(E_{\mathrm{F}}-eV+ϵ){\rho }_T(E_{\mathrm{F}}+ϵ){|M_{\mu \nu }(E)|}^2\mathrm{d}ϵ\end{array}

where f(E)={[1+\exp (E-E_{\mathrm{F}})/(k_{\mathrm{B}}T)]}^{-1} is the Fermi distribution function while {\rho }_S and {\rho }_T denote the Densities of States (DOS) of the sample and tip respectively.

In Bardeen’s approach the wavefunctions of tip and sample are calculated separately and kept fixed in the following current calculations. In realistic situations, however, as an STM tip moves very close to the sample surface or a large bias voltage is applied, the wave functions of both sub-systems will be distorted and so will their eigenvalues and occupation probabilities. As a result, the application of Eq. (13) is inherently limited to STM simulations of low bias voltage and large tip-sample separations. Despite these limitations, the model has proved to be successful in a wide range of applications and has become a fairly reliable method in interpreting experimental observations.

Let us finally note that the derivation of the Bardeen current formalism [Eq. (13)] is not unique [23]. Based on multiple scattering theory, electron tunneling through a vacuum barrier can be generalized by the Landauer–Büttiker equation [24, 25]:

\begin{array}{}\text{(14)}& I={\displaystyle \frac{4\pi e}{\mathrm{\hslash }}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}dϵ[f({\mu }_S,ϵ)-f({\mu }_T,ϵ)]\cdot \mathrm{Tr}[{\Gamma }_T(ϵ)G^R(ϵ){\Gamma }_S(ϵ)G^A(ϵ)]\end{array}

Here, {\mu }_S and {\mu }_T denote chemical potentials of surface and tip respectively. G^R and G^A are the retarded and advanced Green’s function of the vacuum barrier. {\Gamma }_S,{\Gamma }_T represent the surface and tip contacts, which are defined by the difference between retarded and advanced selfenergies of the corresponding electrodes:

\begin{array}{}\text{(15)}& {\Gamma }_S=\mathrm{i}[{\mathrm{\Sigma }}_S^R-{\mathrm{\Sigma }}_S^A],{\Gamma }_T=\mathrm{i}[{\mathrm{\Sigma }}_T^R-{\mathrm{\Sigma }}_T^A]\end{array}

Palotás and Hofer have demonstrated that Eq. (14) and Eq. (13) are equivalent at zero bias if the Green’s function of the vacuum barrier is approximated to zero order [26]:

\begin{array}{}\text{(16)}& G_{(0)}^{R(A)}({\mathit{r}}_1,{\mathit{r}}_2,ϵ)=G_S^{R(A)}({\mathit{r}}_1,{\mathit{r}}_2,ϵ)+G_T^{R(A)}({\mathit{r}}_1,{\mathit{r}}_2,ϵ)\end{array}

In higher order approximations of the Green’s functions via the Dyson equation tip-sample interactions and multiple electron pathways can also be taken into account [27–30] In this case the tunneling process from the surface to the STM tip can in principle be described with arbitrary precision. However, the electron transport into the leads of the STM circuit is not part of this multiple scattering approach.

3 Complete electron tunneling model for STM simulations

In practical applications of the Bardeen and TH methods the wave functions, overlap matrix elements, and state densities entering Eq. (13) are obtained from DFT calculations where the surface and tip (if the tip is modelled explicitly) are represented as finite clusters or slabs. The fact that the real system contains infinite (macroscopic) electrodes is then modelled by artificially broadening the discrete energy levels to obtain a continuous DOS.

Rather than using wave functions, the transport theory can be formulated in terms of non-equilibrium Green’s functions (NEGF) which allow for inclusion of infinite electrodes thus removing the need for artificial level broadening [31–35]. As discussed in the introduction, the standard implementation of NEGF-transport is, however, unsuited for STM simulations where high efficiency for large and complex systems is crucial. Below we describe how the NEGF formalism can be used in conjunction with a Wannier function basis set to achieve highly efficient and accurate STM simulations taking the bulk nature of the electrodes fully into account.

We begin by making the following assumptions similar to the ones underlying Bardeen and TH theory:

1)The tip and surface are weakly coupled so that the effective potential of the combined tip+surface system can be expressed as the sum of the effective potentials of the separate systems. This assumption also implies that the current can be calculated to lowest order in the tip–surface coupling.

2)Electron–electron interactions in the tunneling process are negligible so that the independent particle approximation is valid.

3)Tip and sample are in their respective thermal equilibrium at all times, i.e., injected hot electrons are instantly thermalized.

4)The voltage drop due to the finite bias is located in the vacuum region and does not affect the coupling matrix elements.

Initially, the tip and sample surface are decoupled and in their thermal equilibria with chemical potentials {\mu }_T and {\mu }_S, respectively. In STM operation, as they are brought close to each other, a net tunneling current starts to flow provided that {\mu }_T\ne {\mu }_S. After a transient period the current flow becomes steady. By introducinga basis set consisting of localized basis functions, the Hamiltonian of the combined system can be written as:

\begin{array}{}\text{(17)}& H=H_S+H_T+H^{(t)}\end{array}

Here, H_S and H_T are the Hamiltonian operators for the tip and sample, while H^{(t)} represents the coupling between the two subsystems. In second quantization they are expressed as:

\begin{array}{}\text{(18)}& H_S={\displaystyle \sum _{i,j\in S}}{H}_{ij}{d}_i{d}_j^{†}\end{array}

\begin{array}{}\text{(19)}& H_T={\displaystyle \sum _{i,j\in T}}{H}_{ij}{c}_i{c}_j^{†}\end{array}

\begin{array}{}\text{(20)}& H^{(t)}={\displaystyle \sum _{i\in T,j\in S}}(V_{ij}{c}_i^{†}{d}_j+V_{ij}^{\ast }{d}_j^{†}{c}_i)\end{array}

The assumption of orthogonality of the basis functions is adopted here to simplify the presentation. In practice, however, most localized basis sets (including the Wannier functions used in the present work) are non-orthogonal. For an extension of the second quantization and NEGF-transport formalism to non-orthogonal basis sets we refer the reader to Ref. [22].

Note that V_{ij} is the coupling matrix element for the full Hamiltonian of the combined system, V_{ij}=⟨{\varphi }_i|H|{\varphi }_j⟩, i.e., including the potential due to both tip and sample. It is equivalent to Bardeen’s transition matrix element and expresses the transition probability from the jth basis function of the surface to the ith basis function of the tip. The current flowing between the two systems is given by the change of the number of electrons in, e.g., the tip, I(t)=-e\frac{\mathrm{d}}{\mathrm{d}t}⟨N_T⟩, where N_T={\sum }_i{c}_i^{†}{c}_i is the number operator. The time derivative can then be evaluated from the equation of motion:

\begin{array}{}\text{(21)}& \begin{array}{cc}I\hfill & =-\mathrm{i}{\displaystyle \frac{e}{\mathrm{\hslash }}}⟨[H,N_T]⟩(t)\hfill \\ \hfill & =\mathrm{i}{\displaystyle \frac{e}{\mathrm{\hslash }}}{\displaystyle \sum _{i,j}}[V_{ji}⟨c_i^{†}(t)d_j(t)⟩-V_{ij}^{\ast }⟨d_j^{†}(t)c_i(t)⟩]\hfill \end{array}\end{array}

From this expression we identify the lesser Green’s functions, and . To lowest order in the tip–surface interaction, V_{ij}, we have the following relation for the (contour-ordered) Green’s function

\begin{array}{}\text{(22)}& G_{ST}(\tau ,{\tau }^{\prime })={\displaystyle {\int }_C}d{\tau }_1{\tilde{G}}_{SS}(\tau ,{\tau }_1)V_{ST}{\tilde{G}}_{TT}({\tau }_1,{\tau }^{\prime })\end{array}

\begin{array}{}\text{(23)}& G_{TS}(\tau ,{\tau }^{\prime })={\displaystyle {\int }_C}d{\tau }_1{\tilde{G}}_{TT}(\tau ,{\tau }_1)V_{TS}{\tilde{G}}_{SS}({\tau }_1,{\tau }^{\prime })\end{array}

where C denotes the Keldysh contour and matrix multiplication in the basis function indices is implied. Here, {\tilde{G}}_{SS} and {\tilde{G}}_{TT} represent the equilibrium Green’s functions of the isolated semi-infinite surface and tip respectively. Applying the Langreth rules for converting the above expressions to the real time we obtain [36]

In steady state, the Green’s function will only depend on the time difference t-t^{\prime }. This allows us to make a Fourier transform of the current formula from the time domain to the frequency domain:

As and in Eq. (24) and Eq. (25) are lesser Green’s functions in equilibrium, they can be expressed in terms of the spectral function [27]:

where the spectral function A=\mathrm{i}[{\tilde{G}}^R-{\tilde{G}}^A]. As a consequence, Eq. (26) can be rewritten as:

\begin{array}{}\text{(29)}& I={\displaystyle \frac{2e}{h}}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}}d\omega [f_S(\omega )-f_T(\omega )]\cdot \mathrm{Tr}[V_{ST}{A}_{TT}(\omega )V_{TS}{A}_{SS}(\omega )]\end{array}

If the tip–sample coupling is weak enough, so that the wavefunctions and eigenvalues of the two subsystems are not distorted, the effects of a finite bias V then enters solely through a shift in energy of the surface eigenstates relative to the tip eigenstates by an amount of eV. With this approximation, the current formula becomes

\begin{array}{}\text{(29)}& I={\displaystyle \frac{2e^2}{h}}{\displaystyle {\int }_{E_{\mathrm{F}}}^{E_{\mathrm{F}}+eV}}\mathrm{Tr}[V_{ST}{A}_{TT}(\omega -eV)V_{TS}{A}_{SS}(\omega )]d\omega \end{array}

Because the spectral function A_{TT} and A_{SS} are independent of the tip–sample configuration, they need to be calculated only once in STM simulations. The retarded Green’s function of, e.g., the tip is given by

\begin{array}{}\text{(30)}& {\tilde{G}}_{TT}^R(\omega )={[(\omega +{\mathrm{i0}}^{+})S_T-H_T-{\mathrm{\Sigma }}^R(ϵ)]}^{-1}\end{array}

where S_T and H_T are the overlap and Hamiltonian matrices of the “active” part of the tip electrode (blue region in Fig. 3) and {\mathrm{\Sigma }}^R is a self-energy describing the coupling to the semi-infinite periodic bulk. The latter can be determined iteratively by the decimation technique [37].

As the tip scans across the sample surface, the current change is determined by the changes in V_{TS}=V_{ST}^{†} which are given by

\begin{array}{}\text{(31)}& V_{TS,ij}=⟨{\varphi }_i^T|-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2m}}{\nabla }^2+V_{\mathrm{eff}}(\mathit{r})|{\varphi }_j^S⟩\end{array}

where the effective potential is

\begin{array}{}\text{(32)}& V_{\mathrm{eff}}(\mathit{r})=V_{\mathrm{eff},S}(\mathit{r})+V_{\mathrm{eff},T}(\mathit{r})\end{array}

and

\begin{array}{}\text{(33)}& V_{\mathrm{eff},\alpha }(\mathit{r})=V_{\mathrm{Hartree},\alpha }(\mathit{r})+V_{\mathrm{ion},\alpha }(\mathit{r})+V_{xc,\alpha }(\mathit{r})\end{array}

with \alpha =S,T.

The derivations of Eq. (30) and Eq. (13) are based on the same approximations, and it is straightforward to show that the two conductance formulas are equivalent. The spectral representations of the lead Hamiltonians are

\begin{array}{}\text{(34)}& H_S={\displaystyle \sum _k}{ϵ}_k^S{d}_k^{†}{d}_k\end{array}

\begin{array}{}\text{(35)}& H_T={\displaystyle \sum _k}{ϵ}_k^T{c}_k^{†}{c}_k\end{array}

Here, {ϵ}_k represents the eigenvalues of the isolated sample or tip Hamiltonians. The spectral functions are also diagonal and are given by [38]:

\begin{array}{}\text{(36)}& A_{SS}(k,\omega )=2\pi \delta (\omega -{ϵ}_k^S)\end{array}

\begin{array}{}\text{(37)}& A_{TT}(k,\omega )=2\pi \delta (\omega -{ϵ}_k^T)\end{array}

Substituting Eq. (36) to Eq. (30) we recover Bardeen’s formula:

\begin{array}{}\text{(38)}& I={\displaystyle \frac{4\pi e^2}{\mathrm{\hslash }}}{V}_{\mathrm{bias}}{\displaystyle \sum _{k\in T,k^{\prime }\in S}}{|V_{k,k^{\prime }}|}^2\delta (E_F-{ϵ}_k^T)\delta (E_F-{ϵ}_{k^{\prime }}^S)\end{array}

When the assumptions of the Tersoff–Hamann model are used, the tunneling matrix element is proportional to the sample wavefunction at tip position z. The equation above then leads to the Tersoff–Hamann current formula:

\begin{array}{}\text{(39)}& I\propto {\displaystyle \sum _{k\in S}}{|{\psi }_k(z)|}^2\delta (E_F-{ϵ}_k^S)\end{array}

As illustrated, the characteristic working principles and limitations of Eq. (30) are identical to those of Bardeen’s approach, but the description has been extended to a complete transport process by taking the coupling to the macroscopic electrodes fully into account. However, it should be noted that the NEGF formalism requires a set of basis functions localized in the transport direction in order to achieve the division into the interaction region and leads, see Fig. 3. At present, the most common choice to fulfill this requirement in first-principles calculations are atomic orbitals. However, recent work by Garcia etal. [39] has demonstrated that a strictly localized basis, optimized from bulk calculations, is insufficient to describe the surface properties of the sample, such as the electronic structure of the surface states and the decay of wavefunctions into the vacuum region. By contrast, first-principles calculations based on plane waves and a grid basis are generally more accurate, but they do not lend themselves directly as basis sets for transport calculations. An ideal approach is to generate a set of localized orbitals (Wannier functions) from the extended eigenstates obtained from a plane wave or grid based calculation [18–21]. In this way accurate transport calculations can be carried out at very low computational cost [14, 17].

4 Applications

4.1 Experimental background

STM has been used for the study of buried dopants in silicon [40–43]. However, STM lithography [44] also allows the possibility to position dopant atoms in silicon with near atomic precision [45]. Recent studies have adopted this technique to the fabrication of atomically controlled dopant devices in silicon [8, 46, 47]. In these experiments, the clean Si(001)-(2×1) surface is exposed to atomic hydrogen so that all surface dangling bonds are saturated by the formation of Si–H bonds, see Fig. 1(a). In a second step, the STM tip is employed to selectively desorb hydrogen atoms from the monohydride Si(001) surface. The recreated Si dangling bonds then react with dosed phosphine (PH₃) from the gas phase, leaving the surrounding hydrogen terminated area intact. Annealing the surface at around 350℃ results in the incorporation of P atoms from phosphine into the Si subsurface. The precision lithography of 3 adjacent dimers individual for placement of P atoms has also been demonstrated previously [45]. Once dosed and annealed to incorporate the dopants, the surfaces are then encapsulated with epitaxial silicion to activate the dopant.

The electronic property of the buried devices can then be investigated with STM techniques [43]. Because the calculation of the tunneling current for this kind of systems should include a complete electron transportprocess (semi-infinite substrate-interface-vacuum-interface-semi-infinite tip), the prevalent Bardeen approach becomes unreliable.

In order to demonstrate the potential of this new STM simulation scheme, we take a single P-doped Si(001)-(2×1):H surface as a test system. Figure 5 shows empty states and filled states STM imaging of the same area of the sample [43]. The P dopant is known to be buried several atomic layers below the surface. As can be seen, the surface area above the P atom is imaged as a rectangular protrusion at all biases. The bright peak extends over two dimer rows in the [110] direction and about three dimers in the [\overline{1}10] direction. This feature has been frequently observed for buried dopant atoms beneath Si(001)-(2×1):H surfaces, and can therefore be used as a benchmark to identify the location of the dopant atoms. Interestingly, the apparent height of this bright feature decreases as the voltage value increases from 1.0 V to 1.6 V and fade out at –2 V bias. An STM simulation of this bias-dependence of the imaging contrast is thus a suitable challenge to the new method.

4.2 DFT calculations

The DFT calculations for the Si surface are performed with the Vienna ab-initio Simulation Package (VASP) [49–52] and the Grid-based Projector-Augmented Wave (GPAW) code [53, 54]. STM simulations using the TH and Bardeen approach are carried out using the BSKAN code, which is described in Ref. [55]. In these calculations the LDOS and overlap of vacuum wavefunctions of the tip and sample at the separation surface are calculated numerically based on the VASP results. The electronic structure calculations used as input to the new NEGF-STM scheme are performed with GPAW. In all calculations we use the LDA exchange-correlation functional. We use a 400 eV plane wave cut off in the VASP calculations and a grid spacing of 0.2 \mathring{\mathrm{A}} in the GPAW calculations. Because of the large size of the system the first Brillouin zone was sampled at the \Gamma point only.

Figure 6(a) shows the supercell used to represent the P-doped Si film. The P atom is placed 6.6\mathring{\mathrm{A}} below the top H atoms. The top three layers of Si atoms and all H atoms were fully relaxed until the residual force was smaller than 0.01\mathrm{eV}/\mathring{\mathrm{A}}. The non-periodic part of the semi-inifite surface lead is illustrated by Fig. 6. A couple of convergence layers of silicion atoms (in red color) are included to assure the convergence of the potential inside the surface. All the atoms above convergence layers are referred as the non-periodic part of the surface lead (see Fig. 3). The convergence layers are then cut out and the surface atoms are extended with an infinite number of principal layers from silicon bulk calculations. The STM tip is modeled by a W(111) film composed of three periodic layers along the [111] direction and a small pyramid attached to one side, see Fig. 7. Atomic relaxations in the tip calculation are limited to atoms in the surface layer.

4.3 Bardeen and Tersoff–Hamann models

The STM simulation using the TH model is shown in Fig. 8. The position of the buried P atom is marked with a blue circle. At short tip–sample distances (around 3.5\mathring{\mathrm{A}}), the TH model leads to surface topographies at all biases [Fig. 8(a)–(c)]. However, when the tip-sample distance increases to about 5.5\mathring{\mathrm{A}}, the dopant induced modification on surface states becomes visible [Fig. 8(d)–(f)]. At +1.0 V, the H atom above the phosphorus appears as a small peak. Its surrounding area is shown as depressions; however, the dimers at the corners of the unit cell are imaged as the highest protrusions. When the bias changes to +1.4 V, we see that the symmetry of the contrast changes dramatically. The surface area directly above the P atom is now the highest peak. Its neighbouring Si-H from the same silicon dimer is imaged as a hole while the adjacent dimer row appears a valley. At –2.0 Vvoltage, however, the contrast of the calculated STM image remain unchanged. Figure 9 shows STM calculations at the same biases (+1.0 V, +1.4 V and -2.0 V) with the Bardeen model. Not surprisingly, the obtained images qualitatively agree with the TH results at larger tip-sample distances. However, the simulated STM images from the TH and Bardeen methods do not fit the experimental observations, as the surface atoms near the blue circle should appear as bright protrusions at all these biases.

From the DFT calculations we know that the dopant state is in the band gap of the Si(001)-(2\times 1):H surface positioned about 0.3 eV below the conduction band edge. The band decomposed charge density calculation allows a view of the dopant states in the unit cell and quite helpful in understanding the calculated STM images:

\begin{array}{}\text{(40)}& \rho (\mathit{r})={\displaystyle \sum ^N}{|\psi (n,\mathit{r})|}^2\end{array}

Here, N represents the number of bands in the selected energy window and n denotes band index. Figure 10(a) and (b) show the calculated band decomposed charge density of the P states in the band gap of undoped Si(001)-(2\times 1):H. As can be seen, the states of the extra electron from P dopant are well confined below the vacuum–sample interface and therefore do not directly contribute to the Bardeen’s tunneling current. Figure 11(a) shows the site Projected Density of States (PDOS) of top-most Si and H atoms of the doped and undoped sample. The Fermi level of the undoped surface is used as a reference. As can been seen, the PDOS of the Si and H atoms on the doped surface have some states in the energy window of 0.2–1.5\mathrm{eV} relative to the Fermi level. These states arise because of the long range effects of the P dopant. The STM contrast obtained from the Bardeen method in Fig. 9 is thus attributed to the indirect effect of the buried P atom. The mismatch between STM observation and the Bardeen model arises from the incomplete calculation of the electron transport process. As shown in Fig. 11(b), the dopant atom has a much stronger effect on the Si atoms close to it. An immediate and important consequence is that when the electron propagates from surface into bulk, the scattering probability is strongly dependent on the local environment below the surface. The calculation of electron transport inside the sample is thus the dominating step for evaluating the tunneling current.

4.4 NEGF-STM simulations

In this section we apply the NEGF-STM method to the P-doped Si(001)-(2\times1):H system. In order to fulfill the basic approximation of weak sample–tip coupling, the apex of the W(111) tip was always kept at least 5\mathring{\mathrm{A}} above the top H atoms. The atom-centered Wannier functions are generated from Bloch states of GPAW DFT calculations. All the eigenstates below an energy of E_{\mathrm{F}} + 3.5 eV are included in the construction of the Wannier basis. The locality of the Wannier functions is exploited by representing each basis function in a small box of dimensions 10\times 10\times 10\mathring{\mathrm{A}}, see Fig. 12. Outside the box the Wannier function is set to zero. We have checked that the results are converged with respect to the size of the localization box. The Hamiltonian and overlap matrix elements of two orbitals can be evaluated by considering only the overlapping area of the corresponding localization boxes, see Fig. 12. As the tip scans over the surface only a small number of matrix elements, corresponding to localization boxes with non-zero overlap, needs to be calculated. The Green’s functions of the tip and sample, {\tilde{G}}_{TT} and {\tilde{G}}_{SS}, are calculated at the beginning of the calculation using standard techniques, see e.g. Ref. [14].

The calculated STM images are shown in Fig. 13. As can be seen, in the doped system, the surface atoms above the P atom appear as protrusions at all biases, marked with the red boxes in Fig. 13(a)–(d). As our unit cell includes only two dimers rows and four dimers in each dimer row, the contrast of the characteristic peak is not very obvious. However, comparing the STM simulations of undoped Si(001)-(2×1):H, seen in Fig. 13(e)–(h), the doped surfaces appear brighter at all biases. The image contrasts inside the red boxes in Fig. 13(a)–(d) are smaller than those in the undoped Si(001)-(2×1):H surface. This is also consistent with experimental results. Of importance, the increased image contrast has also been seen at the edges of doped surface unit cell, indicating the calculated characteristic protrusion is of the same dimension as observed experimentally. Figure 13(i) shows a height-profile along the line indicated in Fig. 13(a). The value of apparent height in the line-profile is obtained by subtracting the average height of clean surfaces at the same biases. The decrease of apparent height with increasing bias voltage from +1.0 V to +1.6 V has also been reproduced. In addition, Fig. 13(i) shows that at –2.0 V, the apparent height of the bright feature is almost the same as for the undoped surface. This also agrees well with the experimental observations.

In order to understand the origin of the characteristic protrusion for Si encapsulated dopant atoms, we performed band decomposed charge density calculations of the P-doped Si(001)-(2×1):H system at the energy window from 0 eV to +1.0 eV, +1.6 eV and –2.0 eV relative to the surface Fermi level. As shown in Fig. 14, at small biases, the dopant states provide the majority of states for the tunneling current inside the surface. This enhancement is imaged as a protrusion in STM measurements. However, as the absolute value of bias voltage increases, more and more Si states get involved in the electron transport, the contribution from dopant states becomes less important. Correspondingly, the apparent height of the characteristic protrusion decreases. For the characteristic protrusions, both experiments and theory show that one dimer row appears higher than the other. This is to some extent due to the relaxation of surface atoms. However, PDOS calculations reveal that the electronic properties of these surface atoms are not identical either. The electronic effects can therefore also be important in interpreting the contrast difference of the two dimer rows. As we have shown, the calculation of STM images based on non-equilibrium Green’s functions reproduces all features of the experimental data. This new STM simulation scheme is thus capable of addressing electron transport in STM junctions comprehensively and even allows to address the hitherto difficult to analyze transport on doped semiconductor surfaces.

5 Summary

We have presented an accurate and efficient STM simulation method based on non-equilibrium Green’s functions (NEGF) and a Wannier function basis set. The scheme has been compared to the more traditional STM schemes, i.e., the Tersoff–Hamann and Bardeen’s models, with particular focus on doped semiconductor systems. The main novelty of the NEGF-STM method is the proper account of the semi-infinite (macroscopic) electrodes which provides the correct continuous spectrum at the surface thus removing the need for an artificial level broadening inherent in practical applications of the Tersoff-Hamann and Bardeen’s models. In order to show the difference between this new method and the Bardeen model, we performed STM simulations on a single P-doped Si(001)-(2×1):H surface. In this system, the important transport region that determines the STM contrast is not the separation vacuum but the interface below the sample surface. Consequently, the Bardeen model fails to reproduce the experimental observations. The new STM simulation scheme, by contrast, evaluates electron transport in the whole circuit. As a result, the simulated STM images completely agree with the bias dependent experimental data. The scheme, introduced in this paper, thus substantially improves the state of the art in STM theory, and should become the de facto standard of simulations within the next years.

References

[1] G. Binnig, H. Rohrer, Ch. Gerber, and E. Weibel, Phys. Rev. Lett., 1982, 49: 57

[2] G. Binnig, H. Rohrer, Ch. Gerber, and E. Weibel, Phys. Rev. Lett., 1983, 50: 120

[3] J. A. Heinrich, C. P. Lutz, J. A. Gupta, and D. M. Eigler, Science, 2002, 298: 1381

[4] C. Joachim, J. K. Gimzewski, and A. Aviram, Nature, 2000, 408: 541

[5] H. C. Manoharan, C. P. Lutz, and D. M. Eigler, Nature, 2000, 403: 512

[6] S.-W. Hla, L. Bartels, G. Meyer, and K.-H. Rieder, Phys. Rev. Lett., 2000, 85: 2777

[7] S.-W. Hla and K.-H. Rieder, Ann. Rev. Phys. Chem., 2003, 54: 307

[8] M. Fuechsle, S. Mahapatra, F. A. Zwanenburg, M. Friesen, M. A. Eriksson, and M. Y. Simmons, Nature Nanotechnology, 2010, 5: 502

[9] J. Tersoff and D. R. Hamann, Phys. Rev. B, 1981, 31: 805

[10] J. Tersoff and D. R. Hamann, Phys. Rev. Lett., 1985, 50: 1988

[11] W. A. Hofer, G. Ritz, W. Hebenstreit, M. Schmid, P. Varga, J. Redinger, and R. Podloucky, Surf. Sci. Lett., 1998, 405: L514

[12] J. Bardeen, Phys. Rev. Lett., 1961, 6: 57

[13] W. A. Hofer and J. Redinger, Surf. Sci., 2000, 447: 51

[14] K. S. Thygesen and K. W. Jacobsen, Chem. Phys., 2005, 319: 111

[15] W. A. Hofer, A. S. Foster, and A. L. Shluger, Rev. Mod. Phys., 2003, 75: 1287

[16] Z. T. Deng, H. Lin, W. Ji, L. Gao, X. Lin, Z. H. Cheng, X. B. He, J. L. Lu, D. X. Shi, W. A. Hofer, and H. J. Gao, Phys. Rev. Lett., 2006, 96: 156102

[17] A. Calzolari, N. Marzari, I. Souza, and M. B. Nardelli, Phys. Rev. B, 2004, 69: 035108

[18] G. H. Wannier, Phys. Rev., 1937, 52: 191

[19] N. Marzari and D. Vanderbilt, Phys. Rev. B, 1997, 56: 12847

[20] K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen, Phys. Rev. Lett., 2005, 94: 026405

[21] K. S. Thygesen, L. B. Hansen, and K. W. Jacobsen, Phys. Rev. B, 2005, 72: 125119

[22] K. S. Thygesen, Phys. Rev. B, 2006, 73: 035309

[23] C. J. Chen, Introduction to Scanning Tunnelling Microscopy, New York: Oxford University Press, 1993

[24] S. Datta, Electronic Transport in Mesoscopic Systems, Cambridge: Cambridge University Press, 1995

[25] M. Büttiker, Y. Imry, R. Landauer, and S. Pinhas, Phys. Rev. B, 1985, 31: 6207

[26] K. Palotás and W. A. Hofer, J. Phys.: Condens. Matter, 2005, 17: 2705

[27] A. S. Foster and W. A. Hofer, Scanning Probe Microscopy, Spring Series in NanoScience and Technology, Springer, 2006

[28] W. A. Hofer and A. J. Fisher, Phys. Rev. Lett., 2003, 91: 036803

[29] W. A. Hofer and A. Garcia-Lekue, Phys. Rev. B, 2005, 71: 085401

[30] W. A. Hofer, A. Garcia-Lekue, and H. Brune, Chem. Phys. Lett., 2004, 397: 354

[31] C. Caroli, R. Combescot, P. Nozieres, and D. Saint-James, Journal of Physics C, 1971, 4: 916

[32] T. E. Feuchtwang, Phys. Rev. B, 1974, B10: 4135

[33] T. E. Feuchtwang, Phys. Rev. B, 1974, 10: 4121

[34] T. E. Feuchtwang, Phys. Rev. B, 1976, 13: 517

[35] Y. Meir and N. S. Wingreen, Phys. Rev. Lett., 1992, 68: 2512

[36] H. Hauge and A. P. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors, Springer Series in Solid-State Physics, Springer, 1996

[37] F. Flores, F. Guinea, C. Tejedor, and E. Louis, Phys. Rev. B, 1983, 28: 4397

[38] K. Flensberg and H. Bruus, Many-Body Quantum Theory in Condensed Matter Physics, Chapter 8, New York: Oxford University Press, 2004

[39] S. Garcia-Gil, A. Garcia, N. Lorente, and P. Ordejon, Phys. Rev. B, 2009, 79: 075441

[40] L. Liu, J. Yu, and J. W. Lyding, Appl. Phys. Lett., 2001, 78: 386

[41] L. Liu, J. Yu, and J. W. Lyding, IEEE Trans. Nanotechnol., 2002, 1: 176

[42] G. W. Brown, H. Grube, and M. E. Hawley, Phys. Rev. B, 2004, 70: 121301

[43] L. Oberbeck, N. J. Curson, T. Hallam, M. Y. Simmons, and R. G. Clark, Thin Solid Films, 2004, 464: 23

[44] J. W. Lyding, T. C. Shen, J. S. Hubacek, J. R. Tucker, and G. C. Abeln, Appl. Phys. Lett., 1994, 64: 2010

[45] S. R. Schofield, N. J. Curson, M. Y. Simmons, F. J. Rueβ, T. Hallam, L. Oberbeck, and R. G. Clark, Phys. Rev. Lett., 2003, 91: 136104

[46] F. J. Ruess, L. Oberbeck, M. Y. Simmons, K. E. J. Goh, A. R. Hamilton, T. Hallam, S. R. Schofield, N. J. Curson, and R. G. Clark, Nano Lett., 2004, 4: 1969

[47] A. Fuhrer, M. Fchsle, T. C. G. Reusch, B. Weber, and M. Y. Simmons, Nano Lett., 2009, 9: 707

[48] J. L. O’Brien, S. R. Schofield, M. Y. Simmons, R. G. Clark, A. S. Dzurak, N. J. Curson, B. E. Kane, N. S. McAlpine, M. E. Hawley, and G. W. Brown, Phys. Rev. B, 2001, 64: 161401(R)

[49] G. Kresse and J. Hafner, Phys. Rev. B, 1993, 47: 558

[50] G. Kresse and J. Hafner, Phys. Rev. B, 1994, 49: 14251

[51] G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6: 15

[52] G. Kresse and J. Furthmüller, Phys. Rev. B, 1996, 54: 11169

[53] J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen, Phys. Rev. B, 2005, 71: 035109

[54] J. Enkovaara, et al., J. Phys.: Condens. Matter (in press)

[55] W. A. Hofer, Progr. Surf. Sci., 2003, 71: 147

Rights & permissions

©Higher Education Press and Springer-Verlag Berlin Heidelberg 2010