Three-dimensional approach to scanning tunneling spectroscopy and application to Shockley states

The problem of the interpretation of scanning tunneling spectroscopy (STS) data is analytically solved using a three-dimensional (3D) transfer Hamiltonian approach. We present an analytical model capable of including both the electronic structure of the sample and the symmetry of the tip states (s, pz, dz2, …) and we discuss the role of these 3D aspects in tunneling. We applied this model to the case of Shockley states. This system, allowing a full analytical treatment, led us to a detailed simulation and comprehension of the tunneling process. A procedure for the recovery of the sample local density of states from STS measurements is then proposed and applied to both the simulated and the experimental STS data of Shockley states. Comparing this approach with other methods proposed in the literature, the importance of considering the 3D aspects in treating and interpreting STS data is demonstrated.


Introduction
The investigation of the local, surface-projected electron density of states (DOS) at the nanometric scale and the possibility of accessing both the occupied and the unoccupied states are two of the most attractive features of scanning tunneling spectroscopy (STS) [1]. Over the last few years, different approaches have been followed to solve the challenging problem of interpreting the raw STS quantities and connecting them to the physical properties of surfaces.
The tunneling current is the result of a non-equilibrium process taking place on a surface in the presence of a tip (with a given electronic structure), and its general relationship with the surface DOS is indeed quite complicated. In principle, making use of the theory of the non-equilibrium Green's function [2], a full computational approach can be developed to simulate the process [3]. Various attempts for more simplified treatments exist. One of the most successful approaches consists in describing the tunneling in the one-dimensional (1D) Wentzel-Kramer-Brillouin (WKB) approximation. From the earliest works [4], it has been clear that both the tunneling current I and, in a more direct way, the differential conductivity dI /dV provide information on the electron DOS of the surface. However, such quantities generally contain other spurious contributions due both to the presence of the tip and to the tunneling process itself. As a consequence, the identification of dI /dV as a function of the applied bias V with the surface DOS is valid only in a short-energy range around the Fermi level. Stroscio et al [5] proposed to normalize the differential conductivity to the total conductivity (dI /dV )/(I /V ) to qualitatively remove the tunneling process contribution. Ukraintsev [6] proposed a more consistent normalization procedure based on the analytical evaluation of the 1D-WKB transmission coefficient of the tunneling barrier. A more refined treatment was provided by Koslowski et al [7], who proposed an analytical normalization procedure based on a suitable combination of dI /dV and I , to quantitatively recover the surface DOS. The possibility of recovering the sample DOS by solving iteratively the corresponding Volterra equation of the second kind in the Neumann scheme was also suggested [7]. Wagner et al [8] directly approached the same problem numerically inverting the integral relation between the sample DOS and I (which is a Volterra equation of the first kind) and applied it to the case of organic molecule spectra, neglecting possible tip effects. Passoni et al [9] evaluated the effectiveness of these methods on experimental data and considered the effect of the presence of a non-constant tip DOS in dI /dV spectra. The possibility of removing such tip effects has also been demonstrated in principle, exploiting iterative solution methods [10,11]. Finally, Ziegler et al [12] proposed a method to treat and interpret the spectra taken by varying the tip-sample distance in the so-called constant current mode. One of the most attractive features of the 1D-WKB approach is that it 'bypasses' the non-equilibrium situation, connecting in a simple and direct way measured quantities with physical properties. On the other hand, several important aspects cannot be treated within the 1D-WKB approximation. In particular, 3D effects have been included in the WKB framework considering the problem of tunneling through a 3D extended tunneling barrier only within the zero bias approximation [13], while the contributions of a non-vanishing applied bias, of electron states located in different regions of the sample Brillouin zone and of the symmetry of the tip states cannot be easily included. As a consequence, a 3D-WKB theory suitable for STS is not available yet.
In a different framework, the theory of tunneling can be treated within the so-called transfer Hamiltonian (TH) approach, as was originally done by Bardeen [14]. This approach, within its limits of validity [15], opens the possibility for an analytical treatment of the 3D effects.

3
From a historical point of view, these aspects were included in scanning tunneling microscopy (STM) descriptions first by Tersoff and Hamann [16] and then by Chen [17][18][19], in the low bias approximation, with the purpose of explaining the origin of atomic resolution in STM. The extension of this model to finite bias has been recently proposed with the aim of providing a solid theoretical framework for the interpretation of STS measurements [20].
We followed the possibilities offered by this approach to present a complete 3D analytical model of tunneling and a specific STS normalization procedure, applied in particular to the physical system consisting of Shockley states with a parabolic dispersion in the presence of a tip state with well-defined angular symmetry. In section 2, we present the analytical description of the tunneling process in a 3D-TH approach, providing an expression for I and dI /dV that is able to take into account both the band structure of the sample surface states and the tip state angular symmetry. In section 3, we consider the case of an ideal Shockley surface state as a sample in the STS measurement, showing 3D-related features in simulated tunneling current. The presented analytical formulation is exploited in section 4, where we propose a normalization procedure for the recovery of the sample local density of states (LDOS). The normalization method is applied to the simulated Shockley state spectra and compared with existing procedures [5][6][7]. The reliability of the method is also evaluated by testing sensitivity to the model parameters (e.g. estimated tip-sample distance and tip symmetry) and eventually applying the method to experimental STS spectra taken on Cu(111). The concluding remarks are presented in section 5.

Theory
We start by summarizing the formulation of the tunneling problem in the TH approximation. Under the quite general assumptions of independent electrons for both the sample and the tip and of weak coupling between the two physical systems, we can define the LDOS ρ S(T) (E) = γ δ(E − E γ ) and the LDOS L S(T) (E, r) = γ |ψ γ (r)| 2 δ(E − E γ ), where E γ and ψ γ are the γ th eigenvalue and eigenfunction of the unperturbed single-electron Hamiltonian, respectively. In the following, we will use the notation γ → µ for the tip states and γ → ν for the sample states. In the 1D treatment, a difference between LDOS and DOS does not explicitly appear and these two quantities are used almost equivalently. On the other hand, they have to be distinguished whenever we consider a real 3D system. One of the most attractive features of the 1D-WKB approximation is that it leads to a formulation of the tunneling current in the form of a convolution integral of tip and sample DOSs, weighted by the WKB tunneling transmission coefficient [21]. Such an expression is achieved after several simplifying assumptions and it is not valid in the general case, as we will show in the following. The TH formulation for the tunneling current in the zero-temperature, elastic regime [14,22], using dimensionless units where V is the applied voltage, W F is the effective work function and e, m e are the modulus of electron charge and mass, respectively), is given by Here, the probability of a transition between the sample state with the energy E ν = E and the tip state with the energy E µ = E − V is described by means of the Fermi golden rule through the matrix elements M ν,µ , which depend on the corresponding wave functions ψ µ , ψ ν [14].
In the general situation, the evaluation of M ν,µ in (1) does not lead to a simple connection between the current and the sample DOS or LDOS [15,20]. If we expand the tip wave functions into spherical harmonics [18], it is possible to write these matrix elements in the general form [20] where r 0 is the center of curvature of the tip apex andÔ µ (E, V ) is a differential operator that acts on the spatial coordinates of the sample wave functions ψ ν and depends on both the tip properties and the applied finite bias. If we assume that (a) the sample wave functions are the solutions of a Schrödinger equation that is separable into Cartesian coordinates (x, y for the in-plane parallel ( ) directions, which are associated with an eigenvalue E , and z for the out-of-plane orthogonal (⊥) direction, which is associated with an eigenvalue E z = E − E ). (b) The potential barrier between the two systems can be well approximated by a (biasdependent) effective rectangular barrier of height equal to W F + V /2; we can express the sample wave function in the vacuum region in the form where A ν (x, y) and χ ν ⊥ (E z , V, z) are the eigenfunctions for the parallel and perpendicular equations which depend on the quantum numbers ν ⊥ and ν , respectively; k S = √ 1 + V /2 − E z is the decay factor for the sample state in vacuum and C ν ⊥ (E z , V ) is a normalization parameter. Although for some special cases such a parameter could be calculated, its dependence on V is in general not known and some assumptions are therefore needed. Here, we neglect the quantity ∂|C ν ⊥ | 2 /∂ V and assume that the electric field affects only the out-of-plane component of the eigenfunction; we can write Let us now consider, in the spherical harmonics expansion, the tip states s, p z and d z 2 that are generally responsible for the main contribution to the tunneling current [18]. The matrix elements (2) corresponding to these states contain differential operators acting only on the z-coordinate and can be expressed in the functional form [20] where l and m are the quantum numbers related to the tip spherical harmonics expansion and α l m (E, E z , V ) depends on the particular tip state Here, k T (E, V ) = √ 1 + V /2 − E is the decay constant for the tip wave functions and a s , a p z , a d z 2 are constants related to the relative weight of the corresponding l m state in the spherical harmonics expansion. Now, introducing the l m -partial tip DOS ρ T,l m , obtained by summing over the tip states with the symmetry l m , we can rewrite the tunneling current as a sum of contributions I l m coming from each tip state symmetry l m = s, p z , d z 2 : Equation (9) can be expressed in a more transparent form by introducing the partial sample In (10), the tunneling current is expressed as an integral including L S,ν evaluated at the center of curvature of the tip apex. Usually, it is more interesting to consider this quantity evaluated at the sample surface. Making use of (4), it follows that identifying in this way the sample perturbed partial LDOS in vacuum as a product of an unperturbed partial LDOSs L S,ν (E, r ) evaluated at the surface and a tunneling decay function that depends on the bias and the tip-sample distance. This allows us to write (10) as where we introduce the function Equation (12) shows a number of similarities to the common 1D-WKB formula. It relates the tunneling current to the tip DOS and sample unperturbed partial LDOSs, each weighted by a generalized transmission coefficient τ l m . However, even though (12) has the same structure as that of the 1D formulation of the tunneling current, τ l m here depends not only on the total energy E of the surface states and the tip-sample distance z, but also on the out-of-plane energy E z and the symmetry of the tip apex l m , in this way including the relevant 3D features of the tunneling process. Comparing with the 1D-WKB transmission coefficient [8], τ l m still depends exponentially on z and V . However, the exponential term contains only the E z energy instead of E, which is the physical quantity contributing to the decaying behavior in vacuum of ψ ν [16,23]. Moreover, the 1D coefficient T does not contain the term α 2 l m , which arises naturally in this TH formulation and allows us to include the assumed symmetry of the orbitals at the tip apex.
Since the α l m coefficients are functions of V , E and E z , τ l m always differs from the 1D-WKB transmission coefficient. It is therefore important to understand how the transmission coefficients τ l m depend on the relevant quantities, in particular on the out-of-plane sample energy E z (see figure 1(a)). In figure 1(a), the total energy is fixed at a certain positive value, equal to the applied voltage (E = V ), and the behavior of τ l m as a function of E z is shown. The out-ofplane energy E z cannot overcome the total energy E and E z = E corresponds to a condition for which the in-plane motion energy is zero (i.e. the¯ point in the projected Brillouin zone). It can be observed that the transmission coefficient is overwhelmed by the exponential decay factor e −2k s z that strongly damps the tunneling probability of the states with low E z (the ones that lie far away from the¯ point) [23]. However, even though this effect is the dominant one, differences are visible when we compare the behavior of τ l m associated with different tip states. Looking at the functions α 2 l m in figure 1(b), it is clear that for p z and d z 2 states the tunneling probability is enhanced for surface states with lower E z or in other words, for surface states far from thē point, the probability does not depend on E z for s tip states. Furthermore, the transmission probability is always enhanced for the p z states with respect to the s states, whereas the d z 2 states damp the contribution when approaching¯ . Therefore, considering sample states with the same total energy, the tunneling is indeed dominated by those that are closer to the¯ point, but the symmetry of the tip state in general modulates the contributions of states belonging to different regions of the Brillouin zone. This aspect has to be properly considered when facing the issue of quantitatively extracting the sample LDOS from STS experiments.
We now consider the form of dI /dV , which represents the most relevant quantity in STS experiments, by directly differentiating (12): where B t,l m and B s,l m are given by In (14), the first term contains the most direct information on the surface electron properties. It represents a sum of every sample surface state having a total energy E = V 7 weighted by the transmission coefficient τ l m , which accounts both for the specific decaying behavior in vacuum of the sample wave functions and for the coupling with the tip states. The term B t,l m (in (15)) depends on the derivative of ρ T with respect to the energy and has to be considered only when the tip DOS is non-constant [9]. Experimentally, it is common to try to modify the tip apex electronic configuration by applying a strong voltage pulse during the measurement to achieve a sufficiently smooth tip DOS around the Fermi level. This modus operandi usually justifies the possibility of neglecting the contribution described in (15). In a recent work [9], different effects due to the presence of a non-constant tip DOS have been theoretically evaluated in the 1D framework and measured by experiments. Here, we want to focus mainly on the role of the symmetry of the tip apex orbitals in tunneling, and we will not further consider the effects due to its energetic structure, thus neglecting the term B t,l m in the following, leaving this for future work.
The term B s,l m (in (16)) involves a sum containing the derivative of the transmission coefficient τ l m , given by the following relations, where and consequently, by applying the mean value theorem [7] it is possible to write The termsβ l m ,ν andγ ν can be calculated for the case of constant tip and sample DOS by following the approaches proposed in the literature [7,9]. The term B s,l m has therefore been decomposed into a weighted sum of partial currents I l m ,ν , which represent the contributions to the total current of each sample state ν, tunneling into the tip state with the symmetry l m . It is a smooth term representing a background contribution to dI /dV , which has to be properly considered for a quantitative recovery of the surface LDOS, as will be discussed in section 4. In (10)- (14), a simple dependence of I on the total sample LDOS is by no means evident. In order to obtain a very clear connection between the tunneling current and the sample LDOS, we introduce a further assumption: The tunneling current I in (12) thus explicitly contains the term of interest L S (E, r ): and consequently, (14) becomes where B s,l m in (16), under condition (20), becomes Equation (22) now directly contains information on the sample LDOS L S at the energy E = V . L S is multiplied by the generalized transmission coefficientτ l m , evaluated at E = V and superposed on the background term B s,l m .
The assumption expressed in (20) is thus pointing to the possibility of a simplified analytical treatment and it is worthwhile to understand the conditions under which it holds. Looking at (6)-(8) and (13), we see that (20) is rigorously correct only if all the states described in the system have the same constant out-of-plane energy E z , which can occur in several cases of physical interest (see section 3). In the general case, the energy E z = E − E is different for different sample states. Consequently, the contribution of each sample state to the tunneling current cannot be simplified by considering only the quantity L S (E), and (12) has to be evaluated considering the full band structure of the investigated system.
If the sample electronic structure is in some way known (e.g. by means of ab initio calculations or analytical modeling), it is possible to use (9) or (10) to simulate an STS measurement depending on the assumed symmetry of the tip apex orbitals. This approach will be presented in section 3. On the other hand, if one is interested in recovering the surface LDOS from experimental STS data, (21) and (22) can be used for this purpose, as will be shown in section 4.

Tunneling spectroscopy of Shockley states
To provide an example of physical interest, we consider the case of noble metal surfaces such as Cu(111) and Au(111), which are characterized by a surface state confined in a few monolayers from the surface termination [24]. The presence of a projected gap at¯ produces a weak interaction between the surface and the bulk states [25] and, in these conditions, the surface states are well described by an ideal 2D free-electron gas regarding the in-plane motion. The possibility of treating the electronic properties of this system using an analytical formulation makes it a good candidate for applying the theoretical approach presented in section 2 and simulating the results of an STS measurement. For this system, a parabolic relation between energy and momentum is expected; neglecting all the deviations that affect the real system, we express such a relation in the form where E S is the bottom energy of the parabolic band (for Au and Cu, this value is estimated to be 460 and 430 meV below E F [26,27]), respectively, k is the momentum related to the in-plane motion and m * is the electron effective mass of the state (0. 23  Au and Cu, respectively [26,27]). It is worth noting that the out-of-plane component of the energy E z does not depend on the momentum in this case, so that E z = E S for every k . For this system, the use of (20)- (21) is then fully justified. The DOS is obtained by integrating over the momentum space. For a full 2D translational symmetry, a step function with the discontinuity placed at E = E S is obtained. In STS, this step is usually visible as a jump in dI /dV signal for V = E S [25][26][27]. In the presence of a scattering center like point defects or steps, the translational symmetry is partially broken (scattering sites interact with the electrons of the surface and generate the well-known phenomenon of standing waves [28]) and the proper physical quantity to be considered is the LDOS instead of the DOS. The LDOS is modulated both in space and in energy, as can be observed routinely in low-temperature STM experiments.
In order to include such effects, we describe a model system of a 2D Shockley state, with an ideal step barrier along the x-direction. In this configuration, the eigenfunctions of the system are [29] where g k x (x) = L −1/2 sin(k x x); h k y (y) = L −1/2 e ik y y ; χ(z) = Ce −k S z , and k 2 x + k 2 y = k 2 . The corresponding LDOS is The oscillations of the LDOS in energy are then described by the Bessel function J 0 (see figure 2). We apply the tunneling theory described in section 2 to this system to perform a full analytical simulation of the tunneling current. Simulated I -V characteristics at typical tunneling parameters are shown in figure 3(a), considering the cases of s, p z and d z 2 tip states. As the common behavior of all these curves, the absolute value of the tunneling current is larger in the negative bias region, whereas in the positive region it approaches a zero slope limit. Due to the integral nature of the current, we can observe only a very weak contribution coming from the oscillations of the LDOS. In figure 3 the oscillations in energy. The simulated dI /dV s agree well with typical experimental data on Shockley surfaces, such as that in figure 6 acquired on Cu(111), in particular regarding the monotonically decreasing behavior and the sharp onset. The maxima of these oscillations in the dI /dV curves are close to the L S ones but a shift towards lower energy is present. In this sense, our simulation supports the well-known problems in extracting physical quantities of the sample directly from the dI /dV data, without further treatment. More precisely, we want to point out that the shape of the dI /dV is both qualitatively and quantitatively different from that of the sample LDOS. From a comparison of figures 2 and 3(b), we can see that (i) the differential conductivity shows a decreasing background distortion and it can even become negative after a certain value of the bias ( figure 3(b)), and (ii) for negative biases, dI /dV continuously increases even where the LDOS evaluated for E = V is zero. Furthermore, the behavior of dI /dV is affected in a non-obvious way by the symmetry of the tip apex orbitals. In order to better understand the origin of these spurious contributions, we now interpret the shape of the dI /dV curves in the light of (22). Let us consider the decreasing behavior of the dI /dV signal. This effect can be usually ascribed to the energy behavior of the tunneling barrier and has been exhaustively explained in the framework of 1D-WKB theory [8,9]. However, in the case of Shockley states, this effect is not expected simulating the tunneling using a 1D model [7], so that for this specific case its origin has to be ascribed to the 3D features included in the transmission coefficient. In order to elucidate this point, the quantityτ l m (E, V, E z , z)| E=V is shown in figure 3(c) as a function of V for a constant E z . We can see that in any case it shows completely different behavior compared with the 1D-WKB transmission coefficient [6]. For the Shockley states, the tunneling probability decreases when the applied bias increases. This is due to the fact that each state with a total energy E has a fixed out-of-plane energy E z , and increasing the applied voltage only has the effect of increasing the barrier height, which determines the observed decrease in the transmission coefficient and the damping of the differential conductivity. On the other hand, for the general case of surface or surface-projected states, it is likely that a higher total energy E results in a higher out-of-plane energy E z and a higher tunneling probability. Therefore, increasing the applied voltage not only increases the barrier height but also allows tunneling among states with a higher energy and tunneling probability; thus the resulting voltage behavior is the opposite. Equation (22) is also helpful in understanding why we observe a non-vanishing dI /dV below the edge of the Shockley state. If we consider the term B s,l m of this expression, we see that it is proportional to the tunneling current, which, being non-vanishing at negative biases, provides a finite contribution also in the energetic interval in which the surface has no states.
We now want to highlight the effect of the symmetry of the tip apex orbitals. From figure 3(b), we see that the s, p z and d z 2 tip symmetries generate three different dI /dV curves with different amplitudes. In particular, the s and p z cases are very similar, whereas for d z 2 the amplitude of dI /dV is reduced. The differences in their amplitude can be once again understood by looking at the behavior ofτ l m in figure 3(c). The effect of the tip apex orbital symmetry is contained inτ l m , which depends on the tunneling matrix elements. Since the coefficient τ d z 2 is always smaller than the other two coefficients, this results in decreased values of tunneling current and differential conductivity for the d z 2 states, as visible in figure 3(b). To summarize the previous discussion, we observed that the main features of the Shockley state can be directly understood by looking at the dI /dV curve. However, various spurious contributions from the tunneling process arise and they have to be properly removed for a quantitative recovery of the sample LDOS.

Recovery of the sample local density of states
Let us now discuss the crucial topic of how to extract quantitative information on the sample LDOS from the tunneling current. If the tunneling is dominated by one tip state symmetry in the energy range of interest (I l m ≡ I ), an explicit expression for the sample LDOS in terms of experimental and modeling parameters can be obtained combining (22) and (23) for systems satisfying the conditions (20): Equation (27) is expressed in terms of the tunneling current I , in a similar way to that suggested by Koslowski et al in the 1D case [7]. The parametersβ l m andγ are functions of V . Equation (27) represents an explicit expression to obtain the sample LDOS as a function of the measured quantities dI /dV and I and simple analytical functions. In the general case, we may expect that more than one tip state contributes to the tunneling, so that the terms ρ T,l m ,τ l m and β l m should consequently be substituted by proper mean quantities. However, the assumption of one symmetry tunneling represents a simple and usually the best solution from a practical point of view.
We are therefore led to evaluate the effectiveness of (27) on the physical system discussed in section 3. For this purpose, a comparison of (27) with the normalization methods proposed in the literature by Stroscio et al [5], Ukraintsev [6] and Koslowski et al [7] will be presented below, followed by an analysis of the sensitivity to the modeling parameters (i.e. tip-sample distance z and tip orbital symmetry l m ). Then, (27) will be also applied to experimental data acquired regarding Cu(111). In these investigations, both the presented method and those from the literature have been applied to the I and dI /dV curves simulated in the case of the s tip. Except normalization to total conductivity I /V , all these methods require an estimate of the tunneling parameters such as tip-sample distance and effective work-function [9]. For the sake of simplicity, here these parameters have been directly taken from the simulations. As far as the Koslowski method is concerned, we directly use the 3D-WKB extended method (equation (13) of [7]), specifically suited for the Shockley state system.
The comparison is shown in figure 4. All the previous methods preserve an overall background distortion represented by a negative slope in all the recovered spectra. This distortion also affects the positions of the maxima of the oscillation in the spectra, producing a shift towards lower energy values. Moreover, some other shape distortions are visible in the spectra normalized using the total conductivity method and the Ukraintsev method. To explain such effects, we can observe that these procedures have been developed in a 1D framework, while this peculiar system presents important features that can be taken into account only in a 3D treatment. Even using Koslowski normalization, which is a 3D-WKB extension of the corresponding 1D-WKB method, it is not possible to completely remove the overall distortion, mainly because a proper inclusion of the tunneling matrix elements between the sample and the tip is missing. From figure 4, it is clear that the best result is achieved by applying (27).
The proposed method is also capable of properly removing the effects of different tip states. In figure 5(a), we show the results of applying (27) to the three dI /dV and I obtained by simulating the tunneling through the s, p z and d z 2 tip states making use of the corresponding β, γ and τ functions. Even though, below the Shockley edge, the recovered curves show some small differences, all three normalized curves are coincident above this edge and in full agreement with the sample LDOS, thus showing the capability of the procedure to remove the influence of tip apex orbitals symmetry and achieve the correct result. It might appear that the applicability of this method relies on a detailed knowledge of the tip apex orbital symmetry and the tip-sample distance, which are often not achievable from a practical point of view. In order to evaluate possible errors introduced by wrong estimates of these experimental parameters, we consider the behavior of the proposed procedure when a 'wrong' estimated tip state is used, as shown in figure 5(b). Here, a recovered dI /dV curve related to a p z state tip is shown, together with the application of (27) also assuming other tip symmetries. Some background distortion can be introduced if the tip apex orbital symmetry differs from that responsible for the tunneling process, but such an effect is anyhow clearly less pronounced when compared to that obtained using conventional 1D methods. In other words, it is better to normalize using (27) with a trial matrix element (possibly, even an incorrect one) instead of applying a method that does not consider the tip apex symmetry at all. We similarly consider the problem of a wrong estimation of the tip-sample distance. In figure 5(c), the dI /dV curve is simulated by considering an s state tip located 10 Å from the surface, and (27) was applied using three different estimated distances. Once again, the overall effect of using a wrong parameter is a background distortion, particularly visible in the change of slope of the recovered LDOS at voltages above the surface state onset. Since both the tip state and the tip-sample distance parameters act in a similar way on the background of the recovered spectra, it might be possible to obtain similar results by applying (27) using different combinations of estimated parameters. As an example of this, we show in figure 5(d) the results of applying (27) to the simulated STS curve with l m = s and z = 10 Å using both the correct parameters and two different combinations of estimated l m and z, namely l m = p z , z = 11 Å and l m = d z 2 , z = 13 Å. In all these three cases, the shape of the recovered LDOS is very close to that of the ideal surface LDOS (see figure 2 for a comparison), and only small errors are introduced because of using the improper set of parameters. We now conclude this section by applying the procedure to experimental STS data acquired regarding the Cu(111) surface. Measurements were made at 7 K using an Omicron lowtemperature STM on a conventionally prepared Cu(111) single crystal [27] and using a bulk Cr tip [30]. In figures 6(a) and (b) the dI /dV and I curves are presented. One can apply (27) in order to extract the measurement parameter l m and z from the STS data. Since real measurements may contain several sources of spurious contributions not considered in this work (e.g. surface-projected bulk band states and non-constant tip DOS effects [9]), it is likely that complete agreement between recovered experimental and theoretical LDOSs could not be expected over the whole energy range of the experiments. Therefore, we focus our attention on the energy region close to the Shockley onset, which represents the most reproducible and reliable feature in STS measurements of this surface. We varied the set of parameters l m and z until finding satisfactory agreement between theoretical and recovered LDOS in the range of −0.5-0 V. Figure 6(c) shows the result achieved using l m = s and z = 10 Å. In agreement with what was obtained on simulated spectra, a similar result can be also achieved using other sets of parameters, such as l m = p z , z = 11 Å and l m = d z 2 , z = 13 Å. Although none of these values can be excluded a priori, we are here led to choose the s tip case as the most probable one, being the only one which corresponds to the typical z range of 5-10 Å according to the set point usually adopted in STS measurements for this class of surfaces [12,31,32]. Comparing figures 2 and 6(a), we see two effects that do not appear in our simulated spectra outside the reference range, namely a wide feature below the Shockley edge and an increasing background of the recovered LDOS curve as a function of the applied voltage above +0.1 V. These features are not due to our procedure since we found that they are stably detectable also employing other l m values and varying z by ±2 Å. Since no sample states, neither from the surface nor from the projected bulk band, are expected within the energy range of −0.8-0.45 eV [33], we safely attributed the feature below the Shockley onset to a non-constant tip DOS effect [9]. The rising signal above Fermi energy suggests instead a different origin, probably connected with a change of the decaying behavior of the sample states in vacuum as a function of energy. It is worth recalling that Cu surface states are located within the surface projected bulk band gap below Fermi energy, whereas, with increasing momentum as a function of energy, they cross the bulk band edge just above the Fermi level [34], losing in this way its pure surface character and becoming instead a surface resonance [29]. As a consequence, the decaying behavior in vacuum is partially modified by the bulk contributions and the assumption that E z = E S is not valid in this situation. This might therefore explain the different behavior of the recovered L S as a function of energy with respect to pure Shockley states LDOS. It is remarkable that such a feature is particularly evident when STS data are treated by the proposed procedure, since (27) relies on the presence of a pure surface state and not on possible bulk band state or surface resonances, so that possible deviations from the ideal free-electron behavior are particularly emphasized.
The presented analysis suggests two possible ways of applying (27) to real situations. First, one can be interested in normalizing STS spectra from unknown surfaces using some reasonable trial parameters. In this case, we showed that it is possible to achieve a proper recovery of the LDOS by using different sets of parameters, which makes the procedure even more reliable. On the other hand, the proposed procedure could be applied to STS data from well-known surfaces (such as Shockley states) with the aim of extracting tip apex orbital symmetry and experimental parameters. This can be implemented by tentatively varying the parameters until the recovered LDOS (or alternatively a selected part of the recovered LDOS in a certain energy window) is very similar in shape to the ideal form of the surface LDOS. Since this result can be achieved through more than one combination of estimated parameters, a unique set could not be univocally found. However, it is also possible that some of the extracted parameter sets could be excluded if some values are not physically reasonable (e.g. too short/long tip-sample distances or unexpected tip apex symmetry for the used tip material), achieving in this way a unique result.

Conclusions
In this paper, we have presented a 3D-TH model of tunneling that is capable of including both the tip angular symmetry and the sample band structure. We found that the tunneling current can be represented as an energy convolution between the tip DOS and the sample LDOS weighted by an effective transmission coefficient τ . In contrast to what occurs in the 1D-WKB description, this coefficient depends on the 3D properties of the tip and the sample, in particular on the out-of-plane sample state energy (E z ) and the tip state angular symmetry (s, p z or d z 2 ). On the basis of these findings, we proposed a specific normalization procedure to properly recover the sample LDOS from an STS measurement. We considered the case of a Shockley state, which is analytically developed to also compare this normalization method with other relevant existing approaches. This analysis clearly shows the importance of including the 3D effects in the description of tunneling in order to better quantitatively recover the sample LDOS from experimental measurements.
The presented theory shows several advantages with respect to the other 1D-WKB-based models. First, the theory can be used to simulate the dI /dV curves starting from a theoretical sample electronic structure in order to reveal how an STS measurement should appear. This would be of great help in overcoming the usual difficulties faced in comparing the ab initio simulations of surface electronic properties with STS data. With respect to WKB theories, the decaying behavior in vacuum of states belonging to different regions of the Brillouin zone can now be properly considered, together with effects introduced by the tip apex orbital symmetry. When the starting point is the experimental data, the proposed normalization procedure can be exploited to recover (i) the sample LDOS from STS measurements assuming a constant tip DOS or (ii) the tip apex orbital symmetry and electronic structure when STS is performed on well-known surfaces. In principle, the presented method for the recovery of tip and sample LDOS could also be exploited in connection with iterative methods, as proposed in the 1D-WKB framework [7,10], e.g. by using a set of measurements taken at different tip-sample distances [11], providing also in this case a correct inclusion of the 3D effects.