A simple model of molecular imaging with noncontact atomic force microscopy

Using functionalized tips, the atomic resolution of a single organic molecule can be achieved by noncontact atomic force microscopy (nc-AFM) operating in the regime of short-ranged repulsive Pauli forces. To theoretically describe the atomic contrast in such AFM images, we propose a simple model in which the Pauli repulsion is assumed to follow a power law as a function of the probed charge density. As the exponent in this power law is found to be largely independent of the sample molecule, our model provides a general method for simulating atomically resolved AFM images of organic molecules. For a single perylene-tetracarboxylic-dianhydride (PTCDA) molecule imaged with a CO-terminated tip, we find excellent agreement with the experimental data. Our model eliminates the need to take into account the full tip and sample system and therefore reduces computational cost by three orders of magnitude.

3 approaches, which mostly employ first-principles calculations of the complete system consisting of the sample and the tip, end up with a very high computational cost. Here we propose a model that takes into account the specific probe atom or probe molecule, but reduces the computational cost by about three orders of magnitude. With this model, the Pauli repulsive contributions to the interaction energy, the force and the frequency shift are calculated directly from the charge density of the sample molecule. We show that our model achieves good quantitative agreement with experimental results.

Pauli repulsion of two hydrogen-like atoms
The Pauli exclusion principle states that no two electrons may occupy the same quantum state simultaneously [20]. When two electrons of the same spin state intersect in space, their orthogonalized wavefunctions take on larger slopes to keep their overlap zero. The increase in slope causes an increase in kinetic energy, which reflects the Pauli repulsion. To understand this, we first consider two hydrogen-like atoms with the same electron spin and an atomic distance R. Their unperturbed wavefunctions ψ 1 and ψ 2 in spherical coordinates are where r is the radial distance, θ is the azimuthal angle, φ is the polar angle, Z n is the nuclear charge and atomic units are used. The orthogonal wavefunction ψ 2 can be analytically calculated using the Gram-Schmidt orthogonalization: As wavefunction ψ 1 is unchanged, the increase in kinetic energy E kin is solely due to the increase in kinetic energy of wavefunction ψ 2 : withT being the kinetic energy operator. For large atomic distances R, the increase in kinetic energy is where the exponential term which is equal to the charge density dominates. In a system with two different atoms, equation (5) becomes a function of both charge densities and some additional polynomial terms of the distance. Approximately, the increase in kinetic energy still decreases exponentially with the distance R, but less rapidly than the charge density.
For a system of two atoms with different nuclear charges Z n 1 and Z n 2 , as is typically the case in experiments, the increase in kinetic energy for large distances R becomes a much more 4 complicated relation: The increase in kinetic energy is a mixture of exponential functions and polynomials which cannot be simplified. However, for a system of two atoms with different nuclear charges, the increase in kinetic energy also decreases mostly exponentially for larger distances R. Therefore, we propose a simple power law for the increase in kinetic energy: It solely depends on the sample density ρ s and the two parameters A and B, with the exponent B being smaller than 1. With this we separated the sample and probe atoms with the probe atom entering the two parameters. In the following, we will show that this power law is a good model for the increase in kinetic energy which stems from the overlap of the sample and probe wavefunctions. Furthermore, we determine the two parameters for a couple of simple atomic systems and a system of two molecules: a perylene-tetracarboxylic-dianhydride (PTCDA) sample molecule and a CO tip molecule. This molecular system is then compared with experimental images.

Pauli repulsion of two atoms
Next, we examine the increase in kinetic energy for a simple system of two isolated atoms: the first atom is the sample atom and the second is the probe. The wavefunctions were calculated employing density functional theory (DFT) [21] and the highly optimized code CPMD 2 . We applied the Perdew-Burke-Ernzerhof exchange-correlation functional [22] and ab-initio normconserving pseudopotentials [23]. The wavefunctions were expanded into plane waves [24] with an energy cutoff of 300 Ry. First, the wavefunctions of the sample and probe atom were computed separately. For the atomic species in which the p-states are not completely filled, e.g. the C atom, the three p-states were equally filled, e.g. for the C atom with 2/3, 2/3 and 2/3. In this way, the atoms are spherically symmetric. The wavefunctions of the sample and the probe atom were then orthogonalized to each other, and the increase in kinetic energy was calculated. Note that the actual increase in kinetic energy of a diatomic system would be lowered because the orthogonalized wavefunctions could change to obtain a lower total energy. The kinetic energy would then decrease, and other energy contributions such as the exchange correlation energy would increase in the same order of magnitude [19]. Here, we neglect this for simplicity, and later see that even with this assumption the calculated images agree very well with the experimental ones.
We first consider a C atom as the sample and an O atom as the probe. In figure 1, the increase in kinetic energy E kin from the overlap of the two atomic wavefunctions is shown as a function of the distance of the two atoms. For comparison, the charge density ρ s of the sample C atom as a function of the radial distance r is also shown. Then we fit our model E model kin from equation (7) to the calculated values for the increase in kinetic energy E kin using the sample density ρ s . The fit was performed for distances R ranging from 2 to 5 Å. This interval includes the range of distances from 3.0 to 4.5 Å, where atomic-resolution images can be obtained as we estimated from recent experimental works [7]. Furthermore, the differences of the logarithm of the increase in kinetic energies and the logarithm of the fit function were minimized. This takes into account the exponential behavior. We obtain the fit parameters A = 1324 eV and B = 0.79. Because our model is a power law, the fit parameters depend on the unit of the charge density. In the following, all charge densities are in 1/Å 3 and all energies in eV. The resulting model for the increase in kinetic energy is shown in figure 1. When comparing directly calculated energies E kin with the model energies E model kin , they agree excellently over an energy range of almost three orders of magnitude for distances longer than 2 Å.
The same fit procedure for the model was executed for four different sample and probe atoms: H, C, N and O. These four atom types are those most commonly found in organic molecules. The results for the fit parameters A and B are given in table 1. Within one row, i.e. for a fixed probe atom, parameter A increases with the nuclear number. This increase reflects the fact that the Pauli repulsion, e.g. the overlap of the sample and probe wavefunctions, increases with the nuclear number. In contrast, the exponent B decreases with the nuclear number and is smaller than 1. The Pauli repulsion, as it stems from exponentially decaying wavefunctions,

Pauli repulsion of two molecules
In the same way as for the atoms, we calculate the increase in kinetic energy of two approaching molecules: PTCDA, which is the sample molecule, and CO, which is the probe molecule. As in previous work [7,19,25], we include only the PTCDA molecule and the CO molecule in our calculations and neglect both the substrate and the metallic part of the tip. We performed calculations with a metal tip behind the CO molecule and found that the Pauli repulsion remained unchanged. If the sample molecule interacts only weakly with the substrate, as in the case of NaCl, the presence of the substrate adds just a constant attractive background to the tip-sample interaction and can be neglected [7,19]. For chemisorbed systems it would be computationally too expensive to calculate the images directly. This should be possible with our simple model: it could be applied to the total charge density calculated for the molecule on the substrate. A further approximation in our simple model is that the CO molecule is fixed and oriented perpendicularly to the plane of the PTCDA molecule, with the oxygen atom pointing towards the PTCDA. Therefore, our model cannot account for any relaxations. Atomic relaxations, especially bending of the CO molecule at the tip, causes lateral distortions of the images [26,27]. Calculations of these relaxations are not included. However, calculated images can be easily compared to the measured distorted images which show the molecules laterally expanded. The lateral distortions of the experimental images due to the CO bending can be reduced with a simple method [27]. At the smallest experimentally accessible tip-sample distances (at approximately Z = 3 Å), the effects of the CO bending become increasingly important. To understand all the details of the contrast formation in this regime the CO bending has to be taken into account in the calculations. However, this renders calculations of entire molecular images computationally very expensive and this is beyond the scope of this paper. Furthermore, our simple model neglects electrostatic forces. Neglecting electrostatic forces is justified as it was demonstrated that they have no significant influence on the appearance of atomically resolved AFM frequency shift images of a molecule with a non-vanishing quadrupole moment [12].   The charge density of the sample PTCDA molecule is shown in figure 2 for three different characteristic positions above the molecule: above the C atom (x = 1.43 Å, y = 0 Å) close to the center of the PTCDA, above the corner O atom and above the center of the PTCDA. The charge density above the C atom always remains the largest for all distances, whereas the charge density above the O atom has a large value for small distances and decays faster. It is even lower than the charge density for the center site for distances larger than 2.15 Å. The calculated increase in kinetic energy E kin and the increase in kinetic energy from the model E model kin are shown as functions of the molecular distance Z . The molecular distance Z is defined as the distance between the PTCDA plane and the O atom of the CO molecule. As the PTCDA consists mostly of C atoms, we perform the fit for one position, above the C atom, to obtain the parameter A and the exponent B of the model. We consider the same range of distances as for the atoms.   figure 4 for a more quantitative comparison. The density is plotted again for the two perpendicular distances z = 2.46 Å and z = 3.73 Å. The directly calculated increase in kinetic energy and our model are plotted for the two corresponding molecular distances Z = 2.46 Å and Z = 3.73 Å. The profile lines are through the center y = 0 Å and through the outer O atom y = 2.24 Å. When changing the distance, there is a difference of two orders of magnitude in the density and one order of magnitude in the increase in kinetic energy. Still our model fits quite well with the directly calculated increase in kinetic energy. However, as the model is a simple power law of the density, the model inherits the curve progression of the density. This leads to a noticeable difference between the model and the directly calculated increase in kinetic energy. The model seems to be more corrugated which can also be seen in the images of figure 3. The reason is that the model does not take into account the size of the probe when scanning laterally. In principle, better agreement of the images could be obtained by convoluting the model image with the probe, which would lead to an averaging out of the model image. Note that although the charge density exhibits most of the features visible in the images of the increase in kinetic energy, it cannot be compared quantitatively.

Calculation of the frequency shift from the increase in kinetic energy
In the last step, we calculate the frequency shift using our model to compare it directly with the experiment. Just like in classical mechanics, in quantum mechanics the virial theorem relates the kinetic energy with the potential energy. In quantum mechanics the Pauli repulsion is a non-conservative force with a non-homogeneous potential. Therefore, the relation between the kinetic energy and the potential energy becomes more complicated [28,29]. This relation becomes accessible for a diatomic system because there is only one degree of freedom [29,30]. For our system which is a dimolecule system and also exhibits only one degree of freedom, the Pauli repulsion energy can also be easily approximated from the increase in kinetic energy: Inserting our model from equation (7) into equation (8) and taking the second derivative, we obtain the following for the frequency shift: The frequency shift can therefore be easily computed directly from the charge density of the sample. The result is shown in figure 5 together with the experimental frequency shift. As our model only considers the Pauli repulsion and omits the attractive van der Waals and electrostatic interactions, the calculated frequency shift is always positive. In comparison, the experimental frequency shift is always negative and exhibits a halo from the attractive van der Waals interactions. In the region inside the halo where repulsive contributions dominate, the experimental image compares qualitatively quite well with the calculated image: the five C rings are visible, whereas the O atoms are not. The π-electrons of the C atoms are further extended compared to O atoms which can be seen in the charge density. Furthermore, both the experimental and calculated frequency shifts exhibit a similar corrugation of approximately 4 Hz. Taking all into account, the calculated frequency shift compares qualitatively and quantitatively quite well with experiments. The attractive interactions could be added to the model using an empirical potential, but would excessively complicate our model, or could possibly be removed from the experimental images.

Summary
The calculation of the image is performed in four steps: determining the sample density ρ s from a single calculation of the sample molecule, calculating the increase in kinetic energy for one lateral position and several perpendicular distances Z , fitting the two parameters A and B, and then calculating the whole image. We applied the same procedure to a pentacene molecule and to graphene. B = 0.80 and then determining the parameter A with a single calculation considering the full tip-sample system will result in good images for other aromatic molecules and possibly also for other inert organic molecules. We showed the relations between charge density, increase in kinetic energy and the experimental measured frequency shift. We corroborated that the charge density compares qualitatively well to the measured images. To obtain a correct distance dependence and furthermore allow a quantitative comparison, we proposed a simple model based on a power law. This model is thought to be a fast and easy accessible tool for the interpretation of AFM images. With this, the number of calculations is reduced from a couple thousand to a few 12 self-consistent DFT energy calculations. Therefore, the total computational cost can be reduced by three orders of magnitude when calculating the frequency shift images.