# STM contrast of a CO dimer on a Cu(111) surface: A wave-function analysis

## Abstract

We present a method used to intuitively interpret the STM contrast by investigating individual wave functions originating from the substrate and tip side. We use localized basis orbital density functional theory, and propagate the wave functions into the vacuum region at a real-space grid, including averaging over the lateral reciprocal space. Optimization by means of the method of Lagrange multipliers is implemented to perform a unitary transformation of the wave functions in the middle of the vacuum region. The method enables (i) reduction of the number of contributing tip-substrate wave function combinations used in the corresponding transmission matrix, and (ii) to bundle up wave functions with similar symmetry in the lateral plane, so that (iii) an intuitive understanding of the STM contrast can be achieved. The theory is applied to a CO dimer adsorbed on a Cu(111) surface scanned by a single-atom Cu tip, whose STM image is discussed in detail by the outlined method.

###### pacs:

68.37.Ef, 33.20.Tp, 68.35.Ja, 68.43.Pq## I Introduction

Scanning tunneling microscopy (STM) is a mature technique to reveal atomic structures on surfaces. In addition to images of atomic arrangements, detailed measurements can be made of molecular orbitals Repp et al. [2005]; Gross et al. [2011], electronic structure, vibrational Stipe et al. [1998] and magnetic excitations Hijibehedin et al. [2007]. With atomic force microscoscopy (AFM), the interaction between tip and sample is characterized, which allows for examination of the tip-apex structure Ternes et al. [2008]; Emmrich et al. [2015a]; Welker and Giessibl [2012]; Emmrich et al. [2015b]. The tip-apex structure may strongly influence the tunneling process, both when making a quantitative comparison between experimental and theoretical STM images Gustafsson et al. [2017], and when investigating the inelastic tunneling signal from an adsorbate species Okabayashi et al. [2016].

To model STM images, the approximation methods of Bardeen Bardeen [1961] and Tersoff-Hamann Tersoff and Hamann [1985] have been widely used. The Tersoff-Hamann approach, which can be derived from the Bardeen’s approximation with an -wave tip Hofer et al. [2003], has provided an intuitive understanding of STM experiments with non-functionalized STM tips Persson and Ueba [2002]. For CO-functionalized STM, the Bardeen method Paz and Soler [2006]; Rossen et al. [2013]; Bocquet and Sautet [1996]; Teobaldi et al. [2007]; Zhang et al. [2014], the Chen’s derivative rule Chen [1990]; Mándi and Palotás [2015]; Mándi et al. [2015], and Landauer-based Green’s function methods Cerdá et al. [1997a, b] include the effects of more complicated tip states.

Our previous work on this topic Gustafsson and Paulsson [2016]; Gustafsson et al. [2017] concerns first-principles modeling based on Bardeen’s approximation, using localized-basis DFT. The wave functions close to the atoms were calculated in a localized basis set using Green’s function techniques, whereafter they were propagated into the vacuum region in real space by utilizing the total DFT potential. We concluded that calculations in the -point () may, for many systems, qualitatively reproduce the STM contrast of the k averaged calculations, whereas a quantitative comparison to experiments always seems to require averaging over the lateral reciprocal space Gustafsson et al. [2017].

In this paper we further develop the theoretical method focusing on providing a simpler interpretations of the STM contrast. The numerous wave functions given by DFT in the vacuum region obscures the interpretation and here we have developed a unitary transformation to find the dominating tip- and substrate wave function combinations (henceforth denoted - combinations) that give the largest contribution to a specific STM image. We also present a simple formula, relating the real-space tip- and substrate wave functions to the transmission probability. To illustrate the method, we analyse the STM contrast of a CO dimer Heinrich et al. [2002]; Zöphel et al. [1999]; Niemi and Nieminen [2004] adsorbed on a Cu(111) surface scanned by a single-atom Cu tip.

## Ii Methodology

### ii.1 Transmission coefficient

In our previous development Gustafsson and Paulsson [2016]; Gustafsson et al. [2017] of calculating STM images from localized-orbital DFT, we used the Bardeen’s approximation Bardeen [1961] for the real-space propagated wave functions in the vacuum region. The wave functions are calculated on an equidistant discrete lattice utilizing a simple nearest-neighbour finite difference (FD) scheme. In this approximation, the coupling strength operator (matrix), , between individual lattice planes along the direction of transport () therefore is a diagonal matrix ^{1}

(1) |

This equation can be derived from the more familiar expression for the transmission coefficient (henceforth omitting the superscript ), , which relates the Green’s function of the device region (i.e., the substrate wave functions when connected to the tip part), to the isolated tip part by the broadening ,

(2) | ||||

(3) |

from which Eq. (1) follows. In the derivation, the broadening of the electronic states from the tip side, , is expressed by utilizing the partial spectral function, , so that , which in turn is composed by the totally reflected tip wave functions, (see Ref. Paulsson and Brandbyge [2007] for further details). This means that the total transmission is proportional to the sum of the squares of the overlap of all - combinations, due to the utilized FD approximation for the wave functions.

Upon using the Bardeen’s approximation we introduce a separation surface Gustafsson and Paulsson [2016]; Gustafsson et al. [2017], , in vicinity of the middle of the vacuum gap (where the total DFT potential reaches its maximum) and set the potential away from this slice as the average of the potential at . The wave functions from the substrate, , and the tip, , are then computed equivalently, but separately from each side of the vacuum region.

An important note is that the substrate wave functions used in the Bardeen formula and in Eq. (1) are exactly the same, whereas the tip wave functions in Eq. (1), , are totally reflected at the separation surface, in contrast to the Bardeen formula. The isolation of the tip side is modeled by introducing an impenetrable barrier at the separation surface when propagating the wave functions from the tip side. This means that the totally reflected tip wave functions, in all essentials, are identical to the normally propagated ones, apart from the amplitude, which depends on the lattice constant . This subtle difference in the calculation of the wave functions enables usage of Eq. (1) instead of the Bardeen formula, and the two methods give identical results with same memory consumption in the same CPU time. However, calculation of the derivatives of the wave functions (needed in Bardeen’s approximation) is not necessary in Eq. (1), and the suggested formula provides shorter and simpler derivations upon considering the forthcoming unitary transformation of the wave functions.

### ii.2 Unitary transformation of wave functions

Below we have approximately 1000 - combinations, and to interpret the calculated STM contrast we describe a simple method to perform a unitary transformation of the wave functions. The main idea is to maximize on the separation surface (and similarly for the tip wave functions), so that the important wave functions from each side of the vacuum region become ordered by their amplitude, which in turn is nearly proportional to their importance as tunneling channels. Maximization with respect to the conductance, i.e., , may be appropriate when considering a fixed tip position, where each - combination gives a single-valued conductance. However, since we consider a tip scanning over the whole substrate, the former maximization is more relevant in the present context.

The th unitary transformed wave function can therefore be written as a linear combination of the accessible propagated wave functions,

(4) |

subjected to the the constraint for the expansion coefficients, , where is the number of considered wave functions from each side. By means of the Lagrange multiplier method Arfken et al. [2012]; Bertsekas [1982] within the principal component analysis Pearson [1901]; Jolliffe [2002], the equation to solve (for each ) reads,

(5) |

upon considering transformation of the substrate wave functions. This equation reduces to a standard eigenvalue problem,

(6) |

where is the matrix containing all overlap integrals between individual wave functions. The desired expansion coefficients, , are therefore directly obtained by diagonalization of , Eq. (6), so that the numerical implementation is trivial. This procedure is performed similarly for the totally reflected wave functions from the tip side, so that a different transformation matrix is obtained for these wave functions, unless the substrate and tip slabs are identical.

An important feature of the outlined method is that the sum of the overlaps of all - combinations squared, Eq. (1), is invariant under the unitary transformation, Eq. (4). This property may conveniently be proved backwards, by expansion of the unitary transformed substrate wave functions,

(7) | |||

(8) |

where the orthogonality of the eigenvectors, , is exploited, so that the mix terms vanish. The same procedure applied to the unitary transformed tip wave functions, , completes the proof. This shows that the total transmission, Eq. (1), (for a specific point) remains invariant under a unitary transformation, and ultimately that the k averaged transmission coefficient is unaffected by such a transformation. Hence, the unitary transformation matrices, , are found, such that , leaving the physical properties unchanged.

This procedure often seems to result in that wave functions of the same lateral symmetry, i.e., highly correlated variables, are collected and thereby reduce the dimensionality of the problem. This means that a smaller subset of wave functions are needed in the transmission matrix to obtain an accurate description of the STM image.

### ii.3 Computational details

We use the Siesta Soler et al. [2002] DFT code to geometry optimize a slab consisting of eight Cu layers, where each layer has Cu atoms with a nearest-neighbour distance of 2.57 Å, so that the lateral cell dimensions are Å. In the calculations presented below we use a Å core-core distance along between the very apex of the Cu tip and the oxygen atoms (see Fig. 1), which assures a negligible interaction between the tip and the substrate Gustafsson and Paulsson [2016], so that the substrate wave functions are unaffected when simulating the scanning of the STM tip over the surface.

The Siesta calculations are performed using the Perdew-Burke-Ernzerhof (PBE) parametrization of the generalized gradient approximation (GGA) exchange-correlation functional Perdew et al. [1996], double- (single-) zeta polarized basis set for C, O (Cu) atoms, a 200 Ry real space mesh cutoff, and points in the surface plane. The CO molecules are adsorbed on two adjacent top site of the otherwise clean Cu(111) substrate surface, and on the opposite side of the vacuum region a pyramid tip consisting of four Cu atoms are attached to a similar surface, see Fig. 1. The geometry optimization concerns the adsorbate species, the two top substrate layers, and the tip atoms (forces less than 0.04 eV/Å). Nine additional Cu layers are thereafter added (three at the bottom of the substrate, and six above the tip slab), and the wave functions are calculated by Transiesta Brandbyge et al. [2002], Inelastica Frederiksen et al. [2007], and our recent STM model Gustafsson and Paulsson [2016]; Gustafsson et al. [2017].

The k grid used in the STM simulation calculation consists of k points that are homogeneously distributed in reciprocal space, and shifted so that an odd number of k points includes the point. In the following analysis, the wave functions obtained from a -point STM calculation are discussed. It is therefore important that the considered k averaged STM image agrees qualitatively with the -point calculation, unless there is interest in investigating several k points.

## Iii Results

### iii.1 Calculated STM contrast of a CO dimer

A single CO molecule prefers to stand upright when adsorbed on a top site of a Cu(111) surface, i.e., the C-O bonding axis is perpendicular to the surface. Such an adsorbate exhibits a solid radially symmetric conductance dip when scanned by a pure Cu tip Gustafsson and Paulsson [2016]; Gustafsson et al. [2017].

The present CO dimer, i.e., two CO molecules adsorbed on adjacent top sites of a Cu(111) surface, reveals a bright spot centered in a slightly elongated surrounding dip, when scanned in constant-current mode by an -wave tip Heinrich et al. [2002]. A similar feature is also experimentally observed for multiple CO monomers on a Cu(211) surface Zöphel et al. [1999]. Furthermore, the CO dimer exhibits slightly tilted bonding angles of the C-O bonding axes, as the oxygen atoms push apart, which has an evident effect of the conductance perpendicular to the surface Niemi and Nieminen [2004]. We confirm this feature [Fig. 1(c)] in constant-height mode, when scanning the molecule by a four-atom pyramidal Cu tip at tip height 7.5 Å, averaging over a relatively dense () k grid. The CO bonding length is 1.17 Å, and the tilt angles of the C-O bonding axes deviate from standing perpendicular to the Cu(111) surface, in agreement to previous studies Persson [2004]. We have also noticed that the calculated STM image for this system weakly depends on the tip height. By lowering the tip by 1.5 Å, the central bright spot becomes more pronounced, whereas elevating the tip by 1.5 Å, lowers the central bright spot.

For this, and other systems Gustafsson et al. [2017], a -point STM calculation [Fig. 1(a)] qualitatively reproduces the k averaged STM image [Fig. 1(c)], which allows to restrict the forthcoming wave function analysis to the point. The main difference between these calculations is an overall increment of conductance, as well as a larger contribution from the Cu substrate surface atoms, when performing averaging in reciprocal space. This feature is visualized in Fig. 2, where the conductance in vicinity of the point is smaller than closer to the edge of the surface Brillouin zone. The unit [pA/V] for the conductance used in Fig. 1 is acquired by conventionally relating the transmission probability to the conductance with the formula , where is the conductance quantum, and is the transmission coefficient calculated by Eq. (1). According to Fig. 2(b), we conclude that also k points further away from the Brillouin zone center give -point-like STM images. Therefore, by choosing a k grid that excludes the point does not have an impact on the k averaged STM contrast. For instance, a coarser k grid () that excludes the point, gives an almost identical STM image compared to the grid used here.

### iii.2 Unitary transformation of wave functions

Despite that the considered STM image is already calculated and briefly discussed in reciprocal space, the origin of the STM contrast still needs to be interpreted. Each k-point STM image in Fig. 2 is a result of approximately 1000 () individual - combinations used in Eq. (1). The significance of these combinations to the STM contrast may be examined directly from the propagated wave functions. However, the unitary transformation of these wave functions, described in Sec. II.2, yields a more transparent understanding of the STM contrast.

In Fig. 3(a) and Fig. 3(b) the real part ^{2}

When performing a unitary transformation of these wave functions, the ordering approximately becomes proportional to the maximum amplitude of the wave function, as well as their significance in the tunneling probability for a specific - combination. This is visualized in Fig. 3(d) and Fig. 3(e), and clearly reflected in the transmission matrix, Fig. 3(f), where large transmission coefficients are observed primarily for the low-number wave functions from each side. The magnitude of a specific element in the transmission matrices is defined as the total transmission of the given - combination. Furthermore, wave functions with the same symmetry in the lateral plane are bundled up in a single wave function with the same symmetry. For instance, the most pronounced unitary transformed Cu-tip wave function, Fig. 3(e), is one single wave, which can be compared to the normal Cu-tip wave functions, Fig. 3(b), where several wave functions with this symmetry are present. A cutoff is used in the plot range for the wave functions in Fig. 3, in order to highlight also the less important wave functions. For instance, the doubly degenerate lateral waves from the tip are essentially unimportant in respect to the STM contrast, see Fig. 4(e) and Fig. 4(f).

As in the case of a single CO molecule adsorbed on a Cu(111) surface, we have previously shown that the nonzero amplitude of certain substrate wave functions away from the molecule is crucial to explain the conductance dip over the molecule Gustafsson and Paulsson [2016]; Gustafsson et al. [2017]. That is, if a wave function over the substrate has a significant (constant) value that yields a large conductance even when the tip is laterally placed far away from the molecule, this might give a conductance dip over the molecule, as in the case with the CO monomer on the Cu(111) surface. This sign change in the amplitude is evident also for the present CO dimer upon noticing the nonzero amplitude of the first two unitary transformed substrate wave functions [Fig. 3(d), wave function 0 and 1], which, by interference, give rise to a large current also when the tip is not centered over the molecule. This feature is the main reason to the slightly elongated conductance dip in the STM image, see Fig. 1(b), whereas its central bright spot is assigned to the large amplitude centered in the first unitary transformed wave function of the substrate [Fig. 3(d), wave function 0].

### iii.3 Reduction of tunneling channels

Another important characteristic of a unitary transformation of the wave functions, is that it enables to decrease the number of necessary - combinations that are must be included to recover the STM contrast, as a consequence of the ordering by the amplitude in the vacuum gap for wave functions of the same symmetry. This is demonstrated in Fig. 4.

We compare to the correct -point STM image by introducing a cutoff conductance, so that, upon defining X%-accuracy, where the chosen cutoff gives X% of the total current of the correct -point STM image. For instance, when using the ordinary wave functions [Fig. 3(a) and Fig. 3(b)], one needs 152 - combinations [14% out of the 1089 () combinations], to obtain a 99% accuracy. Achieving the same accuracy with the unitary transformed wave functions [Fig. 3(d) and Fig. 3(e)] only requires 26 combinations [2.3% out of 1089 channels]. In the latter case, the majority of the current is carried by one single -type tip wave function, as expected by a single-atom Cu tip.

Reducing the number of channels further (by increment of the cutoff conductance), by compromising slightly with the quality of the STM image [Fig. 4(b) and Fig. 4(e)] shows the same tendency, where the current is solely carried via the Cu-tip wave. Notice that such a calculation results in an STM contrast that could equivalently be obtained by the Tersoff-Hamann approximation Tersoff and Hamann [1985], due to the pure -wave character of the tip. By using an even larger cutoff, so that only two - combinations are used [tip wave function 0 ( type) Fig. 3(e), and two substrate wave functions, 0 and 1, Fig. 3(d)], may reproduce the -point STM image qualitatively, so that the origin of the STM contrast is easy to interpret, see dashed line in Fig. 4(f). The latter calculation, however, only gives a 59% accuracy, which originates from significantly less contribution from the Cu surface.

## Iv Summary

We have shown that a unitary transformation, by means of the method of Lagrange multipliers, offers a simplified picture when resolving the origin of the STM contrast in terms of individual tip-substrate wave function combinations in the middle of the vacuum region. The unitary transformed wave functions become nicely ordered by (i) their amplitude in the vacuum region, which approximately (ii) determine their importance in the tunneling process, and (iii) are merged into single versions of wave functions of similar symmetry in the lateral plane. The method significantly reduces the number of tunneling channels, and thereby opens up the possibility to make an intuitive and detailed analysis of individual wave function combinations of a specific lateral symmetry. We have further presented an alternative simple formula originating from the Green’s function formalism, to describe the transmission probability accounting for multiple-tip-state STM modeling when considering real-space wave functions in the vacuum gap.

## V Acknowledgements

The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at Lunarc. A.G. and M.P. are supported by a grant from the Swedish Research Council (621-2010-3762).

### Footnotes

- In contrast to localized-basis calculations, where the coupling strength depends on the tip position, which therefore varies when the tip is shifted laterally over the substrate.
- The imaginary part is negligible in the point for large vacuum gaps.

### References

- J. Repp, G. Meyer, S. M. Stojković, A. Gourdon, and C. Joachim, Phys. Rev. Lett. 94, 026803 (2005).
- L. Gross, N. Moll, F. Mohn, A. Curioni, G. Meyer, F. Hanke, and M. Persson, Phys. Rev. Lett. 107, 086101 (2011).
- B. C. Stipe, M. A. Rezaei, and W. Ho, Science 280, 1732 (1998).
- C. F. Hijibehedin, C.-Y. Lin, A. F. Otte, M. Ternes, C. P. Lutz, B. A. Jones, and A. J. Heinrich, Science 317, 1199 (2007).
- M. Ternes, C. P. Lutz, C. F. Hirjibehedin, F. J. Giessibl, and A. J. Heinrich, Science 319, 1066 (2008).
- M. Emmrich, M. Schneiderbauer, F. Huber, A. J. Weymouth, N. Okabayashi, and F. J. Giessibl, Phys. Rev. Lett. 114, 146101 (2015a).
- J. Welker and F. J. Giessibl, Science 336, 444 (2012).
- M. Emmrich, F. Huber, F. Pielmeier, J. Welker, T. Hofmann, M. Schneiderbauer, D. Meuer, S. Polesya, S. Mankovsky, D. Ködderitzsch, et al., Science 348, 308 (2015b).
- A. Gustafsson, N. Okabayashi, A. Peronio, F. J. Giessibl, and M. Paulsson, Phys. Rev. B 96, 085415 (2017).
- N. Okabayashi, A. Gustafsson, A. Peronio, M. Paulsson, T. Arai, and F. J. Giessibl, Phys. Rev. B 93, 165415 (2016).
- J. Bardeen, Phys. Rev. Lett. 6, 57 (1961).
- J. Tersoff and D. R. Hamann, Phys. Rev. B 31, 805 (1985).
- W. A. Hofer, A. S. Foster, and A. L. Shluger, Rev. Mod. Phys. 75, 1287 (2003).
- B. Persson and H. Ueba, Surf. Sci. 502-503, 18 (2002).
- O. Paz and J. M. Soler, Phys. Status Solidi B 243, 1080 (2006).
- E. T. R. Rossen, C. F. J. Flipse, and J. I. Cerdá, Phys. Rev. B 87, 235412 (2013).
- M.-L. Bocquet and P. Sautet, Surf. Sci. 360, 128 (1996).
- G. Teobaldi, M. Penalba, A. Arnau, N. Lorente, and W. A. Hofer, Phys. Rev. B 76, 235407 (2007).
- R. Zhang, Z. Hu, B. Li, and J. Yang, J. Phys. Chem. A 118, 8953 (2014).
- C. J. Chen, Phys. Rev. B 42, 8841 (1990).
- G. Mándi and K. Palotás, Phys. Rev. B 91, 165406 (2015).
- G. Mándi, G. Teobaldi, and K. Palotás, Prog. Surf. Sci. 90, 223 (2015).
- J. Cerdá, M. A. Van Hove, P. Sautet, and M. Salmeron, Phys. Rev. B 56, 15885 (1997a).
- J. Cerdá, A. Yoon, M. A. Van Hove, P. Sautet, M. Salmeron, and G. A. Somorjai, Phys. Rev. B 56, 15900 (1997b).
- A. Gustafsson and M. Paulsson, Phys. Rev. B 93, 115434 (2016).
- A. J. Heinrich, C. P. Lutz, J. A. Gupta, and D. M. Eigler, Science 298, 1381 (2002).
- S. Zöphel, J. Repp, G. Meyer, and K. H. Reider, Chem. Phys. Lett. 310, 145 (1999).
- E. Niemi and J. A. Nieminen, Chem. Phys. Lett. 397, 200 (2004).
- M. Paulsson and M. Brandbyge, Phys. Rev. B 76, 115117 (2007).
- G. B. Arfken, H. J. Weber, and F. E. Harris, Mathematical Methods for Physicists, 7th Edition (Elsevier, 2012), ISBN 978-0-12-384654-9.
- D. P. Bertsekas, Constrained Optimization and Lagrange Multiplier Methods (Academic Press, 1982), ISBN 1-886529-04-3.
- K. Pearson, Philosophical Magazine 2, 559 (1901).
- I. T. Jolliffe, Principal Component Analysis, Second Edition (Springer-Verlag New York, Inc., 2002), ISBN 0-387-95442-2.
- J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal, J. Phys.: Condens. Matter 14, 2745 (2002).
- J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- M. Brandbyge, J.-L. Mozos, P. Ordejón, J. Taylor, and K. Stokbro, Phys. Rev. B 65, 165401 (2002).
- T. Frederiksen, M. Paulsson, M. Brandbyge, and A.-P. Jauho, Phys. Rev. B 75, 205413 (2007).
- M. Persson, Phil. Trans. R. Soc., Lond. A 362, 1173 (2004).