Scanning tunneling spectroscopy of unoccupied surface

resonances at free-electron-like metal surfaces.

T. Fondén^{†}, S. Papadia^{∗}, and M. Persson^{§}

Department of Applied Physics, Chalmers University of Technology, S-412 96 Göteborg, Sweden

Abstract

We extend and detail a previously developed model for formation of electronic resonances at free-electron-like metal surfaces, in order to calculate scanning tunneling (ST) spectra. The effect of the tip is mimicked by inclusion of an external field, self-consistently, in a jellium description of the surface potential. The lattice-induced corrugation of the potential is included perturbatively via a pseudopotential. We compare our calculated ST spectra for Al(111) with experimental spectra for that surface and conclude that a peak occuring below the metal vacuum level is a “crystal-derived” resonance, in the sense that lattice effects are crucial for its manifestation.

SHORT TITLE: STS of resonant electron states at metal surfaces. PACS numbers: 73.20.At, 73.40.Gk

Typeset using REVTEX

I. INTRODUCTION

The investigation of unoccupied electron states at surfaces is crucial for our understanding of phenomena taking place at surfaces, such as scattering of electrons, atoms and ions [1–4] and vacuum tunneling [5–8]. A particular type of such states is the image-potential-induced series, produced when an electron is trapped in a potential well defined by the image-potential part of the surface barrier on the vacuum side, and a projected band gap on the crystal side [9]. Experimentally, image-potential states have been studied mainly by lowenergy electron diffraction (LEED) and inverse photoemission spectroscopy (IPE) [10–12] and preferably at surfaces where a projected band gap exists at or close to the energy of the series. Theoretically, the existence of a high substrate reflectivity, as e.g. given by the presence of a band gap, has been crucial for the building up of an image state as described in multiple-scattering models [1,2]. At energy ranges outside a metal band gap, the image states broaden into resonances and eventually cease to be resolvable by experimental means.

There are though experimental evidence and theoretical support that resonances persist at energies surprisingly far from band gaps. IPE spectra [12] of a free-electron like metal surface with no projected band gaps in the probed region – Al(111) – revealed a spectral feature close to the vacuum level which was suggested to be an image resonance. Several authors provided calculations in favor of that explanation [13–15] with different viewpoints. The interpretation of IPE features close to the vacuum level in terms of the surface density of states (SDOS) was questioned by Schaich and Lee [16] who calculated spectra directly and demonstrated that, even for a pure jellium potential with no SDOS features, the IPE spectra showed several features below the vacuum level. These structures in IPE spectra were interpreted by Schaich and Lee as due to a matrix element effect, similar to the threshold effect in photoemission, with no direct coupling to the SDOS of the metal. Thus, the interpretation of IPE features close to the vacuum level in terms of SDOS was questioned. The calculations by Schaich and Lee [16], motivated further complementary studies of the unoccupied electronic states of free-electron-like metals.

The scanning tunneling spectroscopy (STS) technique has been suggested as a possibility of obtaining information about the surface electronic structure [17–20]. This spectroscopy has been proven to be able to give information for surface states on Graphite [19] and Si [20]. It has also been asserted that STS provides identification of the image state spectrum on metals [5,21]. In such experiments, the normalized conductance is recorded, often in a constant-current mode where the bias is varied under the condition of keeping the tunnel current constant. These curves are then related to the electronic states of the substrate within different theories.

We will in the following present a model for interpretation of conductance curves, which is an extension of the model in Ref. [14], and we demonstrate how peaks in the STS curves are connected to SDOS features of the substrate. Application of the model on Al(111) and comparison between STS and theoretical calculations demonstrate several resonances in the SDOS. In particular, there is one resonance below the metal vacuum level which shows up in STS because of lattice effects – it is a “crystal-derived” resonance. From extrapolation to the field-free case (no STM tip present), we can identify this resonance as a “crystalderived” image resonance. In previous analyses by us [14,22], we assigned this resonance to a feature seen in IPE spectra from the close agreement in their energy positions. Our later analysis has revealed an error in the evaluated phase shift within our model. In this paper, we correct the phase shift. The consequences are that this resonance appears at a lower energy with corresponding corrections of its position in STS spectra, but the basic conclusion that without a resonance there is no characteristic tunneling feature below the vacuum energy is still valid. Finally, the crystal reflectivities obtained in our model and in a two-band model [15] are compared and it is found that the two-band model underestimates the crystal reflectivity at energies far from the considered band gap.

The paper is organized as follows. In Sec. II, we present the different steps in our model. First, the construction of the potential between tip and substrate in the tunneling geometry is described in subsection (II A); the effects of the scattering from the ion-core lattice are included in the expressions for the current density and in the expression for the substrate density of states, described in subsections II C and II D, respectively. We present and discuss our results in Sec. III by first comparing our calculated conductance curves with our calculated SDOS, obtained for different initial conditions in experiments; then we compare the calculated conductance to experimental conductance curves for Al(111). A summary is given in Sec. IV. In the Appendix, we present the detailed derivation of the crystal reflectivity and a comparison is made to the reflectivity obtained in a two-band model.

II. MODEL

Our calculations are based on a quasi-one-dimensional model where the starting point for the modelling of both the sample and the tip is the jellium model for a metal surface. This is motivated by the use of rather blunt tips in STS-experiments. The calculated current densities in such a one-dimensional model should thus give a fair description of the actual experimental situation. In the following, we describe our model for the tip-sample interface potential and the calculation of the reflectivity from the sample crystal potential in detail. In addition we show how the tunneling current and the surface-induced DOS are calculated using the obtained interface potential and crystal reflectivity. If not explicitly stated otherwise, Hartree atomic units are used throughout where h̄ = m = e = 1 (m is the electron mass).

A. Potential construction

We model the interface potential for the scanning tunneling experiment on Al(111) for large initial separations between the tip and sample of typically 15 - 25 Å. In this situation, we expect to have a linear region in the interface potential for a constant external electric field F. The interface potential, V_{I}(z), may then be decomposed into a tip and a sample part which are treated separately. This decomposition is not justified at short distances, e.g., below 5 Å where the barrier do not have a linear part, as can be seen in the work by Lang [23]. The potentials, V_{T} (z) and V_{S}(z) of the isolated tip and sample, respectively, are both modelled by the potential V_{0}(z; ±F_{ext}) of a semi-infinite jellium (z → −∞ in the bulk) in external uniform electric fields ±F_{ext}, with an electron-gas density parameter, r_{s} = 2.07a_{0}, appropriate for aluminium. The potential V_{0}(z; F_{ext}) is calculated self-consistently, in the density functional scheme [24] as applied by Schreier and Rebentrost [25], using the boundary condition that the electric field F should be uniform and equal to 2F_{ext} far from the surface. The details of the approximation for the exchange-correlation part of the potential are given at the end of this subsection. For the sample, V_{S}(z) = V_{0}(z; F_{ext}) gives rise to a repulsive and confining potential while for the tip, V_{T} (z) = V_{0}(z; −F_{ext}) gives rise to a barrier that the electrons can tunnel through. The interface potential is now obtained by matching V_{S}(z) and V_{T} (z) at a point z_{m} in the linear asymptotic region as,

V_{I}(z) = \{ \begin{matrix}V_{S}(z), & z < z_{m} \\ V_{T} (s − z) + V + Φ_{T} − Φ_{S} & z > z_{m} \end{matrix}  (2.1)

where Φ_{S} and Φ_{T} are the work functions of the substrate and tip, respectively, V is the bias, and the tip-sample distance s, measured relative to the jellium edges of the tip and sample respectively, is determined by the condition that V_{I}(z) should be continuous at z_{m}. This matching gives rise to an implicit relation between the distance s, the field F and the bias V. In the particular experimental setup, the tip is made of Ir that has a work function of ∼ 1.5 eV larger that Al, but it is very likely that it is covered by Al atoms [26]. As a consequence, the contact potential in the experiments may vary between 0 and 1.5 eV. In the discussions, we will use the lower limit as an example, corresponding to an Al tip, but we will , when comparing to experiments, also present results that account for the maximum work function difference of 1.5 eV, corresponding to an Ir tip. This is done by choosing Φ_{T} − Φ_{S} = 1.5 eV in Eq. (2.1) [27]. We will see that the important issue is to have a good model of the substrate. A typical interface potential shape obtained as in Eq. (2.1), with Φ_{S} = Φ_{T}, is shown in Fig. 1. Note that we have an extended linear part in the interface potential due to the relatively low bias V = 4.0 eV and large tip-surface separation s=17 Å. The electric field F = 0.20 eV/Å, is close to what is expected from an assumption based on an extended region of a uniform electric field between the tip and the sample, V/s= 0.24 eV/Å.

In order to model the (static) image-barrier effects, it is necessary to include the longrange image potential of a metal surface. We proceed as proposed by Serena, Soler, and Garcia [28] as follows. The local density approximation (LDA) for the exchange-correlation potential V_{xc} is used in the bulk and up to the location of the image plane. Beyond that, a nonlocal form for the potential is used, of the form

V_{xc}^{NL}(z) = −\frac{1 − [1 + b(z − z_{IP})/4]e^{−b(z−z_{IP})}}{4(z − z_{IP})} ,  (2.2)

where z_{IP} is the location of the image plane and b is a parameter chosen to give the LDA value for V_{xc}^{NL} at z_{IP}. The location of the image plane is also calculated self-consistently in that the position of the centroid of the excess charge, due to the applied external field, is taken as z_{IP} at every iteration. When self-consistency is reached, the value of z_{IP} is stable. This approach yields values for z_{IP} that are essentially the same as Lang and Kohn’s calculated values using linear response of a jellium to a test charge [29]. If exchange effects are important, there will be a modification to the results from this treatment [30]. As demonstrated by Eguiluz et. al. [30], a “Kohn-Sham”, electron “feels” an image plane location closer to the surface than a test charge does (z_{0} = 0.72 a.u. instead of 1.46 a.u for r_{s} = 2.07). These effects are even more pronounced for larger r_{s} values. This suggests that there are still improvements to be made on our model for the surface barrier potential.

B. Inclusion of sample ion-cores

The scattering of the electrons from the Al-sample ion-cores is handled perturbatively using pseudopotential theory. We confine ourselves to the situation where we are in an energy range far from band gaps. Such a situation and the fact that Al can be represented by a weak pseudopotential, suggest a treatment of the scattering of the electrons from the ion-cores using the Born approximation. The detailed derivation of the reflectivity, r_{c}, in this approximation is presented in the Appendix, where we also discuss the connection between our approach and the standard two-band model. The resulting Born expression for the reflectivity of a semi-infinite crystal measured relative to the first lattice plane (see Fig. 2) is given by,

r_{c}(k′, k) = −i\frac{v(k′ − k)}{|a × b| |k_{z}| (1 − e^{i(k−k′)·c})} .  (2.3)

Here k and k′ are the wavevectors of the incident and scattered wave, respectively, v(k′ − k) is the form factor, i.e the Fourier transform of the atomic pseudopotential, a and b are the surface basis vectors, and c connects an atom to a neighbouring atom in the next atomic layer further into the crystal. In the case of Al(111), the area of the surface unit cell is given by, |a × b| = \sqrt{3}a^{2}/4 and the distance [31] d between crystal planes parallel to the surface is d = a/\sqrt{3}, where a=7.7 a.u is the bulk lattice constant of Al. The surface periodicity gives the restriction k′_{∥} = k_{∥} + g_{∥}, where g_{∥} is a surface reciprocal lattice vector. There is no such restriction for the momentum transfer perpendicular to the surface. Note that the result for r_{c} in Eq. (2.3) neglects possible modifications of the form factor for the atoms in the surface region and also neglects the possibility of (geometrical) surface relaxation [32]. Furthermore, the behaviour of the denominator with momentum transfer explicitly shows that this expression breaks down when the transfer comes close to a Brillouin zone boundary where a band gap opens up.

In order to estimate the reflectivity, we have chosen a pseudopotential form factor which was originally fitted to the measured phonon dispersion of bulk aluminum [33]. The justification for the use of such a pseudopotential is that we are sufficiently close in energy to the Fermi level. This pseudopotential reproduces accurately the width of the band gap in the (111) direction located about 5 eV below ϵ_{F} [34]. In the tunneling situation we find that the effective tunneling occurs from states with perpendicular energy close to the Fermi level, i.e., k_{∥} ∼ 0. Considering then specular reflection of an electron propagating normal to the surface, k_{∥} = k′_{∥} = 0 and −k_{z} = k′_{z} = k, the pseudopotential form factor does not vary much for energies around the vacuum level ϵ − ϵ_{F} ∼ 4 eV, where 2k/k_{F} ≈ 2.3. For this situation, we have from Eq. (2.3) that,

\begin{aligned}|r_{c}| & = & \frac{|v(2k)|}{2|a × b| |sin(kd)|} \\ ϕ_{c} & = & π . \end{aligned}  (2.4)

Here, and in the following, the phase shift ϕ_{c} is taken with respect to the jellium edge z_{0} and includes thus an extra phase due to the electron motion between the first crystal layer and the jellium edge. In previous analyses [14,22], ϕ_{c} was erroneously evaluated to ∼ 4π, which has consequences on the energy of the resonances but not on their existence. The value of |r_{c}| is around 0.1 as obtained from (2.4), while the actual value of ϕ_{c} is determined from the sign of v(2k)/sin(kd); v(2k) is positive and sin(kd) is negative in the described situation. The value |r_{c}| = 0.1 will be used throughout in our calculations. Finally, we would like to stress, as shown in the Appendix, that a simple two-band model underestimates |r_{c}| in an energy region far away from a band gap.

The range of applicability of the result in Eq. (2.3) is not restricted to aluminium but can be applied to other free-electron-like materials. The condition that has to be fulfilled in order for the Born expression to be valid is that we apply it in an energy range sufficiently far from projected band gaps. However, if the effect of the crystal reflectivity is to be noticeable, we need a reflectivity comparable to the one for aluminium. We find similar and even larger reflectivities e.g for Pb and Ba. In the case of Pb, there is a band gap with a width of ∼ 3 eV at 3.5 eV below the Fermi energy, so it may be out of the range of the application of the Born approximation. In the case of Ba(100), we have equally favourable conditions as in the case of Al(111) so there is as chance to experimentally resolve “crystal-derived” image resonances on that surface.

C. Tunneling current

We describe in the following the calculation of the tunneling current and the conductance in our interface potential and show how the effects of the sample ion-cores can be accounted for in these calculations by using the calculated crystal reflectivity r_{c}.

The tunneling current density J_{0} for two free-electron systems in a planar geometry is given by the expression [35],

J_{0} = \frac{1}{2π^{2}} [∫_{ϵ_{F T} −V}^{ϵ_{F T}} dϵ_{z}(ϵ_{F T} − ϵ_{z})|T_{0}(ϵ_{z})|^{2} + V ∫_{U_{T}}^{ϵ_{F T} −V} dϵ_{z}|T_{0}(ϵ_{z})|^{2}]  (2.5)

where V is the bias, T_{0}(ϵ_{z}) is the transmission amplitude of electrons through the barrier with an energy, ϵ_{z}, normal to the surface. Further, ϵ_{F T} = ϵ_{F S} + V and ϵ_{F S} are the Fermi energies of the tip and the sample, respectively, and U_{T} is the energy of the tip band bottom. All energies are measured relative to the vacuum level, ϵ_{vac} = ϵ_{F S} + Φ_{S}, of the field-free sample, where Φ_{S} is the work function of the sample. The first term gives the contribution to J_{0} from tip electrons tunneling into unoccupied states of the sample and the second term the net contribution to J_{0} from tip and sample electrons tunneling between occupied states. Before we show how the effects of the crystal potential influence the tunneling current, we will first show that, in the situation of interest, the dominant contribution to J_{0} comes from electrons with ϵ_{z} within a few tenths of an eV from ϵ_{F T} (k_{∥} ∼ 0). In particular the contribution from the second term in (2.5) to J_{0} is negligible.

The transmission amplitude, T_{0}(ϵ_{z}), has been determined from the electron wavefunctions in V_{I}(z) by a direct numerical integration of the one-electron Schrödinger equation based on a modified midpoint method and Richardson extrapolation [36]. The boundary conditions in a potential like in Fig. 1 for electron wave functions through the barrier are

Ψ(z) = \{ \begin{matrix}e^{−ik_{1}z} + r(ϵ_{z})e^{ik_{1}z} & for z → ∞ \\ t(ϵ_{z})e^{−ik_{2}z} & for z → −∞ \end{matrix}  (2.6)

where r(ϵ_{z}) is the amplitude of the reflected wave on the tip side and t(ϵ_{z}) the amplitude of the transmitted wave on the substrate side. We have that the momentum k_{1} = \sqrt{2(ϵ_{z} − U_{T})} and k_{2} = \sqrt{2(ϵ_{z} − U_{S})} where U_{S} is the energy of the substrate band bottom. The transmission coefficient is then given by |T_{0}(ϵ_{z})|^{2} = \frac{k_{2}}{k_{1}} |t|^{2}.

At the tip-sample distances, s, and biases, V, appropriate for the scanning tunneling experiment, we obtain typically very small transmission probabilities for the tip-electrons. For instance, under the conditions prevailing in Fig. 1, we find that |T_{0}(ϵ_{F T})|^{2} ∼ 10^{−9}. Such a small |T_{0}(ϵ_{F T})|^{2} suggests that the magnitude can be understood from a single bounce WKB-approximation for the penetration of the electron-wave through the barrier. In this approximation, the transmission probability through the barrier is given by,

|T_{0}^{WKB}(ϵ_{z})|^{2} ∼ e^{−W} ,  (2.7)

where

W = 2 ∫_{z_{1}}^{z_{2}} \sqrt{2(V (z) − ϵ_{z}})dz  (2.8)

and z_{1} and z_{2} are the classical turning points in the tip-sample interface potential V_{I}(z) of an electron with energy ϵ_{z} normal to the surface. The behaviour of W can be understood from a trapezoidal model for the barrier potential, which gives

W = \{ \begin{matrix}\frac{4\sqrt{2}}{3F}(Φ_{T}^{3/2} − (Φ_{S} − V)^{3/2}) , & V < Φ_{S} \\ \frac{4\sqrt{2}}{3F}Φ_{T}^{3/2} , & V > Φ_{S} \end{matrix}  (2.9)

where Φ_{T} is the workfunction of the tip. In the (field emission) region V > Φ_{S}, the tunneling distance, s_{T} ≡ z_{2} − z_{1}, varies with F, while in the region V < Φ_{S}, s_{T} is constant and equal to the tip-sample distance s. An estimate of the magnitude of the transmission probability based on (2.7) – (2.9) with Φ_{T} = Φ_{S} and for conditions as in Fig. 1 gives |T_{0}^{WKB}(ϵ_{F T})|^{2} ∼ 10^{−10} which is about the same magnitude as the calculated value.

The WKB-result in (2.7) and (2.8) shows also directly that |T_{0}(ϵ_{z})|^{2} will increase rapidly with ϵ_{z} due to a decreasing magnitude of the negative kinetic energy in the barrier. A Taylor expansion in ϵ around ϵ_{F T} of the action integral in (2.8) gives,

|T_{0}^{WKB}(ϵ_{z})|^{2} ∼ |T_{0}^{WKB}(ϵ_{F T})|^{2}e^{− \frac{(ϵ_{F T} −ϵz)}{Δ}} ,  (2.10)

where

Δ = \frac{1}{2 ∫_{z_{1}}^{z_{2}} \frac{dz}{\sqrt{2(V (z)−ϵ_{F T}})}} .  (2.11)

An application once again of the trapezoidal model to V_{I}(z) gives the estimate,

\frac{1}{Δ} ∼ \{ \begin{matrix}\frac{4}{\sqrt{2}F} (Φ_{T}^{1/2} − (Φ_{S} − V)^{1/2}), & V < Φ_{S} \\ \frac{4}{\sqrt{2}F} Φ_{T}^{1/2}, & V > Φ_{S} \end{matrix}  (2.12)

Under the conditions in Fig. 1, we find from the estimate in (2.12) a small value for Δ, about 0.1 eV, due to a large s_{T} ≈ 17 Å, which suggests that the tunneling is dominated by electron in states with ϵ_{z} close to ϵ_{F T}. Correspondingly, k_{∥} is close to zero. The prefactor ϵ_{F T} − ϵ_{z} in the first term of (2.5) vanishes exactly at the Fermi-level due to the vanishing density of states at ϵ_{z} = ϵ_{F T} which makes the tunneling to be dominated by states with ϵ_{z} ≈ ϵ_{F T} − Δ. The calculated (ϵ_{F T} − ϵ_{z})|T_{0}(ϵ_{z})|^{2}, shown in Fig. 3, clearly demonstrates a behaviour of rapid exponential decrease with decreasing ϵ_{z} over just a few tenths of an eV. We estimate Δ to 0.15 - 0.2 eV. Hence, the second term in (2.5) is negligible and the first term can accurately be represented by the approximation

J_{0} ≈ \frac{eΔ^{2}}{2π^{2}} |T_{0}(ϵ_{F T} − Δ)|^{2} ,  (2.13)

where Δ is now obtained from the exponential decrease of the calculated |T_{0}(ϵ_{z})|^{2}. In practice, we have already shown that Δ is so small that it can be ignored in the argument of the transmission coefficient in (2.13). Furthermore, the results from the simple trapezoidal model, Eq. (2.12), shows that Δ will not be exponentially dependent on s, like the transmission probability, and can be assumed to be constant in the calculation of J_{0}.

We demonstrate further how the effects of the sample ion-cores modify the tunneling current. In principle, we need to calculate the current carried by the Bloch wave propagating into the sample. This is done here by noting that, from particle conservation, this current is given by 1 − |r_{c}|^{2} times the current of the wave incident on the crystal. The multiple scattering of the electron wave between the crystal potential and the barrier potential (see Fig. 2) results in an amplification of the wave incident on the crystal with the factor 1/(1 − r_{b}r_{c}). Hence, the transmission coefficient in this situation is simply obtained by multiplying |T_{0}|^{2} with,

R_{tunn}(ϵ) = \frac{1 − |r_{c}|^{2}}{|1 − r_{b}r_{c}|^{2}}  (2.14)

where r_{b} is the reflection amplitude of an incoming wave from the sample on the barrier.

The final result for the tunneling current density J that we use in our calculations is obtained by collecting the results in (2.13) and (2.14) and is given by,

J ≈ \frac{Δ_{0}^{2}}{2π^{2}} R_{tunn}(ϵ_{F S} + V)|T_{0}(ϵ_{F S} + V)|^{2}  (2.15)

where Δ_{0}^{2} is assumed to be constant and Δ is neglected in the argument of the transmission coefficient. The current density J in (2.15) depends on two independent control parameters – the tip-sample distance s and the bias V – through the interface potential and the argument ϵ_{F S} + V.

Scanning tunneling experiments often record the conductance versus bias, \frac{dI}{dV} − V, in a constant-current mode, i.e., the bias is varied under the condition that the current is constant. Since we are only able to calculate current densities in our planar model, we make a comparison between the observed \frac{dI}{IdV} − V in a constant-current mode and the calculated \frac{dJ}{JdV} − V in a constant-current-density mode. The latter quantity has been calculated in the following way; (i) for initial V and J which fixes F and s, \frac{dJ}{dV} is calculated using (2.15), upon a small change of V by δV; (ii) V is increased and F and s are adjusted so that J across the new barrier is the same within 3% of the original J; (iii) \frac{dJ}{dV} is then calculated as in (i) for this new V. A smaller tolerance than 3% in the current gives no change in the resulting curves and the functional dependence of V_{I}(z) on F is obtained by an interpolation between potentials calculated for a set of different fields.

D. Surface density of states

Since we are in regime where the dominant contribution to the tunneling comes from states with ϵ_{z} close to ϵ_{F T} = ϵ_{F S} + V, we attempt to relate \frac{dJ}{JdV} to a surface DOS of the sample in the presence of a positive external uniform field. In order to do that, we need first to define the surface DOS for electrons with k_{∥} = 0 and show how it is influenced by the crystal reflectivity.

In a situation where we are concerned with the scattering of electrons that are incident normal to the crystal, the propagation is described by a one-dimensional potential, V_{S}(z). The surface density of states n_{s}(ϵ) of the surface region z > z_{0}, is then defined as an integral over the local DOS,

n_{s}(ϵ) ≡ ∫_{z_{0}}^{∞} ρ(z; ϵ)dz,  (2.16)

where the local DOS, ρ(z; ϵ), is determined from the retarded Green function, g(z, z′; ϵ), of V_{S}(z) as,

ρ(z; ϵ) = −\frac{1}{π} Im g(z, z; ϵ) .  (2.17)

Note that we have explicitly excluded the spin-degeneracy in (2.16). We evaluate ρ(z; ϵ) in (2.17) using a standard representation of the Green function in terms of solutions ψ_{<}(z) and ψ_{>}(z) to the one-electron Schrödinger equation in V_{S}(z), as

g(z, z′; ϵ) = \frac{2ψ_{<}(z_{<})ψ_{>}(z_{>})}{W}  (2.18)

where z_{<} = min(z, z′), z_{>} = max(z, z′). The Wronskian W = ψ_{<}(z)\frac{dψ_{>}(z)}{dz} − ψ_{>}(z)\frac{dψ_{<}(z)}{dz} is independent of z. The wave functions satisfy the boundary conditions that ψ_{<}(z) should give rise to an outgoing Bloch wave in the crystal (z → −∞) and ψ_{>}(z) should decay to zero on the barrier side (z → ∞). In the infinitesimal region around z_{0} , these two wavefunctions are determined from the crystal and barrier reflectivity, respectively, and can be expressed as,

\begin{matrix}ψ_{<}(z) = e^{−ik(z−z_{0})} + r_{c}e^{ik(z−z_{0})} \\ ψ_{>}(z) = e^{ik(z−z_{0})} + r_{b}e^{−ik(z−z_{0})}. \end{matrix}  (2.19)

The corresponding Wronskian W is equal to 2ik(1 − r_{b}r_{c}) and ρ(z; ϵ) can be shown after some straightforward algebra to be given by,

ρ(z; ϵ) = R(ϵ)ρ_{0}(z; ϵ)  (2.20)

where

R(ϵ) = \frac{1 − |r_{c}|^{2}}{1 + |r_{c}|^{2} − 2|r_{c}| cos(ϕ)} ,  (2.21)

and ϕ = ϕ_{b} + ϕ_{c}. Note that R(ϵ) is independent of z. Here ρ_{0}(z; ϵ) is the local DOS when r_{c} = 0 and is given by,

ρ_{0}(z; ϵ) = \frac{1}{2πk} |ψ_{>}(z)|^{2}  (2.22)

Similar expressions for ρ(z; ϵ) have also been given by Radny [37] using a semiclassical approach. In the evaluation of the integral in (2.16) over ρ(z; ϵ), we use a method introduced by Echenique and Pendry [1,9] which is based on introducing an infinitesimal absorption in the Schrödinger equation. An application of this method to the integral over |ψ_{>}(z)|^{2} gives

\frac{1}{k} ∫_{z_{0}}^{∞} dz|ψ_{>}(z)|^{2} = \frac{dϕ_{b}}{dϵ} (ϵ) + \frac{sin ϕ_{b}(ϵ)}{2ϵ}  (2.23)

The term \frac{sin ϕ_{b}}{2ϵ} comes from the interference between the incident and scattered wave from the barrier [38] and is dependent on the position of z_{0}. Using the result in (2.23) and the definition (2.16) of n_{s}(ϵ) we obtain, in the absence of a crystal reflectivity, r_{c} = 0,

n_{s}^{0}(ϵ) = \frac{1}{2π} [\frac{dϕ_{b}}{dϵ} (ϵ) + \frac{sin ϕ_{b}(ϵ)}{2ϵ}] ,  (2.24)

and in the situation of a non-zero crystal reflectivity,

n_{s}(ϵ) = R(ϵ)n_{s}^{0}(ϵ) .  (2.25)

Hence, all the effects of the crystal reflectivity both on surface and local DOS are just included in R(ϵ) as defined in (2.20). Note that we have a similar result for the effects on the transmission coefficient in (2.14) through the factor R_{tunn}(ϵ). This factor differs only from R(ϵ) by the fact that r_{b} is defined in (2.14) for V_{I}(z) while in (2.21) it is defined for V_{S}(z). The differences between r_{b} are very small in these two situations since the transmission coefficient is very low through V_{I}(z).

In the calculations of n_{s}(ϵ) and n_{s}^{0}(ϵ), the barrier phase shift, ϕ_{b}, and its energy derivative are obtained from the wavefunctions calculated numerically from the one-electron Schrödinger equation in V_{S}(z). The magnitude of the crystal reflectivity |r_{c}| and its phase shift, ϕ_{c}, have been extracted from Eq. (2.4).

III. RESULTS AND DISCUSSION

We begin by presenting tunneling spectra, \frac{dJ}{JdV} calculated in a constant-current-density mode in the tip-sample interface potential. These spectra are then related to the surface electronic structure of the sample potential in the presence of a static external field by examining the calculated surface density of states, n_{s}(ϵ) for electrons with k_{∥} = 0. In particular, we study the effects of the crystal reflectivity on the tunneling spectra and compare to surface density of states. Finally, we make a direct comparison of the calculated tunneling spectra with experimental data.

A. Calculated tunneling spectra and surface DOS

In Fig. 4, we present the effects of the crystal reflectivity on the \frac{dJ}{JdV} −V curves calculated in a constant J mode for a field strength that will be shown in the next subsection to be appropriate when making a comparison with experiments. The approximation (2.15) that we use in the calculation of \frac{dJ}{JdV} −V curves in the constant J mode corresponds to a constant |T (ϵ_{F T})|^{2}. At the lowest bias of about V =2.8 eV, we have |T (ϵ_{F T})|^{2} = 10^{−7} and s= 15 Å. In order to keep |T (ϵ_{F T})|^{2} constant at higher V, s has to increase from 15 to 42 Å while F increase (in an oscillatory manner) from 0.17 to 0.21 eV/Å when V increases from 2.8 to 6.8 eV. This behaviour is in line with what is expected from the WKB treatment of the trapezoidal model introduced in section II A. According to (2.7), |T (ϵ_{F T})|^{2} is governed by the tunneling distance s_{T} when V > Φ_{S} and s has to increase almost linearly with V in order to keep s_{T} constant. This variation of s with V also makes F stay more or less constant. We note that the largest variation of F in our calculations occur at low biases, V ≤ Φ_{S}. Finally, we note that this behaviour of s_{T} and F with V in the constant J mode is consistent with the approximation of keeping Δ constant according to the WKB-result in (2.12). Moreover the effects on J due to the broadening and shift from Δ ≈ 0.2 eV which is neglected in (2.15) is much smaller than the energy scale for the structures in the calculated \frac{dJ}{JdV} − V curve which are ∼ 1 eV.

The oscillatory structure in the calculated \frac{dJ}{JdV} − V curve in Fig. 4 for V > Φ_{S} is present independent of crystal effects of the sample. This kind of oscillatory structure has been identified previously both in previous scanning tunneling experiments [5,6] and in theoretical models similar to our model [7,21,39]. The issue of the characterization of this oscillatory structure in terms of “field-induced” resonances has been brought up earlier [5–7,21] with different viewpoints. Binnig et al [5] and Garcia et al [21], interpreted all peaks as being a Stark shifted Rydberg series, whereas Becker, Golovchenko, and Swartzentruber [6] presented their results as standing waves formed between the probe and the sample, as Gundlach resonances [40], having nothing to do with surface electronic states. The origin and behaviour of this structure will be examined and discussed again in more detail in connection with the presentation of the results of surface DOS.

The most important point of our work is the peak at V below Φ_{S} in the calculated \frac{dJ}{JdV} − V curve which is present only when the effects of the crystal reflectivity are included. This peak has a completely different origin than the oscillatory structure present for biases V larger than Φ_{S}. It is a “crystal-derived” barrier resonance formed by a reflection back and forth of the electron wave from the crystal and the barrier. The precise position in V of this resonance below Φ_{S} is dependent on the behaviour of ϕ_{c} and ϕ_{b}.

As we have already stressed, our approximate expression for the current density in (2.15) is well-justified and shows that J reveals the sample electronic structure at an energy ϵ_{F S}+V. The reflection coefficient from the barrier part of V_{I}(z) at the energy ϵ_{F S} + V will be very close to the corresponding coefficient for V_{S}(z) with the same F. Thus it is of direct interest to try to relate the \frac{dJ}{JdV} − V curves to the surface DOS, of V_{S}(z) at a fixed F.

In Fig. 5, we present the results for n_{s}^{0}(ϵ) and n_{s}(ϵ) of V_{S}(z) calculated at fixed F without and with the effects of the crystal reflectivity, respectively. The field F has been chosen to cover the narrow range of fields F (v; J) generated in the constant J mode. There is a clear one-to-one correspondence between the peaks in the calculated \frac{dJ}{JdV} − V curves in Fig. 4 and n_{s}(ϵ) in Fig. 5. The oscillatory structure due to the “field-induced” resonances, is also present in the calculated n_{s}^{0}(ϵ) and n_{s}(ϵ) at energies above ϵ_{F S} +Φ_{S}. When the lattice effects are switched on, the peaks above the vacuum level narrow and shift slightly to lower energies and a “crystal-derived” peak occur below the vacuum level due to the maximum of R(ϵ) at that energy.

The “crystal-derived” barrier resonance in n_{s}(ϵ) comes from the multiple scattering between the crystal and the surface barrier which is described by the prefactor, R(ϵ), to n_{s}^{0}(ϵ) in (2.25). We will first discuss the case where F = 0. In a situation where |r_{c}| is close to one, [41] e.g., close to a band gap, R(ϵ) gives rise to narrow Lorentzian resonance peaks, ∼ ((ϵ − ϵ_{m})^{2} + Γ^{2})^{−1}, at energies ϵ_{m} determined by a phase condition corresponding to constructive interference as,

ϕ(ϵ_{m}) ≡ ϕ_{b}(ϵ_{m}) + ϕ_{c}(ϵ_{m}) = 2πm ,  (3.1)

where m is an integer. The full width at half-maximum, 2Γ, of the Lorentzian is given by,

Γ = \frac{1 − |r_{c}|}{\sqrt{|r_{c}|}|ϕ′(ϵ_{m})|}  (3.2)

where ϕ′(ϵ_{m}) ≡ \frac{dϕ}{dϵ}(ϵ_{m}). An image barrier (without an external electric field) is known to give rise to a divergent ϕ_{b} and \frac{dϕ_{b}}{dϵ} at the vacuum level which results in the formation of a well-defined image series with a decreasing width converging to ϵ_{vac}. This behaviour of the phase shift can be understood within the WKB-approximation which is well justified for a 1/z-potential like the image potential. In this approximation, the phase shift is given by,

ϕ_{b}^{WKB}(ϵ) = 2 ∫_{z_{ref}}^{z_{t}} k(z)dz − \frac{π}{2}  (3.3)

where z_{t} is the classical turning point for an electron with energy ϵ in V_{S}(z), and k(z) = \sqrt{2(ϵ − V_{S}(z))}. Explicitly, ϕ_{b}^{WKB} for a −\frac{1}{4z} potential (z_{ref} = 0) is given by [11]

ϕ_{b}^{WKB}(ϵ) = π [\sqrt{\frac{1}{8(ϵ_{vac} − ϵ)}} − 1] .  (3.4)

Note that both ϕ_{b}^{WKB}(ϵ) and \frac{dϕ_{b}^{WKB}(ϵ)}{dϵ} are divergent at ϵ = ϵ_{vac}. From the phase condition in (3.1) and neglecting the energy variation of ϕ_{c}, a Rydberg series is obtained with binding energies (relative ϵ_{vac}) as

ϵ_{m} = −\frac{1}{32(m + a)^{2}} ,  (3.5)

where a = \frac{1}{2} (1 − \frac{ϕ_{c}}{π} ). The phase shift in the Born approximation, ϕ_{c} = π, gives that the lowest member in the image series lies at -0.85 eV below ϵ_{vac}, which is close to ∼ -1 eV that we obtain in our surface barrier model. This difference can be understood from that the image-plane position z_{IP} in our model lies at ∼ 1.5 a_{0} instead of z_{IP} = 0, which gives rise to a larger barrier phase shift ϕ_{b} than does Eq. (3.4).

At values of |r_{c}| well below 0.5 like in our case where we are far from a band gap, the resonance line-shape is no longer like a Lorentzian and R(ϵ) is well approximated by only including a single scattering event from the crystal,

R(ϵ) ≈ 1 + 2|r_{c}| cos ϕ(ϵ)  (3.6)

where the weak energy dependence of |r_{c}| has been ignored. In this situation, R(ϵ) will show an oscillatory behaviour with a series of maxima also at ϵ_{m}.

When turning on the electric field F, V_{S}(z) increases in the surface region which in turn should, according to the result in (3.3) give an overall decrease and convergence of ϕ_{b} to a finite value at ϵ_{vac} since z_{t} is finite for all ϵ. This effect of ϕ_{b} by a non-zero F will limit the number of members in the image-series below ϵ_{F S} +Φ_{S}, i.e. the vacuum level for F = 0, and shift up their energies with increasing F. This behaviour is demonstrated by the results for peak positions versus field strength in Fig. 6 where there exists only one, by us interpreted as a “crystal-derived”, barrier resonance even down to F =0.05 eV/Å whose peak position also increases in energy with F. This result raises immediately the question: what happens to the members of the image series which “disappears” at finite F? The answer is simply that they shift up above ϵ_{F S} +Φ_{S} and change their character into “field-induced” resonances since their energy position will now be determined by the contribution to ϕ_{b} from the linear part of the potential.

We turn to an examination of the nature of the “field-induced” resonances. In this case it is not necessary for their existence to have a finite reflectivity from the crystal. This can be understood from the following. The sharp variation in the surface barrier potential when going from bulk to the linear region of the V_{S}(z), is sufficiently rapid, compared to the wavelength of the electron, to give rise to a non-zero reflectivity of an electron wave incident on this part of barrier from the linear region of V_{S}(z). Such a behaviour cannot be understood in the WKB-approximation as stressed, e.g., by Gundlach [40].

We can make a closer analysis of the behaviour of the “field-induced” resonances by finding a similar R̃(ϵ) to R(ϵ) that describes these resonances. Such a factor can be found from Eq. (2.21) simply by changing the location of z_{0} to a position z̃_{0} in the linear region of V_{S}(z) whereby r_{c} and r_{b} gets modified to r̃_{c} and r̃_{b}. The reflectivity, r̃_{c}, includes now also the reflectivity of the electron waves from the rapid variation of V_{S}(z) around the jellium edge of the sample, earlier included in r_{b}. The result in (2.20) is also valid in this situation and the local DOS can be expressed as,

ρ(z; ϵ) = R̃(ϵ)ρ̃_{0}(z; ϵ) ,  (3.7)

where ρ̃_{0}(z; ϵ) is determined as in (2.22) from a wave function ψ̃_{>}(z) ∝ ψ_{>}(z) but with an incident amplitude of unity at z = z̃_{0} on the barrier. All the “field-induced” resonances in n_{s}(ϵ) will now be described by R̃(ϵ). The calculated n_{s}^{0}(ϵ) in Fig. 5 shows that we do not have any narrow resonances whose width is much less than the level separation. Hence, also in this case we have a small value for |r̃_{c}| due to the reflection of V_{S}(z) at the selvedge and the lineshape is given by an oscillatory function like in (3.6). A small value for |r̃_{c}| in this situation is consistent with the observation that the inclusion of the crystal reflectivity makes the amplitudes of the “field-induced” resonances to increase with more than a factor of two. The non-linear variation of the “field-induced” resonance energies with F as shown in Fig. (6) can be understood from the behaviour of ϕ̃_{b} with respect to F. We evaluate the WKB-approximation in (3.3) for ϕ̃_{b} in a linear repulsive potential V_{S}(z̃_{0}) + F (z − z̃_{0}). The resulting approximate phase shift is then given by,

ϕ̃_{b}(ϵ) ≈ \frac{4\sqrt{2}}{3} \frac{ϵ^{3/2}}{F} − \frac{π}{2} .  (3.8)

If the energy variation of ϕ̃_{c} is neglected, we obtain directly from (3.1) the approximate resonance energies (relative V_{S}(z̃_{0})),

ϵ_{m} ≈ (\frac{3π}{\sqrt{8}})^{2/3} F^{2/3}(m + \frac{1}{4} + α)^{2/3}.  (3.9)

where α is determined from ϕ̃_{c}. For instance, in the case of a triangular well with an infinite potential barrier at z = z̃_{0}, α is equal to −\frac{1}{2}. The non-linear power-law behaviour of ϵ_{m} with F in (3.9) with an exponent less than one fits the calculated variation of the “field-induced” resonances in Fig. 6. Best agreement is obtained for the high-lying resonances.

B. Comparison with experiments

The prime objective with the comparison of the results from our one-dimensional potential model with the measured tunneling spectra on Al(111) is to show that one of the observed peaks in the measured spectra originates from the “crystal-derived” barrier resonance of the sample potential. Since the work function of the Ir-tip is unknown due to its likely covering by Al atoms, we will discuss the experimental data in relation to results both from an Ir- and an Al-tip model.

The experimental tunneling spectra are recorded in a constant current mode under different initial tunneling-gap distances, corresponding to several values of the tunneling current. The relative change in distance with bias is known experimentally, but the initial gap distance cannot be determined. Spectra for current values in the range 4 - 102 pA are shown in Fig.7. Two peaks are clearly discernible in all spectra and at the lowest tunneling current even a third peak can be resolved at higher biases. An increase of the bias from 4 to 7 eV at I ∼ 100 pA corresponds to an increase of about 14 Å [26]. This observation in the field emission region of the tunneling suggests that the field strength is about 0.2 eV/Å. In Fig. 8, we have consequently presented results of \frac{dJ}{JdV} − V in the Ir-tip model for a set of initial field strengths: 0.20, 0.25, and 0.30 eV/Å and for the Al-tip model for a field strength of 0.20 eV/Å. Note that there are no differences in the peak positions for the two models.

The calculated and measured tunneling spectra agree in important ways. The observed peak separations in the measured \frac{dI}{IdV} − V are well reproduced by the narrow range of field strengths of 0.2 to 0.25 eV/Å in the calculated \frac{dJ}{JdV} −V. The calculated data reproduces also nicely the observed shift of the peak positions to higher V with increasing I, corresponding to higher J in the calculations.

In our earlier discussion in the previous section, we have demonstrated that the first peak in the calculated tunneling spectra with a position below Φ_{S} is a “crystal-derived” barrier resonance and that the higher lying peaks derive from “field-induced” resonances of the sample potential. The fact that there is an observed peak below the vacuum level leads us to the interpretation that it is a “crystal-derived” barrier resonance while the other peaks, well above the vacuum level, are due to “field-induced” resonances. The observed position of the “crystal-derived” peak in Fig. 7 seems to occur at a bias ∼ 0.5 V higher than the calculated peak position. In the calculations of the tunneling spectra, we have neglected that the effective tunneling occurs Δ ∼ 0.2 eV below ϵ_{F T} (see Eq. (2.11)). Accounting for that, would shift the calculated peak position up in bias with Δ, yielding better agreement with experiment. However, we would like to stress that the peak position is sensitive to the model potential used, which suggests that further work on the improvement of the model potential both in the surface region and for the crystal would be of great interest. Note that the value of Φ_{S} in the jellium model is about 0.4 eV lower than the measured Φ of Al(111). We compensate for that by making the comparison of the peak positions with respect to V − Φ_{S} instead of V in Fig. 7. Why the peak positions shift to higher V with increasing I is understood in our model from that an increase of J at a fixed V corresponds to a decrease in the tunneling distance and an increase of F. Such an increase of F will increase the resonance positions in n_{s}(ϵ) of V_{S}(z) as shown in Fig. 6 with a concomitant increase of the positions of the corresponding peaks in the tunneling spectra.

The positions of the peaks in the calculated tunneling spectra are insensitive to our modelling of the work function of the tip as can be seen from the comparison of the calculated tunneling spectra in Fig. 8 for the Ir- and Al-tip where Φ_{T} − Φ_{S}=1.5 eV and Φ_{T} = Φ_{S}, respectively. This effect can be understood from that the resonance positions are determined by the interface potential on the sample side which is not very much influenced by the presence of the tip in our situation of small transmission probabilities. The prime effect of the Ir-tip is that there will be a contact potential difference of about 1.5 eV between the tip and the sample even at zero bias. This difference makes the tunneling distances larger for the Ir-tip at the same F and V than for the Al-tip with no contact potential difference.

In this context, we would like to comment on the relevance of dynamical effects of the tunneling through the potential barrier on the peak positions in the tunneling spectra. There have been some model studies on the dynamical effects of the image potential on the tunneling through a potential barrier relevant for STM [42]. These studies have shown that there is a reduction of the contribution of the image potential to the interface potential if the “traversal” time for tunneling τ is smaller than the inverse surface plasmon frequency ω_{sp}^{−1}, i.e, if τω_{sp} < 1. If this regime would apply in our case, a static image potential approximation would then give a too low potential barrier with the implication that the calculated peak positions would appear at a lower bias in the tunneling spectra than in a calculation based on a dynamic image potential. However, we find that for the relevant tunneling conditions that a static approximation of the image barrier potential is well justified in our model. This is based on the WKB result for the “traversal” time [43], τ = \frac{1}{2Δ} and that Δ is ∼ 0.2 eV which makes ω_{sp}τ = \frac{ω_{sp}}{2Δ} ≈ 25 for Al where ω_{sp} ∼ 10 eV. Dynamic effects are thus unlikely to have a major influence on the image potential here.

IV. SUMMARY

We have extended a previously developed model [14] for formation of electronic resonances on free-electron-like surfaces to include static external fields, thus enabling us to compare with scanning tunneling spectroscopy (STS) data from such surfaces [22]. The model makes use of a self-consistent jellium potential under the influence of uniform external fields which mimic the presence of a biased tip. The effect of the weak Al pseudopotential due to the ion-core lattice is included perturbatively. This model is appropriate for application on free-electron-like surfaces at energies sufficiently far from band gaps. The crystal reflectivities obtained in our approach and from the often used two-band model are compared and it is found that the two-band model underestimates the reflectivity of the crystal at energies far from the considered band gap.

We demonstrate that a peak in the calculated conductance curves, below the vacuum level, is a “crystal-derived” resonance since it is not present in the calculated conductance curves when the lattice effects are omitted. The observed peak in the STS spectra below the vacuum level is therefore attributed to this resonance. Moreover, we argue that the precise position of this peak is a test of the surface barrier and the crystal potential. Comparison of the calculated conductance and Al surface density of states (SDOS) shows a one-to-one correspondence between the resonances in the (SDOS) and peaks in the conductance curves. The resonances above the metal vacuum level are due to multiple scattering between the linear part of the barrier and a sharp part of the potential when going from the linear region to the bulk. Extrapolation to zero field demonstrates that the resonance in the SDOS below the Al vacuum level for finite fields is mapped onto a “crystal-derived” image resonance at zero field. We believe that such resonances influence the IPE spectra and thus that IPE features can be interpreted in such terms and not only in terms of matrix element effects as suggested earlier [16].

ACKNOWLEDGMENTS

This work has been supported by the Swedish Natural Science Research Council (NFR). Discusssions with Bob Bartynski are gratefully acknowledged.

APPENDIX: DERIVATION OF THE REFLECTIVITY OF A SEMI-INFINITE CRYSTAL

1. Born approximation

The scattering from a semi-infinite crystal of a nearly free-electron metal may be described by a weak potential V_{c}, and a treatment based on the first order Born approximation which is justified at energies far away from any band-gap. We derive the result for the reflectivity from a semi-infinite crystal with the crystal side on z < 0 using a local V_{c}, but the derivation can also be done for a non-local and energy-dependent pseudopotential.

The first order Born approximation is based on the approximation of replacing the total wave in the integrand of the Lippman-Schwinger equation by the incoming wave as,

Ψ^{Born}(r) = Ψ_{in}(r) + \frac{2m}{h̄^{2}} ∫ d^{3}r′G_{0}(r − r′)V_{c}(r′)Ψ_{in}(r′) ,  (A1)

where Ψ^{Born}(r) is the total wave, Ψ_{in}(r) = e^{ik_{i}·r}, is the incident wave with k_{iz} < 0. The two-dimensional scattering geometry makes it convenient to use a wave-vector expansion parallel to the surface of the retarded free-particle Green function given by,

G_{0}(r − r′) = −\frac{i}{2} ∫ \frac{d^{2}k_{∥}}{(2π)^{2}} \frac{e^{ik_{z}|z−z′|} e^{ik_{∥}·(r_{∥}−r′_{∥})}}{k_{z}} ,  (A2)

where k_{z} and k_{∥} denote the momentum components perpendicular and parallel to the surface, respectively. For the propagating waves, we have k_{∥} < k_{i} and k_{z} = \sqrt{k_{i}^{2} − k_{∥}^{2}} while for the evanescent waves k_{∥} > k_{i} and k_{z} = i\sqrt{k_{∥}^{2} − k_{i}^{2}} .

The crystal potential, V_{c}(r), in (A1) is now expressed as a sum over atomic potentials, v(r), located at atomic sites r_{n} as,

V_{c}(r) = ∑_{n} v(r − r_{n}) .  (A3)

Inserting this sum into (A1) and using the wave-vector representation in (A2) gives the scattered wave outside the scattering potential, z > 0, as,

Ψ_{scatt}(r) = −\frac{im}{h̄^{2}} ∫ \frac{d^{2}k_{∥}}{(2π)^{2}} \frac{e^{i(k_{z}z+k_{∥}·r_{∥})}}{k_{z}} ∑_{n} ∫ d^{3}r′v(r′ − r_{n})e^{i(k_{i}−k)·r′} ,  (A4)

The integral over r′ in (A4) is a Fourier transform of the atomic potentials and the sum over n factorizes into an atomic form factor and a structure factor as,

∑_{n} ∫ d^{3}r′v(r′ − r_{n})e^{i(k_{i}−k)·r′} = v(k − k_{i}) ∑_{n} e^{i(k_{i}−k)·r_{n}} .  (A5)

Next, we make use of the translational symmetry of the semi-infinite system of atoms with one atom in each unit cell as,

r_{n} = m_{1}a + m_{2}b + m_{3}c ,  (A6)

where a and b are surface base vectors for atoms within an atomic layer parallel to the surface, c is a vector connecting an atom to another neighbouring atom in the next atomic layer further into the crystal, and m_{1}, m_{2} and m_{3} are integers with the restriction that m_{3} is a positive integer. The structure factor in (A5) can now be expressed as,

∑_{n} e^{i(k_{i}−k)·r_{n}} = ∑_{m_{1},m_{2}=−∞}^{∞} e^{i(k_{i∥}−k_{∥})·(m_{1}a+m_{2}b)} ∑_{m_{3}=0}^{∞} e^{i(k_{i}−k)·(m_{3}c)}  (A7)
= \frac{4π^{2}}{|a × b|} ∑_{g_{∥}} δ((k_{i∥} − k_{∥}) − g_{∥})\frac{1}{1 − e^{i(k_{i}−k)·c}} .  (A8)

In the last identity, we have performed the summation over the atomic layers and the atoms within each layer, and we have that g_{∥} is a reciprocal surface lattice vector.

The delta functions in (A8) makes the integration over k_{∥} in (A4) straightforward and reduces Ψ_{scatt} to a sum over outgoing plane waves with wavevectors k_{f} satisfying the twodimensional Bragg condition that k_{f∥} − k_{i∥} = g_{∥} and are given by,

Ψ_{scatt}(r) = ∑_{g_{∥}} f(k_{f} = k_{i∥} + g_{∥} + \sqrt{k_{i}^{2} − (k_{i∥} + g_{∥})^{2}}ẑ, k_{i})e^{ik_{f} ·r}  (A9)

where the scattering amplitude f(k_{f}, k_{i}) (or r_{c}) is given by,

f(k_{f}, k_{i}) = −\frac{im}{h̄^{2}} \frac{1}{|a × b|} \frac{v(k_{f} − k_{i})}{|k_{iz}|(1 − e^{i(k_{i}−k_{f})·c})} ,  (A10)

Note that there is a sign-error in the exponent in Eq. (3), corresponding to the expression in (A10), of our previous work [14]. In the more general situation for non-local potentials, where we still can express the crystal potential as a lattice sum over the atomic potentials like in (A3), the derivation still goes through with v(k_{f} − k_{i}) replaced by < k_{f}|v|k_{i} > in (A10).

2. Connection between the two-band model and the Born approximation

We will here show how the reflectivity obtained in a two-band model as (e.g.) applied by Radny to Al [15] is related to the result in the Born approximation. We will consider a situation most appropriate for the tunneling experiment, normal incidence and specular scattering, corresponding to the restriction k_{∥} = 0 and k′_{∥} = 0 for the incident and scattered wave, respectively.

In a two-band model the Schrödinger equation for an electron in an infinite crystal is solved in the basis of two plane waves, e^{−ik_{B}z} and e^{−i(k_{B}−g)z},

Ψ_{B}(z) = ae^{−ik_{B}z} + be^{−i(k_{B}−g)z}  (A11)

where g is a reciprocal lattice vector normal to the surface. We will be primarily concerned with the case, of g = 2π/d where d is the perpendicular distance between the atomic planes. The amplitudes a and b and the electron energy ϵ are determined from,

( \begin{matrix}\frac{1}{2}k_{B}^{2} − ϵ & V_{g} \\ V_{g} & \frac{1}{2}(k_{B} − g)^{2} − ϵ \end{matrix} ) ( \begin{matrix}a \\ b \end{matrix} ) = ( \begin{matrix}0 \\ 0 \end{matrix} ) .  (A12)

where V_{g} is the Fourier component of the crystal potential corresponding to the wave vector g. The solution to the eigenvalue problem in (A12) is determined from the amplitude ratio,

\frac{b}{a} = \frac{2V_{g}}{2ϵ − (k_{B} − g)^{2}}  (A13)

where ϵ is determined from the secular equation to (A12) as,

ϵ = \frac{k_{B}^{2}}{2} + \frac{2 | V_{g} |^{2}}{2ϵ − (k_{B} − g)^{2}}  (A14)

In particular, at the zone-boundary, k_{B} = g/2, a band gap opens up with a width 2|V_{g}| at the energy ϵ_{g} ≡ (g/2)^{2}/2.

The reflectivity, r_{c}, of the semi-infinite crystal is now obtained by matching, at the surface (z = 0), an incident and reflected plane wave in free space, Ψ(z) = e^{−ikz} + r_{c}e^{ikz}, z > 0, to the two-band solution Ψ_{B}(z). In this work we are only concerned with the case where the energy is far away from the bandgap, |ϵ − ϵ_{g}| ≫ 2|V_{g}| where k_{B} ≃ k. In this limit, r_{c} is small and proportional to the amplitude of the wave e^{−i(k_{B}−g)z} in the crystal (r_{c} ≃ (b/a)(g/(2k))),

r_{c} ≃ \frac{1}{V} \frac{v(g)}{k(2k − g)} ,  (A15)

where V_{g} has been expressed in terms of v(g) and the volume, V, of the primitive cell, V =| a × b | d. Using the same value for V_{g} ≈ 0.24 eV as Radny [15] for g = 2π(1, 1, 1)/a of the Al(111) surface, we obtain a reflectivity amplitude of only about 1.1 % at an energy close to the vacuum level, to be compared to the about 10 times larger amplitude obtained in the Born approximation. Note that the pseudo-potential form factor [33], as used in the Born-approximation, reproduces correctly this value for V_{g} which agrees with the measured band gap.

In order to understand the origin of the different magnitudes of r_{c}, we make a comparison of the reflectivity obtained in the two-band model with the result in (A10) from the Born approximation . According to (A10), the Born approximation for the reflectivity is given by

r_{c}^{Born} = −i\frac{1}{| a × b |} \frac{v(2k)}{k} \frac{1}{1 − e^{i2kd}} .  (A16)

This result reduces to the result in the two-band model if the following two conditions are fulfilled: i) |k − g/2| ≪ g/2 and ii) v(2k) ≈ v(g). Condition i) allows a Taylor expansion of the exponential function in Eq. (A16) around k = g/2 which simplifies r_{c} to

r_{c}^{Born} ≃ \frac{1}{V} \frac{v(2k)}{k(2k − g)} .  (A17)

The only difference between the results in the Born approximation and the two-band model, Eqs. (A17) and (A15), respectively, is that in the former case the atomic form factor is probed at the actual momentum transfer, while in the latter case the argument for the form factor is fixed at a momentum transfer corresponding to the reciprocal lattice vector g. Hence the two results agree if the condition ii) is also fulfilled.

The prime reason for the breakdown of the two-band model for describing r_{c} close to the vacuum level of Al(111) is that the condition ii) is violated. v(2k) is about a factor of 6 larger than v(g) at the relevant range of k [33], (k ≈ 3g/4), which accounts for most of the discrepancy between the reults for r_{c} in the two-band model and the Born approximation. Moreover, the restriction of the two-band model to include only two plane waves is simply not justified to describe weak scattering far away from band gaps. This can be understood from the fact that when including additional plane waves from higher order reciprocal lattice vectors in Eq. (A11), the corresponding contributions to r_{c} will be of the form given in Eq. (A15) for energies far way from band gaps. However, these additional contributions can be sizable since each one decays slowly away from the band gap introduced by the corresponding g-vector. The primary aim of the two-band model is rather to describe the opposite situation of strong scattering close to a band-gap.

REFERENCES

† Present address: Institute of Physics, 142 Riia Street, EE-2400 Tartu, Estonia.

∗ Present address: DSM/DRECAM/SRSIM, Centre d’Etudes de Saclay, F-91191 Gif sur Yvette, France. § To whom correspondence should be adressed. Electronic adress: tfymp@fy.chalmers.se [1] P. M. Echenique and J. B. Pendry, J. Phys. C 11, 2065 (1978).

[2] E. G. McRae, Rev. Mod. Phys. 51, 541 (1979).

[3] R. O. Jones, P. J. Jennings, and O. Jepsen, Phys. Rev. B 29, 6474 (1984).

[4] P. Nordlander and J. C. Tully, Surf. Sci. 211/212, 207 (1989).

[5] G. Binnig, K. H. Frank, H. Fuchs, N. Garcia, B. Reihl, H. Rohrer, F. Salvan, and A. R. Williams, Phys. Rev. Lett. 55, 991 (1985).

[6] R. S. Becker, J. A. Golovchenko, and B. S. Swartzentruber, Phys. Rev. Lett. 55, 987 (1985).

[7] J. M. Pitarke, F. Flores, and P. M. Echenique, Surf. Sci. 234, 1 (1990).

[8] R. Garcia, Phys. Rev. B 42, 5479 (1990).

[9] P. M. Echenique and J. B. Pendry, Prog. Surf. Sci. 32, 111 (1990).

[10] M. Weinert, S. L. Hulbert, and P. D. Johnson, Phys. Rev. Lett. 55, 2055 (1985); F. J. Himpsel, Comm. Cond. Mat. Phys. 12, 199 (1986); J. B. Pendry, C. G. Larsson, and P. M. Echenique, Surf. Sci.166, 57 (1986); G. D. Kubiak, Surf. Sci. 201, L475 (1988).

[11] N. V. Smith, Phys. Rev. B 32, 3549 (1985).

[12] D. Heskett, K.-H. Frank, E. E. Koch, and H.-J. Freund, Phys. Rev. B 36, 1276 (1987); D. Heskett, K.-H. Frank, K. Horn, E.E. Koch, H.-J. Freund, A. Baddorf, K.-D. Tsuei, E.W. Plummer, Phys. Rev. B 37, 10387 (1988).

[13] S. A. Lindgren and L. Walldén, Phys. Rev. B 40, 11546 (1989).

[14] S. Papadia, M. Persson, and L.-A. Salmi, Phys. Rev. B 41, 10237 (1990).

[15] M. Radny, Surf. Sci. 231, 43 (1990).

[16] W. L. Schaich and J. T. Lee, Phys. Rev. B 44, 5973 (1991); W. L. Schaich, Phys. Rev. B 45, 3744 (1992).

[17] A. Baratoff, Physica B 127, 143 (1984).

[18] N. D. Lang, Phys. Rev. B 34, 5947 (1986).

[19] A. Selloni, P. Carnevali, E. Tosatti, and C. D. Chen, Phys. Rev. B 31, 2602 (1985).

[20] R. M. Feenstra, J. A. Stroscio, and A. P. Fein, Surf. Sci 181, 295 (1987).

[21] R. García, J. J. Sáenz, J. M. Soler, and N. García, Surf. SCi. 181, 69 (1987).

[22] S. Yang, R. A. Bartynski, G. P. Kochanski, S. Papadia, T. Fondén, and M. Persson, Phys. Rev. Lett. 70, 849 (1992).

[23] N. D. Lang, Phys. Rev. B 45, 13599 (1992).

[24] W. Kohn, and L. J. Sham, Phys. Rev. 140, A1133 (1965).

[25] F. Schreier and F. Rebentrost, J. Phys. C: Solid State Phys. 20, 2609 (1987).

[26] G. P. Kochanski, private communication.

[27] Note that the potential used for the tip always corresponds to a semi-infinite jellium with r_{s} = 2.07a_{0}.

[28] P. A. Serena, J. M. Soler, and N. Garcia, Phys. Rev. B 34, 6767 (1986).

[29] N. D. Lang and W. Kohn, Phys. Rev. B 7, 354 (1973).

[30] A. G. Eguiluz and W. Hanke, Phys. Rev. B 39, 10433 (1989); A. G. Eguiluz, M. Heinrichsmeier, A. Fleszar, and W. Hanke, Phys. Rev. Lett. 68, 1359 (1992).

[31] In the relevant situations, we are concerned with normal incidence and scattering. In this case, c is replaced with its normal component – the distance between layers d.

[32] The measured surface relaxation of Al(111) in [H. B. Nielsen and D. L. Adams, J. Phys. C 15, 615 (1982)] has been found to be only about 1% between the first and second layer and can be neglected.

[33] D. C. Wallace, Phys. Rev. 187, 991 (1969).

[34] H. .J. Levinson, F. Greuter, and E. W. Plummer, Phys. Rev. B 27, 727 (1983).

[35] C. B. Duke, Tunneling in Solids, Solid State Physics, Suppl. 10 (Academic Press, New York, 1969).

[36] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes - The Art of Scientific Computing, (Cambridge University Press, Cambridge 1986).

[37] M. Radny, J. Phys. C 3, 5525 (1991).

[38] It is important to also put k → k+iδ in the flux normalization of Echenique and Pendry , which seems to be omitted in their analysis, otherwise, the second term in Eq. (23) will not be obtained.

[39] J. Bono and R. H. Good, Surf. Sci. 188, 153 (1987).

[40] K. H. Gundlach, Solid-State Electronics 9, 949 (1966).

[41] Note that r_{c} can be calculated in many ways. In the Born approximation, it is not meaningful to assign values for |r_{c}| close to one, whereas it is fully justifiable in other approaches.

[42] M. T. Michalewicz and J. Mahanty, Phys. Lett. A 116, 392 (1986); B. N. J. Persson and A. Baratoff, Phys. Rev. B 38, 9619 (1988).

[43] M. Büttiker and R. Landauer, Phys. Rev. Lett. 49, 1739 (1982).

FIGURES

FIG. 1. Model of sample-tip potential, where the sample is to the left and the tip is to the right. Both tip and sample are represented by a semi-infinite jellium with r_{s} = 2.07 corresponding to aluminum. An external field F is applied (+F on the sample and −F on the tip), mimicking the effect of the applied bias voltage. The difference in Fermi energies of tip and substrate E_{F T} − E_{F S} is equal to the bias voltage V. The two self-consistent solutions for a semi-infinite jellium in the presence of ±F are matched at their respective linear part which here is around 15 a.u. from the jellium edge of the substrate (z = 0). The substrate work function Φ_{S} (defined for F = 0) and the tip band bottom U_{T} are also indicated.

FIG. 2. Schematic multiple scattering geometry. An incident wave Ψ with unit flux towards the crystal is partly reflected at a reference plane. Its reflected part is given by r_{c}Ψ where r_{c} is the (complex) crystal reflectivity. This part is reflected at the barrier reference plane which has a reflectivity r_{b}. The multiple scattered parts are summed to give the total wave Ψ_{tot}.

FIG. 3. Energy variation of the calculated transmission coefficient for an Al “tip”. The bottom of the band is at -11.7 eV.

FIG. 4. Tunneling curves, \frac{dJ}{JdV} – V, for a system of an Ir tip and an Al(111) substrate for an initial field strength of 0.20 eV/Å, with and without inclusion of the crystal pseudopotential. It is only when the crystal reflectivity is taken into account that we get a conductance peak below the metal vacuum level (3.87 eV). The lines connecting the calculated points are cubic spline interpolations.

FIG. 5. Surface density of states of the Al(111) sample in the presence of an external field F = 0.2 eV/Å. The field strength, F is kept constant for all energies. In n_{s}(ϵ), the reflectivity from the sample ion-cores is included, |r_{c}| = 0.1, while in n_{s}^{0}(ϵ) we have |r_{c}| = 0.

FIG. 6. Energies of the resonances in the surface density of states of Al(111) as a function of the applied (static) field strength. The energy position of the “crystal-derived” (image) resonance for the field-free case is also indicated.

FIG. 7. Experimental differential conductance curves for tunneling currents of 4 - 100 pA, recorded in a constant-current mode.

FIG. 8. Calculated conductance (\frac{1}{J} \frac{dJ}{dV}) curves for different initial field strengths, corresponding to different tunnel currents for an Ir tip and an Al tip. The indicated field strengths F apply only to the first point in each curve; the requirement of constant tunnel current throughout the calculations makes the field strength vary slightly with bias V. The dashed lines between the points in the curves are spline interpolations and serve primarily as a guide to the eye.