Electrostatic forces on a CO molecule

(1)

Electrostatic forces on a CO molecule

Simulating an AFM image

Author:

Supervisor:

Date:

Alexandra Wadensjö

Magnus Paulsson

June

12, 2018

Course code: 2FY80E

(2)

Abstract

(3)

1 Introduction

In atomic force microscopy (AFM) a sharp tip is swept a few Ångströms above a substrate. The measured force between tip and substrate is then used to produce images of the substrate with subatomic resolution. What forces are the origin of AFM images? Is it sucient to assume that the main forces are electrostatic forces? Knowing the distribution of electrons in a sample, is it possible to produce a simulated image of the surface that correspond to the AFM image of the same sample? These are the questions that we will investigate in this thesis.

This year, Okabayashi et al. published an article in PNAS [1] in which they discuss the impact the tip has on the oscillation frequencies of a carbon monoxide molecule (CO) on a copper surface. To measure the inuence, they have used a combination of AFM and inelastic tunneling spectroscopy. From this study, we have received high quality AFM images at dierent tip-heights. By analyzing their result, we will try to answer the questions above.

As starting point we will use Maxwell's and Coulomb's laws and deduce the electro-static forces that are present in the AFM experiment. Then, we will compare these forces with the results in the Okabayashi et al. report. The forces that best coincide with the experimental result, will then be applied on the charge distribution obtained from density functional theory and compared with the AFM images.

1.1 AFM

AFM makes it possible to investigate surfaces and single molecules with subatomic res-olution. It was invented by Binnig in 1986, as a development of the scanning tunneling microscope (STM). STM had been introduced by Binnig, Rohrer, Gerber, and Weibel just a few years earlier in 1981 and in 1986 Binnig and Rohrer were rewarded with the Nobel prize [2]. While the STM only operates with conducting substrates, the AFM can produce images on both conducting and insulating substrates.

In STM a sharp tip is brought over the surface with a distance small enough (a few Ångström) such that a tunneling current can ow between tip and sample. With distances in this range, the interaction force between tip and sample is big enough to aect the level of tip. These level displacements can be measured and are used in AFM to translate the forces inicted on the tip by the surface to an image.

One way to get a force sensitive tip, is to use a thin cantilever which acts like a spring with a spring constant k. As a spring, the tip has an eigenfrequency f = 1

2Π

q

k

m. When

a surface force eects the cantilever, it causes a shift in the eigenfrequency ∆f. It is the measured change of ∆f that is transferred to an image. Consequently, the stiness of k is crucial. It needs to be soft such that it is sensitive to the tip-sample force, but it must also be suciently sti to avoid snap into the substrate, i.e. crashing the tip. For the same reason, the mass is also important [3].

(5)

1.2 CO molecule on a copper surface

To investigate our questions, we have the results in the report Vibrations of a molecule in an external force eld by Okabayashi et al. [1]. As comparison material we have received a calculated charge distribution from the thesis supervisor [4]. Using quantum chemistry, one can calculate the charge distribution of the electrons from a given molecule or crystal [5]. This method is applied on a sample constructed by Paulsson to imitate the sample used in the experiment. To simplify the calculations, the Cu(111)-crystal used in the experiments is exchanged to a Cu(100)-crystal to obtain cubic structure.

The data we have got at our disposal, is the geometry of a single carbon monoxide (CO) on a copper plate with an AFM tip above, also constructed of copper. Figure 1 illustrates the geometry of tip and sample. The CO molecule is placed orthogonally to the copper plate on top of a copper atom. Closest to the copper is the carbon atom with the oxygen atom on top. The bond length of the CO molecule is in our sample 1.29 Å.

Figure 1: Geometry of tip and sample. The CO molecule is the grey and red atoms in the center of the gure. The grey atom is carbon and the red atom is oxygen.

Together with all atom positions, we also know the charge distribution of the electrons. We will later apply our ndings from the investigations of the electrostatic forces on the charge distribution and see if we can simulate images that resembles the AFM images in the Okabayashi et al. report.

In addition to the sample data, we also have their published AFM images of the CO molecule as well as plots of the frequency shift and the force by the tip (see gure 2). The z-values in the plots are small since the z-axis in the plots has its origin at "point contact" where the tip is non-oscillating and the tip tunneling conductance would reach the conductance quantum G0 = 2e2/h ∼ (12.906Ω)−1 [1]. The exact height over the CO

(6)

al. report plots with our calculated forces and try to assess the distance between CO molecule and origin more accurately.

(a) (b)

(c) (d)

(7)

2 Simple electrostatic forces

Let us narrow the problem down and only consider a single tip copper atom and how it interacts with a single carbon monoxide (CO) dipole sample. What resolution can we expect with charges of this dimension and with reasonable interatomic distances?

From Coulomb's law, we know that if there are two static point charges, Q and q, the electric eld E caused by Q inicts a force F on q. Also, E is proportional to the charge Q and inversely proportional to the square of the distance, r2. That is

¯

F = q · ¯E = q · Q · ˆr 4πε0· r2

(2.1) where 0 is the vacuum permittivity.

Since the charges are static, Maxwell's laws give that this relation constitutes a con-servative eld - given by an electric potential ¯φ [6]. We can write

¯

E = −∇φ (2.2)

and

¯

F = −q∇φ (2.3)

Thus, if the potential is known we can calculate the force. In the following sections we will dissect what potentials that can be of interest and the electric elds from each potential. Then, we can combine the electric elds with several types of charges and see what forces the combinations give. Having summed up these forces, we can analyze which combination of forces - what net force - that best coincide with the experimental forces in the Okabayashi et al. report.

2.1 Interaction potential between tip and sample

Starting point for this analysis is the potential for a point charge q φpointcharge=

q 4πε0r

(2.4) For the dipole we have two charges +q and −q with the distance d from each other. The potential from the two charges is equal to the sum of the potential for each charge.

φdipole= φ+q+ φ−q = +q 4πε0r+ + −q 4πε0r− = q 4πε0 r−− r+ r+· r− (2.5) When dening the potential for the dipole, it is customary to consider the potential at a distance r much greater than d (see gure 3a). With such a geometry, one can approximate the distance r+ as

r+ ≈ r −

dcosθ

2 (2.6)

Similarly, the distance r− is approximated as

r− ≈ r +

dcosθ

2 (2.7)

Then the expression r−−r+

(8)

q− tip q+ θ ¯ d ¯r− ¯ r ¯ r+

(a) Dipole far from tip

q− tip q+ θ ¯ d ¯ r + ¯d ¯ r

(b) Dipole close to tip

Figure 3: Dipole with distance d between the charges and distance r between dipole and tip. The vector ¯d is here equal to dˆz.

Now, we can put this into eq[2.5] and rewrite the product qdcosθ as the dot product of the dipole moment ¯p and the unit vector ˆr.

φdipole =

¯ p · ˆr 4πε0r2

(2.9) Here, one can argue that r is not much greater than d and that in the frameworks of AFM, d and r are in the same magnitude. In that case, we cannot make approximations as above, but need to keep regard the potential as a compound potential of two point charges. The distances to the tip is then ¯r+ = ¯r − ¯d/2 and ¯r− = ¯r + ¯d/2and we get the

potential φtwocharges = q 4πε0 1 |¯r − ¯d/2| − 1 |¯r + ¯d/2| (2.10) Above the CO molecule where r = z we have

φtwocharges = q 4πε0 d z2_{− d}2_/4 (2.11) To test the two expressions, we put in typical values for the CO molecule obtained from DFT. We let q = 0.2e and d for a CO molecule is 1.29 Å. A plot with these two curves (see gure 4) shows that there is a small dierence between the two potential expressions when positioned straight above, or close to straight above, the dipole. For large x the two curves coalesce. That is - when the distance is increased this dierence decreases, just as expected. At a distance of 2 Å, the dierence at x = 0 is 0.11 V where the compound potential is 1.04 V and the approximated potential is 0.93 V. When the distance is increased to 4 Å, we have that the compound potential is 0.2385 V and the approximated potential is 0.2323 V. Now the dierence is 6.2 mV. We assess these as signicant dierences. Since this is the range that we are interested in, we will further on use the compound potential φtwocharges and treat the dipole as two point charges.

(9)

(a) z = 2Å (b) z = 4Å

Figure 4: Potential in V from a dipole with charge q = 0.2e and length d =1.29Å. Blue is the approximated potential and red is the compound potential. When x >> d the potentials coalesce.

equal to ¯r and ¯r− be equal to ¯r + ¯d. This makes r+ = r and r−= |¯r + ¯d|. We have

φtwocharges = q 4πε0 1 r − 1 |¯r + ¯d| (2.12) If we expand |¯r + ¯d| we get

|¯r + ¯d| = |rˆr + d(cosθˆr + sinθ ˆθ)| = |(r + dcosθ)ˆr + dsinθ ˆθ| (2.13) If the tip is positioned right above the dipole, r− has only a contribution in the

ˆz-direction. That is, r−= z + d which is what we expected. Then, the potential is

φtwocharges = q 4πε0 d z2_{+ zd} (2.14)

2.2 Contributions to the electric eld from the sample

From Maxwell's equations we have that the denition of the electric eld is ¯

E = −(∇φ) (2.15)

For a single point charge, there are no contribution to the potential in the ˆθ- or ˆφ-direction and the electric eld is thus

¯

Epointcharge = −(∇φpointcharge) =

qˆr 4πε0r2

(2.16) When deriving the potential for the dipole (eq[2.12]) we have contributions in both ˆz and ˆθ directions. To simplify, we can split the derivate in two parts. We let the expression of the electric eld stay in this separated form.

(10)

2.3 Theoretical forces between tip and sample

What forces can we assume being parts of the net force? In our simplied example with only our dipole and the tip atom, the contributions to the force is the tip charge and the electric eld from the dipole. Possible electric elds have been examined in the previous section. Let us now consider the charge of the tip. It can have either a point charge, a permanent dipole moment and/or an induced dipole moment. Together there are six distinguished contributions to the net force we are interested in:

1) a force between two point charges Fqq;

2) a force between a dipole sample and a point charge tip Fdq;

3) a force between a point charge sample and a permanent dipole tip Fqd;

4) a force between a dipole sample and a permanent dipole tip Fdd;

5) a force between a point charge sample and an induced dipole tip Fqi;

6) a force between a dipole sample and an induced dipole tip Fdi.

To be able to deduce which of the forces that are of interest, we will go through the distance relation of each force. Later, we will use the distance relations and compare with the power laws that best coincide with the experimental forces. Since we are interested in the distance straight above the CO molecule, we will denote the distance with z.

Fqq and Fdq are simply the products between the electric eld from the sample and

the tip point charge.

¯

Fqq = ¯Epointcharge· qtip (2.18)

¯

Fdq = ¯Edipole· qtip (2.19)

Fqq is decreasing with factor 1/z2 while Fdq is decreasing with the distance with the

order 1/z3_.

For a permanent dipole tip, we have a dipole moment in the tip - a vector. To be able to get the eect of the electric eld on the dipole moment we must take the gradient of the electric eld. Then we have a dot product between two vectors.

¯

Fqd= ∇ ¯Epointcharge· ¯pperm (2.20)

¯

Fdd = ∇ ¯Edipole· ¯pperm (2.21)

The same applies to an induced dipole. Here we assume that ¯pind is proportional to

the electric eld from the sample with a constant α, hC2_m

N i ¯ Fqi = ∇ ¯Epointcharge· ¯pind (2.22) ¯ Fdi = ∇ ¯Edipole· ¯pind (2.23) ¯ pind= α ¯Esample (2.24)

Fqdis, just as Fdq proportional to 1/z3. Fdd is proportional to 1/z4 due to the ∇ ¯Edipole

term. Fqi is proportional to 1/z5because of the ¯pindterm, which holds a 1/z2term. Scalar

multiplied with ∇ ¯Epointcharge this gives an 1/z5-factor. In the same way is the ¯pind term

in Fdi proportional to 1/z3. Scalar multiplied with ∇ ¯Edipole, this gives an 1/z7-factor.

(11)

2.4 Experimental forces between tip and sample

Let us now consider what forces we can deduce from the plots in the Okabayashi et al. report (see gures 2a and 2b). We can state at once that the forces are attractive, and we can also see that there is a signicant background force that probably comes from induced dipoles of the copper plates. Thus, we concentrate on the plot where the background has been subtracted, gure 2b.

From the force plot we can measure the FWHM and see that it alters depending on the height. The curves are narrower at lower heights and widens as the distance gets greater. However, the increase of FWHM seems to stabilize at z > 1.5 Å. After this point the curves coincide at a distance of about 3-4 Å from x = 0. Another dierence is that at z > 1.5 Å the curves seem to have only one extreme point and do not reach positive values.

We estimate the FWHM to about 3.6 Å at z = 170 pm. From the force plot we have also extracted the magnitude of the force at dierent heights (see table 5a). As we have seen in the introduction, the z-values in the plots are low and should be extended by 1-2 Å. In table 5a) we have added these extended values. The range we are interested in is thus from 2.3-4.1 Å

A plot of this force is presented in the Okabayashi et al. report [1], and we have obtained the data from the author. The plot in gure 5b present the extracted gures, with 1Å added to the z-values. In the plot we can see that at about z = 2.5Å the curve's behavior is changed. It is rst at z > 2.5Å that the curve starts to behave as reversed proportional to the distance with some exponent. That resembles what we have seen with the FWHM, which also shifts behavior at approximately the same height. This can perhaps be explained by chemical bondings.[1] At such small distances, the electrons might overlap between tip and sample. Due to the shift of behavior, we concentrate on the ratio from z > 2.5Å.

z [pm] 130 150 170 190 210 z+1Å[pm] 230 250 270 290 310 z+2Å[pm] 330 350 370 390 410 F [pN] -135 -80 -55 -45 -40

(a) (b)

Figure 5: F as a function of the distance z. (a) Table of values of the force F in pN at dierent heights, z. (b) Plot of |F| [pN] as a function of z + 1Å [Å]. Red curve is the extracted values from the force plot in the Okabayashi et al. report [1], and the black curve is from the ratio plot in the same report. The behavior of the black curve shifts at about z = 2.5Å. This can perhaps be explained as chemical bondings. The function is reversed to simplify the reading.

If the force is caused by simple electrostatic interaction, the net force will have some exponent to z, Fnet ∼ 1/z?. As we have seen when summarizing possible theoretical

(12)

k4/z5 and k5/z7. To nd which of them that has the closest resemblance with the report

values, we plot them together, see gure 6.

Figure 6: k1/z2 (green), k2/z3 (cyan), k3/z4 (magenta), k4/z5 (blue) and k5/z7 (red)

plotted with the experimental values (black), z-axis in [Å]. Between 2.5 < z < 2.8Å the experimental values have the closest resemblance to k4/z5 or k5/z7. At z > 2.8Å we have

that k1/z2 and k2/z3 are most similar to the experimental curve.

In the plot the distance to origin is chosen to be 1 Å. We nd that the exponential decrease of the experimental values has no resemblance with any individual inverse power law. If we focus on z > 2.8 the experimental curve has the closest resemblance with k1/z2

or k2/z3. Where z < 2.8, we nd that k4/z5 or k5/z7 follow the experimental curve the

most. The curve k3/z4 do not have any resemblance with the experimental values. When

increasing the distance to origin to 2 Å, we still do not have any resemblance with the experimental curve. However, we still see resemblance with k5/z7 at the lower z-values

and for the higher z-values we have that k2/z3 is overlapping the experimental curve.

It seems like we will have to use a combination of forces, to be able to explain the interaction between tip and sample.

2.5 Comparing experimental forces with theoretical forces.

When comparing our theoretical forces with the experimental forces, we will begin by examine what forces that coincide the most with the experimental values. Then we will compare the FWHM of the forces with what we have been able to deduce from the experimental plots.

To translate what we just have found out by examine the experimental values, we will not have any use for the force proportional to k3/z4. That is, the force between two

permanent dipoles Fdd. The other forces are still interesting. To be able to compare the

forces with the experimental values, we must use representative quantities for charges, dipole moments and distances. We use qpointcharge = 0.2e, qdipole = 0.13e and d = 1.18Å

which seems like a reasonable magnitude estimation from DFT. Since we do not know the charges of the tip, we have tested dierent values with the magnitude of the sample as starting point and tted the charges for each force in order to let them coincide with the experimental data, see gure 7. Fqq had a tip point charge of −0.13e, Fdq had a tip

(13)

a positive ˆz direction. We will at this point leave the forces between the sample and an induced dipole moment out, since they seem to only be important for the lower z-values. Also, when plotting the calculated forces with F as a function of z, we have found that the most reasonable extension to origin in the experimental data is about 1Å. With greater distance to origin we need to use unrealistically large charges on the tip. We will further on assume that 1Å is an accurate position of the origin in the experimental values.

Figure 7: The experimental data together with the forces Fqq (green), Fdq (blue) and Fqd

(red) in [pN]. z-axis in [m].

With these chosen and tted values, we get forces in the same magnitude as the experimental forces. We see that among the forces, Fdq and Fqq follows the curve the

most even though they are not as steep at the lower z-values. Even though one would assume that Fdq and Fqd should overlap, we see that there is a signicant dierence

between the curves. This dierence is due to the treatment of the sample dipole as two point charges. Thus, when d gets much smaller than z, the dierence between Fdq and

Fqd decreases.

Fdq is slightly steeper than Fqq, and Fqq has a tendency of converging earlier than

Fdq. Fqd is converging at too high z-values with our initiated sample charge and is not

interesting for further analyzing. That means that we can now deduce that there is no permanent dipole moment in the tip.

Now we will try to make the force curves steeper at the lower z-values. We add the fast decreasing induced tip dipole Fqi to Fqq (Fsum1) and Fdi to Fdq (Fsum2). This makes the

functions a bit steeper at the lower z-values while it has the slower decreasing behavior of Fqq and Fdq at the higher z-values. By reducing the point charge of the Fqq-function

to −0.12e and set the α-factor to 2 ∗ 10−40_C2_m/N1 _{we get a good t to the experimental}

values, see gure 8a. Regarding Fdq we altered the point charge to −0.25e and set the

α-factor to 10−40C2m/N, see gure 8b. This did not cause such substantial change for the Fdq function, even though it follows the bend in the experimental curve to a higher

degree. It also made it possible to lower the charge of the tip to a more reasonable value. For Fqq, on the other hand, this summation made a signicant impact and the function

now follows the experimental curve at z < 2.8 Å better without losing to much at the higher z-values.

The summed forces still do not follow the experimental curve at the smallest z-values

1_{2 ∗ 10}−40_C2_m/N _{corresponds to 7.8 ∗ 10}−5_e2_{Å/pN. That means that an induced dipole ∼ 0.1 eÅ at}

(14)

(z < 2.6) Å, which might be the lower limit where we can explain the interaction sample-tip with electrostatic forces without chemical bonding.

(a) Fsum = Fqq+ Fqi in green. (b) Fsum= Fdq+ Fdi in blue.

Figure 8: F [pN] as a function of the distance z [m]. The experimental data plotted with two dierent force sums.

We also want to compare the forces' FWHM. Since we are not interested in interactions with a permanent dipole tip, we focus on the forces to a point charge tip and compare them with forces to a point charge tip with an induced dipole moment. In gure 9 the tip charge is swept at constant height along the x-axis while z is increased from 2.5 - 3.5 Å by 0.2 Å per step. We choose to start at this height since it seems like that there are a limit at 2.5 Å and that we cannot explain the experimental values at lower z-values.

We let all sample values be the same as earlier. When plotting the single forces Fqqand

Fdq the tip charge is −0.13e and −0.3e, respectively. For the plots of the summated forces

Fsum1 and Fsum2 the tip point charge is −0.12e and −0.25e, respectively. The α-factor is

set to 2 ∗ 10−40_C2_m/N for F

sum1 and to 10−40C2m/N for Fsum2. When assessing the form

of each plot, we see that both Fdq and Fsum2 has a tendency of taking positive values far

from origin, which the experimental curves do not do. Another resemblance to note, is that the dierent force curves of Fqq and Fsum1 coincide between x = 3 − 4 Å and have

then the same force amplitude independent of the height z, just like the experimental curves. The force curves of Fdq and Fsum2 coincide closer to origin and are then dispersed

again.

At z=2.7 Å we can measure FWHM to: Fqq ∼ 4.13 Å; Fdq ∼ 3.16 Å; Fsum1 ∼

3.84 Å; Fsum2 ∼ 3.13 Å. Compared with the experimental values, where we estimated

FWHM to about 3.6 Å, the closest resemblance is Fsum1. Neither Fqq, Fdq nor Fsum2 is

particularly close, even though they are not that far away. Altogether, the force that has most resemblance with the experimental values is the net force Fsum1. This indicates that

(15)

(a) Fqqin pN. Both tip and sample are treated

as point charges. FWHM is 4.13Å. (b) Fpoint charge. FWHM is 3.16Å.dq in pN. Sample is a dipole and tip is a

(c) Fsum1= Fqq+ Fqi in pN. FWHM is 3.84Å. (d) Fsum2= Fdq+ Fdiin pN. FWHM is 3.13Å.

Figure 9: Force between tip and sample in pN. The tip is swept at constant height along the x-axis. z range from 2.5Å (green) - 3.5Å (black) by 0.2Å per step. FWHM is measured at z = 2.7Å.

2.6 Summary, electrostatic forces

When comparing experimental data with calculated forces, we have found that there is not a single force that can explain the experimental ndings. Instead we must combine dierent forces to compare with the experimental results. When plotting the experimental values next to dierent force combinations, we have obtained the best resemblance from the combination Fsum1 = Fqq+ Fqi and Fsum2 = Fdq + Fdi. That is, the tip behaves like

a point charge with an induced dipole moment, while the sample is either a point charge or a dipole.

We reach the same conclusions when evaluating the shape and the FWHM of the forces. We do not have any single force that is close to the experimental FWHM at about 3.6 Å, and neither is Fsum2. The best match to the FWHM comes from the combined

force Fsum1, which has a FWHM of 3.84 Å compared to 3.6 Å for the experiment. Fsum2

diers also from the experimental values and Fsum1 in the sense that it has a tendency of

reaching positive values at z > 2.5 Å, which neither Fsum1 nor the experimental plots do.

In the section describing the interaction potential we saw that a dipole must be seen as two point charges at these distances. Maybe that can explain why a force dominated by an interaction between two point charges best coincide with the experimental plots?

(16)

moment.

(17)

3 Finite dierence method

Charges act upon each other with forces due to the electric potential that each charge is surrounded with, as seen in previous section. However, in this problem there are many charges close to each other making it complex to determine the potential in each point. A solution is considering the problem not as single charges, but as a charge distribution. When the charge distribution is known, one can calculate the electric eld using one of Maxwell's equations

∇ · ¯E = ρ/0 (3.1)

We know that the electric eld E is the negative of the gradient of the potential φ. Put this in eq[3.1] and we obtain Poisson's equation

− ∇2_{φ = ρ/}

0 (3.2)

If one knows the continuous charge distribution, this partial dierential equation can be solved analytically. Here we have the charge distribution on a discrete grid. We choose to solve the problem with the nite dierence method where the derivatives are approximated with nite dierences. We will in the following sections explain how this method works and how it can be applied when solving Poisson's equation.

3.1 Finite dierence method in 1D

Considering a function A in one dimension. The derivative of A(x) is dened as dA

dx = limh→0

A(x + h) − A(x − h)

2h (3.3)

Instead of letting the distance approaching zero, we let the distance be set to a small distance dx. Then, we have a nite set of x-values xnand the derivative around the point

xi can be approximated by the value of A in xi−1 and xi+1. This is illustrated in gure

10. A(x) x dx Ai−1 Ai Ai+1

Figure 10: Finite dierence method in one dimension. Using the same argument for the second derivative, we have

(18)

which can be simplied to d2A

dx2 ≈

A(x + dx) − 2A(x) + A(x − dx)

dx2 =

= 1

dx2(Ai−1− 2Ai+ Ai+1)

(3.5)

3.2 Finite dierence method in 3D

Transferring the arguments in previous section to three dimensions, the point Aijk is

surrounded with points in ˆi-, ˆj- and ˆk directions. An illustration of this in two dimensions is found in gure 11. y x dx dy _A 11 A21 A31 A12 A22 A32 A13 A23 A33 A14 A24 A34

Figure 11: Finite dierence method in two dimensions. An approximation of the second derivative in three dimensions gives

∇2_{A = δ}2_{xA + δ}2_{yA + δ}2_{zA ≈} 1

dx2(Ai+1,j,k − 2Ai,j,k+ Ai−1,j,k)+

+ 1

dy2(Ai,j+1,k− 2Ai,j,k+ Ai,j−1,k)+

+ 1

dz2(Ai,j,k+1− 2Ai,j,k + Ai,j,k−1)

(3.6)

If the distances dx, dy and dz are chosen to be equal, this expression can be simplied further. In this problem however, the sample dimensions are unsymmetrical, and the small distances are dierently chosen. We let the expression stay in this form.

3.3 Solving a partial dierential problem with the nite dierence

method

In previous sections, we have seen that we want to solve the partial dierential equation −∇2_{φ = ρ/}

0. In Mathematica, this problem can be solved by writing φ and ρ as vectors

(19)

and ρ is written in the same manner.

As stated in eq[3.5], we can approximate the operation ∇2_φ _{in one dimension as}

∇2_{φ =} d2φ d2_x ≈

1

dx2(φi−1− 2φi + φi+1). With φ written as a vector as in eq[3.7], we write

the Laplace operator ∇2 as a nearly diagonal (N × N) -matrix L

L = 1 d2_x                −2 1 0 . . . 0 1 −2 1 0 . . . 0 0 1 −2 1 0 . . . 0 ... ... ... 0 . . . 0 1 −2 1 0 . . . 0 ... ... ... 0 . . . 0 1 −2 1 0 0 . . . 0 1 −2 1 0 . . . 0 1 −2                (3.8)

With this type of construction, the potential is implicitly assumed to obey the bound-ary condition φ0 = φN + 1 = 0.

Now we can write the expression ∇2_φ _as

∇2_{φ = Lφ =} 1 d2_x                −2 1 0 . . . 0 1 −2 1 0 . . . 0 0 1 −2 1 0 . . . 0 ... ... ... 0 . . . 0 1 −2 1 0 . . . 0 ... ... ... 0 . . . 0 1 −2 1 0 0 . . . 0 1 −2 1 0 . . . 0 1 −2                               φ1 φ2 φ3 ... φi ... φN −2 φN −1 φN                (3.9)

3.4 3D nite dierence method

Let us transfer the one-dimensional solution to a problem with a three-dimensional po-tential, Nx× Ny× Nz, where φ can be written as an (Nx× Ny× Nz)-matrix. We need to

rewrite φ = φijk as a vector. To be able to survey the position index, we let the matrix

φijk be unfolded to a vector in following manner

φ =                   φ111 φ112 ... φ11Nx φ121 ... φ1NyNx φ211 ... φNzNyNx                   (3.10)

(20)

To be able to perform the operation in eq[3.6] in a similar way as in one dimension, we need to extend the Laplace matrix to a (NxNyNz × NxNyNz) -matrix. Each row and

column are indexed as φ in eq[3.10]. Since we have chosen to keep the distances dx, dy and dz distinct, we need to implement this in the matrix. To be able to see the indexing scheme in the Laplace matrix we consider a problem where Nx = Ny = Nz = 2. That

is, L as an (8 × 8) -matrix. Also, in this three-dimensional case, we let the boundary conditions be φ0 = φN + 1 = 0.

With intent to illustrate this with a perspicuous matrix, the expressions 1

d2_x, _d12_y and

1

d2_z are replaced by respectively a, b and c and the sum of

1 d2_x+ 1 d2_y + 1 d2_z is replaced with d. L =             −2d a b 0 c 0 0 0 a −2d 0 b 0 c 0 0 b 0 −2d a 0 0 c 0 0 b a −2d 0 0 0 c c 0 0 0 −2d a b 0 0 c 0 0 a −2d 0 b 0 0 c 0 b 0 −2d a 0 0 0 c 0 b a −2d             (3.11)

Note how the o-diagonal alternates due to the boundary conditions. The problem ∇2_φ can now be written as

∇2_{φ = Lφ =}             −2d a b 0 c 0 0 0 a −2d 0 b 0 c 0 0 b 0 −2d a 0 0 c 0 0 b a −2d 0 0 0 c c 0 0 0 −2d a b 0 0 c 0 0 a −2d 0 b 0 0 c 0 b 0 −2d a 0 0 0 c 0 b a −2d                         φ111 φ112 φ121 φ122 φ211 φ212 φ221 φ222             (3.12)

(21)

4 Solving Poisson's equation

We have received a reconstructed electron charge distribution of a similar sample to that used in the Okabayashi et al. report [4]. To be able to simulate AFM images from this data, we want to convert it from charge distribution to a potential. That is, we want to solve Poisson's equation (eq.[3.2]) by applying the nite dierential method as we have seen in the previous sections. To be able to do this operation, we need to process the data. We also need to adjust the Laplace matrix to t our sample. In the following sections we will go through how we have prepared this data to be able to solve Poisson's equation and get the potential of the sample. Finally, we will solve Poisson's equation.

4.1 Electron charge distribution

At our disposal we have the electron charge distribution ρe of the CO molecule and the

closest copper plate atoms corresponding to the copper plate fracture in gure 1. It is a small sample with size 10.32 × 10.32 × 40.40Å and a resolution of 64 × 64 × 250 points. The sample size is chosen to keep the calculations on a manageable level. Already this resolution consumes a lot of the computer resources. We let (Nx, Ny, Nz) = (64, 64, 250)

and (dx, dy, dz) =10.32 N x , 10.32 N y , 40.40 N z .

The sample is not placed in empty space but has periodic boundary condition: x0 =

xNx, y0 = yNy and z0 = zNz. When reaching the boundary, one starts over on the opposite

side. ρe is composed in such way that the copper plate is cut o in the xy-plane eight

layers from the surface creating a gap between two planes. The gap is big enough to later on simulate a tip in between the layers without letting the tip be aected by the upper layer. The two layers are together suciently thick to behave like a metall. Figure 12 illustrates ρe in the yz-plane at the x-value where the CO-molecule is placed. The unit

for the charge distribution is [e/Å3_].

Figure 12: ρe in the zy-plane at the location of the CO molecule [e/Å3]. The CO molecule

is encircled with red.

(22)

4.2 Smearing of proton charges

ρe only includes the charges from the valence electrons and we need to add the positive

charges, ρp. Since the electron charges in ρe comprise only the valence electrons, we

should add the same amount of protons as each atom has valence electrons. For copper that corresponds to 11 protons, oxygen gets 6 protons, and we assign 4 protons to the carbon atom.

Apart from ρe, we have the positions for all atoms. These positions are more accurate

than our grid and the positions on the grid must be rounded. Having accurate positions for the atom cores is crucial. Small misplacements of the atom cores results in a dipole moment that causes a non-zero potential in the gap. We found that rounding to the closest integer gave a satisfying gap potential.

At rst, we tried to allocate the proton charges concentrated as point charges at these positions. It turned out that the numerical solution would not converge with such sharp distribution of charge. Another problem with concentrating the proton charge is that the nite dierence approximation is accurate only when we have a smooth function. To solve these problems, we smeared the proton charge over a few adjacent grid positions.

Statistically, a three-dimensional gaussian distribution seemed as an appropriate choice of smearing technique, with equal distributions in all directions. All atoms where given equal radii and the same standard deviations, regardless of their actual sizes. The attempt when setting the radius and standard deviation, was to overlap the electron distribution. A point where the overlap is important is just over the CO molecule, why there is good point to assess the result of the proton distribution. A radius of 7 grid points and a standard deviation of 2 grid points results in a covering which result in a gap with lowest point at -0.027e. 7 grid points corresponds to 1.13Å and 2 grid points corresponds to 0.32Å. The sum ρe+ ρp constitutes the total charge distribution ρtot. Figure 13 shows

ρtot, unit [e/Å3], at the CO molecule in the zy-plane.

(23)

4.3 Laplace matrix with periodic boundary conditions

The calculation method used to obtain the charge density uses periodic boundary condi-tions. Thus, the potential at the boundaries is not zero but we have periodic boundary conditions. We need to implement this in the Laplace matrix and let the potential at po-sitions x0, y0 and z0 be equal to xNx, yNy and zNz. Considering the end position, it starts

over from the beginning again and vice versa. The operation we perform on the diagonal φi,j,k is still _d−22_x +

−2 d2_y +

−2

d2_z, the dierence appears in the side positions. In the event of

an index exceeding its boundaries, it retakes the counting on the opposite extreme point. For index i that is i − 1 ≤ 0 ⇒ i − 1 = Nx− 1 and i + 1 ≥ Nx ⇒ i + 1 = 1. This also

applies for j and k.

To be able to see the periodicity, we let Nx = 3, Ny = 4 and Nz = 5. That is, the

Laplace matrix is a 60 × 60 -matrix. To increase the transparency, we plot a colored version of the matrix using the same indexing as before (see gure 14). Blue color is the sum −2 d2_x+ −2 d2_y + −2 d2_z. Yellow color is 1 d2_x, light orange is 1

d2_y and the darker orange is

1 d2_z.

Figure 14: Laplace matrix with periodic boundary conditions. Examples of periodic boundary conditions for x, y and z are encircled with red.

4.4 Solving Poisson's equation with Mathematica

To nd the potential for this charge distribution we need to solve Poisson's equation ∇2_{φ = −ρ/}

0. To simplify our calculations we have chosen to write this as a matrix

equation where the Laplace matrix L is the linear map of φ to ρ.

Lφ = −ρ/0 (4.1)

(24)

the equation. Then we have

φ = −L−1(ρ/0) (4.2)

Calculating the inverse of the Laplace matrix is not possible for a (M × M)-matrix if M is larger than ∼ 5000. In our case, the matrix is ∼ (1M × 1M). Instead we use the Mathematica function LinearSolve which can iteratively solve large sparse linear equations. With the input L and ρ it nds the missing φ in the matrix equation eq[4.1].

To solve the equation in this manner, we must rewrite the matrix ρi,j,k to a vector as

we have seen in previous sections. That is

ρ =                   ρ111 ρ112 ... ρ11Nx ρ121 ... ρ1NyNx ρ211 ... ρNzNyNx                   (4.3)

Writing ρ with this index entails that we need to make L a (NxNyNz × NxNyNz

)-matrix. With our sample values, that makes an (1024000 × 1024000)-)-matrix. That is a lot of operations to compute. However, we can lower the amount of operations. As we have seen when plotting the smaller (60 × 60)-Laplace matrix in previous section, the matrix mostly contains zeros. That means that we can use a sparse matrix which only hold the non-zero values and does not use memory for the empty cells. Then the computer only need to perform operations with the active cells.

Now, we can solve Poisson's equation for ρtot. The Mathematica code with the solution

is in the appendix. Figure 15 shows the potential φ at the x-position where the CO molecule is located. With a known potential we can calculate the electric eld and nally try to simulate AFM images.

(a) 3D plot of the potential in [V]. (b) Density plot of the potential in [V].

(25)

5 Simulated AFM images

In this section we will nally simulate a tip being swept over a sample surface and try to further deduce the nature of the sample-tip interaction with these results.

To be able to do this, we have in previous sections solved Poisson's equation for the potential. With a known potential, we can calculate the electric eld at chosen heights. These heights are then the constant heights at which we sweep our tip. Since we have AFM images of the sample at 3.5Å (see gure 2) we will also use that height. Then we can compare maximum force and width and shape of the CO molecule with the experimental images. We calculate the image at another height as well to have more data to compare with. Since 3.5Å is close to the maximum height of the experimental sample, we use the lowest possible tip height as the other height. As we have seen in previous sections there is a limit of about 2.6Å under which we cannot trust electrostatic forces to explain the tip-sample interaction. Therefore, we choose the closest following height - 2.7Å.

Before we concentrate on the tip, we will pay some attention to the electric eld at these heights, see gure 16. From these, we can deduce that we have a positive electric eld above the CO molecule while it is slightly negative over the copper plate. The electric eld has a maximum of 0.2 V/Å over the CO molecule at 3.5 Å and a maximum of 0.5 V/Å at 2.7 Å. Except for the negative forces over the CO molecule we can expect to have positive, though considerably lower, forces over the copper surface.

(a) z = 2.7Å (b) z = 3.5Å

Figure 16: Density plots of the electric eld in [V/Å].

In gure 17 we show how the electric eld is changing with the height z right above the CO molecule. It is plotted together with the power laws k1/z2, k2/z3 and k3/z4.

The electric eld follows the function k2/z3 well, except for smaller z-values. It is unclear

what potential the function corresponds to, but it is probably due to the numerous charges being placed close together. The electric eld from the dipole does not behave as its main component is from a single point charge. It is rather a compound of multiple charges. Not only the dipole charges, but it has also contributions from the surface charges. At long distances it seems like the impact from the surface charges is decreased, while the dipole contribution is more dominant.

(26)

Figure 17: The electric eld as a function of the height z (black) together with the exponential functions k1/z2 (blue), k2/z3 (red) and k3/z4 (magenta). The electric eld is

measured right above the CO molecule. The height is in Å.

start by initiating the point charge values we have found and from these values iterate to nd a reasonable charge for the tip. When we have reached a suitable charge, we will apply an induced dipole moment proportional to the electric eld. After some adjustments we should hopefully obtain an image that resembles the experimental images.

5.1 Finding realistic point charges of the tip

In the previous section Simple electrostatic forces, we concluded that we must expect the tip to be a negative point charge, and that it probably harbors an induced dipole moment. What was a bit uncertain was the magnitude of the point charge and the proportion factor α of the induced dipole moment. The point charge altered between -0.12e and -0.25e depending on whether the potential of the sample comes from a point charge or a dipole. In this section, we will try both point charges and discuss in what range the tip point charge is and we sweep over the surface at both 2.7Å and 3.5Å (see gures 18 and 19).

(a) qtip = −0.12e (b) qtip = −0.25e

Figure 18: AFM images at 3.5Å in [pN]. Tip is a point charge.

(27)

(a) qtip = −0.12e (b) qtip = −0.25e

Figure 19: AFM images at 2.7Å in [pN]. Tip is a point charge.

highly exceed the maximum from the experimental plots. Even at the higher height we exceed the maximum, even though the violation is not as big as at the lower heights. We need a point charge that gives forces lower than the experimental maximum since we also expect the tip to harbor an induced dipole moment. Using a smaller charge of -0.07e we obtain the correct force magnitude (see gure 20). At 3.5Å the force is slightly to lower than maximum and at 2.7Å the force level is quite good.

(a) z = 3.5Å (b) z = 2.7Å

Figure 20: AFM images in [pN]. Tip is a point charge of 0.07e.

(28)

5.2 Induced tip-dipole

We have seen in previous sections, that including an induced dipole moment proportional to the electric eld by some constant α seems important. We have previously estimated this constant to be in the magnitude of about 10−40_C2_m/N _{which corresponds to 3.90 ∗}

10−5e2Å/pN. As we can see in gure 21, this result in too large forces. It also results

in a slightly ellipse shaped form of the CO molecule and four positive areas around the molecule. Maybe, that is the four closest copper atoms beneath the CO molecule that is slightly negative charged due to the CO molecule.

Figure 21: AFM image in [pN], z = 3.5Å. Tip is a point charge of 0.07e and has an induced dipole moment proportional to the electric eld by a factor α = 10−40_.

(a) z = 3.5Å (b) z = 2.7Å

Figure 22: AFM images in [pN]. Tip is a point charge of 0.07e and has an induced dipole moment proportional to the electric eld by a factor α = 5 ∗ 10−42_.

When lowering α to 5 ∗ 10−42_C2_m/N _{(1.95 ∗ 10}−6_e2_{Å/pN), the CO molecule is again}

(29)

At 3.5 Å we know the approximate width of the CO molecule from the FWHM of the force plot in gure 2b. We estimate the width at F < -20 pN to be about 3.5 Å and for F < -10pN the width is about 5 Å. For our AFM image at the same height we measure the core (up to -20 pN) to 2.3 Å and for F < −10 pN the width is 3.9 Å. The measurements from the experimental force plot are unprecise, but we get an indication that the width of the CO molecule in our image is in approximately the right size. We also note that describing the tip as a point probably gives a lower FWHM than a broadened charge distribution would.

5.3 Discussion on the calculated forces

Now when the electric eld is known, and we have estimated realistic parameters for the tip, we can compare our results with the experimental values on the force. We will compare both the width of the curves as well as the plots of the force as a function of the height. Here, we will again compare the point charge tip with the dipole tip to conclude which option that is most realistic.

We start by comparing the width of the curves. In gure 23 the tip is swept along the x-axis. Here the tip has been assigned i) a point charge of -0.07 e (gure 23a); ii) a dipole moment of 0.07 eÅ (gure 23b); iii) an induced dipole moment with an α-factor of 5 ∗ 10−42_C2_m/N (gure 23c). For the point charged tip, the FWHM is 2.90 Å. The

dipole tip has an FWHM of 2.58 Å. The same FWHM applies for the induced dipole tip. All measurements of the FWHM is at z =3.7 Å. For the experimental forces, the FWHM at z = 270pm is much wider - approximately 4.4 Å.

All curves are a bit antisymmetric but are otherwise rather similar to the experimental force plots when z < 170pm. A noticeable dierence between the experimental plot and the calculated plots is, as we have already seen, that there are positive values only in the latter. Compared with the theoretical force plots, they have similar behavior. Here we note that while the theoretical force plots have nicely gathered curves some distance from origin, these curves have a common intersection but are then shattered. This tendency resembles of the forces from a dipole tip and is only seen in the experimental plots for z > 170pm. The point charge tip has however a slower increase of the width than the dipole tip, which makes it slightly more similar to the experimental plots.

When the tip has a point charge and an induced dipole moment we get the plot in gure 24. Here, the FWHM is 2.90 Å again. Still, we see that the width of the force is much narrower than for the experimental plots. In our calculations the tip point charge is concentrated in a point, which is not the case in the experiments where the charges are distributed over a bigger volume. This should result in more narrow curves; the question is how much narrower? If it only eects the width slightly, the dipole tip has still too narrow curves. If it eects the width to a high degree, the dipole tip might be more accurate than the point charge tip.

(30)

(a) Tip is a 0.07 e point charge. FWHM is

2.90 Å. (b) Tip has a 0.07 eÅ dipole moment.FWHM is 2.58 Å.

(c) Tip has an induced dipole moment, α = 5 ∗ 10−42C2m/N. FWHM is 2.58 Å.

Figure 23: Force on the tip due to the calculated sample potential in pN. The tip is swept at constant height along the x-axis. z range from 2.7 Å (green) - 3.7 Å (black) by 0.2 Å per step. FWHM is measured at z = 3.7 Å.

(31)

a factor of 1.6. When these parameters are adjusted, the curves coalesce to a high degree, see gure 26a. If we keep the distance to origin but use the dipole tip instead we need to multiply with a factor 1.5 to let the curves coalesce, see gure 26b. Both the point charge tip and the dipole tip have a rather good t to the experimental values. We note that the point charge tip coalesces with the curve down to z =2.9 Å. Below this height it diers from the experimental values. As we have seen, this is the lowest point below which we cannot explain the values with only electrostatics. The dipole tip diers from the experimental values at this point as well, but is also divergent at z > 3.3 Å. We conclude that the dipole tip does not t as well as the point charge tip.

(a) Height of origin is 0.8 Å. (b) Height of origin is 1.4 Å

Figure 25: Force [pN] as a function of the height z [m]. Experimental values (black) together with the calculated values (red).

(a) Point charge tip. Calculated values are

multiplied with a factor of 1.6. (b) Dipole tip. Calculated values are multi-plied with a factor of 1.5.

Figure 26: Force [pN] as a function of the height z [m]. Experimental values (black) together with the calculated values (red and blue). Height of origin is 1.4 Å. Both tips have an induced dipole moment.

When we tried to nd the proportions for the tip, we assumed that the height to origin was 1 Å. After adjusting the distance to origin to 1.4 Å and multiplying both the tip point charge and the proportion factor α with 1.6 we get the AFM image in gure 27a. The tip point charge is now -0.112e and α is 8 ∗ 10−42_C2_m/N_{. Both values are closer to the}

(32)

pN and 3.9 Å for F < −10 pN. We had that the experimental images have a width of about 3.5 Å for F < −20 pN, and 5 Å for F < −10 pN. Next, in gure 27b we have the AFM image of a dipole tip. The distance to origin is still 1.4 Å but here, we have the multiplied the dipole moment and α with 1.5. Now, the dipole moment is 0.105 eÅ and α is 7.5 ∗ 10−42C2_m/N_{. The dipole moment used for the theoretical forces was 0.24 eÅ.}

In this image we measure the width of the molecule at F < −20 pN to 1.77 Å and at F < −10 pN to 3.39 Å. An even sharper image!

(a) Point charge tip. (b) Dipole tip.

Figure 27: Final AFM images of the CO molecule. z = 3.9 Å.

As we have seen in this section, we have some parameters that we can adjust to obtain forces in the right range. By altering the charge of the tip, we can adjust the force amplitude to a suitable range. When we compare a point charge tip with a dipole tip we see changes in the image resolution. A dipole tip obtains a less wide molecule than a point charge tip, especially at the lower heights. We can also see this dierence when we plot the force as a function of the height. The dipole tip has a steeper function than the point charge tip.

Further, we can change the proportion factor α of the induced dipole moment. By making the factor larger we obtain more narrow force curves when we sweep the tip along the x-axis. It also results in a steeper curve when we plot the force as a function of the height. We have also seen that a too high proportion factor eects the shape of the CO molecule in the AFM image. Finally, we can alter the distance to origin. By changing the distance, we can adjust where on the z-axis the experimental plot is placed. This means that we get another relation between the force and the height while the shape of the relation is still the same.

Is it possible to deduce anything with so many parameters, or are our results too ambiguous? We have two xed points, which we can compare our results with: i) The experimental plots and AFM image; ii) The potential of the DFT charge distribution.

(33)

since that result in strange shapes of the molecule. It can be smaller, but we assess that we have reached a reasonable range. This magnitude on the induced dipole is satisfying and it results in force plots that resembles the experimental plots well.

What is harder to conclude is whether the tip is a dipole or a point charge. When plotting the distance relations of the electric eld we saw that it has an exponent in power law 1/z3 _{or 1/z}4_{. With a point charge tip, we keep this distance relation. For a dipole,}

(34)

6 Conclusion

With AFM images and data from experiments as starting point, we have been able to simulate AFM images using simple electrostatics and more advanced numerical calcula-tions. Closer than z ∼ 3 Å we cannot explain the interaction using electrostatics due to chemical bondings.

By plotting the ratio of the force at dierent heights we have been able to compare its behavior with theoretical forces. The theoretical forces have been obtained by using a simplied model with only the CO molecule and a tip consisting of only a point charge or a dipole. From this comparison, we obtained the most satisfying result with a point charge tip. Also, we deduced that at these short distances, compared with the inter-molecule distance of the CO molecule, we cannot approximate the potential of the molecule as a dipole with a d << r. In addition, we have also seen that to explain the shape of the experimental forces at dierent heights, we must assume that there is an induced dipole moment on the tip.

With a calculated potential obtained by the DFT charge distribution, we reached similar results. The calculated electric eld has a distance relation of the power laws 1/z3

or 1/z4_{. That implies that the force is decreasing slightly too quickly with a dipole tip,}

compared with the experimental forces. By using these properties for the tip and compare with the experimental AFM images, we estimate the tip to have a point charge of about -0.112e. Further, for the induced dipole moment we estimate a proportionality constant α of the order of 8 ∗ 10−42C2_m/N. The conductance quantum is reached at a distance to

origin of approximately 1.4 Å, why this distance has been added to the z-values in the experimental plots in our nal comparisons. With these parameters we have reproduced AFM images that resembles the experimental images.

(35)

References

[1] Peronio A. Paulsson M. Arai T. Okabayashi N. and F. J. Giessibl. Vibrations of a molecule in an external force eld. In: PNAS, vol.115, no.18 (2018).

[2] Quate C. F. Binnig G. and Ch. Gerber. Atomic Force Microscope. In: Physical Review Letters, Vol. 56, No. 9 (1986).

[3] Ph. Schahauser and S. Kümmel. Simulating atomic force microscope images with density functional theory: The role of nonclassical contribution to the force. In: Physical Review, vol. 94 (2016).

[4] M. Paulsson. Charge distribution obtained from density functional theory. 2018. [5] Artacho E. Gale J. D Garcia A. Junquera J. Ordejon P. Sanchez-Portal D. Soler

J. M. The SIESTA method for ab initio order-N materials simulation. In: Journal of Physics: Condensed Matter, vol. 14, No. 11 (2002).

[6] L. A. Engström. Elektromagnetisk fältteori: Från bärnsten till fältteori. 1st ed. Lund: Studentlitteratur, 2000. isbn: 978-91-44-01510-1.

List of Figures

1 Geometry of tip and sample. The CO molecule is the grey and red atoms in the center of the gure. The grey atom is carbon and the red atom is oxygen. . . 2 2 Plots from the report of Okabayashi et al. [1]. a) Constant height

fre-quency shift ∆f. Tip is swept along the x direction while changing the vertical distance z by 5 pm. b) Vertical force acquired by deconvoluting the frequency shift ∆f (x, z) distribution in the z direction. The long-range background force is subtracted. c) Frequency shift image of the CO molecule. d) Current image of the molecule. . . 3 3 Dipole with distance d between the charges and distance r between dipole

and tip. The vector ¯d is here equal to dˆz. . . 5 4 Potential in V from a dipole with charge q = 0.2e and length d =1.29Å.

Blue is the approximated potential and red is the compound potential. When x >> d the potentials coalesce. . . 6 5 F as a function of the distance z. (a) Table of values of the force F in pN

at dierent heights, z. (b) Plot of |F| [pN] as a function of z + 1Å [Å]. Red curve is the extracted values from the force plot in the Okabayashi et al. report [1], and the black curve is from the ratio plot in the same report. The behavior of the black curve shifts at about z = 2.5Å. This can perhaps be explained as chemical bondings. The function is reversed to simplify the reading. . . 8 6 k1/z2 (green), k2/z3 (cyan), k3/z4 (magenta), k4/z5 (blue) and k5/z7 (red)

plotted with the experimental values (black), z-axis in [Å]. Between 2.5 < z < 2.8Å the experimental values have the closest resemblance to k4/z5 or

k5/z7. At z > 2.8Å we have that k1/z2 and k2/z3 are most similar to the

experimental curve. . . 9 7 The experimental data together with the forces Fqq (green), Fdq (blue) and

(36)

8 F [pN] as a function of the distance z [m]. The experimental data plotted with two dierent force sums. . . 11 9 Force between tip and sample in pN. The tip is swept at constant height

along the x-axis. z range from 2.5Å (green) - 3.5Å (black) by 0.2Å per step. FWHM is measured at z = 2.7Å. . . 12 10 Finite dierence method in one dimension. . . 14 11 Finite dierence method in two dimensions. . . 15 12 ρe in the zy-plane at the location of the CO molecule [e/Å3]. The CO

molecule is encircled with red. . . 18 13 ρtot in the zy-plane at the location of the CO molecule [e/Å3]. . . 19

14 Laplace matrix with periodic boundary conditions. Examples of periodic boundary conditions for x, y and z are encircled with red. . . 20 15 The potential φ in the zy-plane at the location of the CO molecule. . . 21 16 Density plots of the electric eld in [V/Å]. . . 22 17 The electric eld as a function of the height z (black) together with the

exponential functions k1/z2 (blue), k2/z3 (red) and k3/z4 (magenta). The

electric eld is measured right above the CO molecule. The height is in Å. 23 18 AFM images at 3.5Å in [pN]. Tip is a point charge. . . 23 19 AFM images at 2.7Å in [pN]. Tip is a point charge. . . 24 20 AFM images in [pN]. Tip is a point charge of 0.07e. . . 24 21 AFM image in [pN], z = 3.5Å. Tip is a point charge of 0.07e and has

an induced dipole moment proportional to the electric eld by a factor α = 10−40. . . 25 22 AFM images in [pN]. Tip is a point charge of 0.07e and has an induced

dipole moment proportional to the electric eld by a factor α = 5 ∗ 10−42_{. . 25}

23 Force on the tip due to the calculated sample potential in pN. The tip is swept at constant height along the x-axis. z range from 2.7 Å (green) - 3.7 Å (black) by 0.2 Å per step. FWHM is measured at z = 3.7 Å. . . 27 24 Force [pN] on the tip due to the calculated sample potential. The tip is

swept at constant height along the x-axis. z range from 2.7 Å (green) - 3.7 Å (black) by 0.2 Å per step. FWHM is measured at z = 3.7 Å. Tip is a point charge of 0.07e and has an induced dipole moment proportional to the electric eld by a factor α = 5 ∗ 10−42_C2_m/N_{. . . 27}

25 Force [pN] as a function of the height z [m]. Experimental values (black) together with the calculated values (red). . . 28 26 Force [pN] as a function of the height z [m]. Experimental values (black)

(37)

A Mathematica code for solving Poisson's equation

* * *

(*Input*)

{Nx,Ny,Nz} = {64,64,250}; (*# gridpoints*)

rhoMatrix; (*charge distribution, with size [Nx, Ny, Nz]*) sampleSize = {10.32, 10.32, 40.40}; (*In Angstrom*)

{dx, dy, dz} = {sampleSize[[1]]/Nx, sampleSize[[2]]/Ny, sampleSize[[3]]/Nz}; (*Distances*)

(*Constants*)

epsilon0 = 8.854*10ˆ(-12); (*Vacuum permittivity [F/m]*) (*Rewriting rho from matrix to vector*)

For[ii = 1, ii <= Nx Ny Nz, ii++, rhoVektor[[ii]] =

rhoMatrix[[Mod[ii - 1, Nx] + 1, Mod[Floor[(ii - 1)/Nx], Ny] + 1, Floor[(ii - 1)/(Nx Ny)] + 1]]

];

(*Constructing the Laplace matrix*)

NN = Nx Ny Nz; (*Size of Laplace matrix*) (*Diagonal and o-diagonal*)

L = SparseArray[{Band[{Nx + 1, Nx}, {NN, NN - 1}, {Nx, Nx}] -> 0., Band[{Nx, Nx + 1}, {NN - 1, NN}, {Nx, Nx}] -> 0,

Band[{1, 1}] -> -2*(1/dxˆ2 + 1/dyˆ2 + 1/dzˆ2), Band[{1, 2}] -> 1/dxˆ2, Band[{2, 1}] -> 1/dxˆ2,

Band[{1, Nx + 1}] -> 1/dyˆ2, Band[{1, Nx Ny + 1}] -> 1/dzˆ2, Band[{Nx + 1, 1}] -> 1/dyˆ2, Band[{Nx Ny + 1, 1}] -> 1/dzˆ2, Band[{1, Nx}, {NN - Nx + 1, NN}, {Nx, Nx}] -> 1/dxˆ2,

Band[{Nx, 1}, {NN, NN - Nx + 1}, {Nx, Nx}] -> 1/dxˆ2, Band[{1, NN - Nx*Ny + 1}, {Nx*Ny, NN}] -> 1/dzˆ2, Band[{NN - Nx*Ny + 1, 1}, {NN, Nx*Ny}] -> 1/dzˆ2}, {NN, NN}];

(*Adding o-diagonal for one side of the diagonal, transpose for other side*) L2 = SparseArray[Flatten@ Join[ Table[ Band[{iz Nx Ny + (Ny - 1) Nx + 1 , iz Nx Ny + Ny Nx + 1 }, {1 + iz Nx Ny + (Ny - 1) Nx , 1 + iz Nx Ny + Ny Nx } + Nx - 1] -> -1/dyˆ2, {iz, 0, Nz - 2}], Table[

Band[{1 + Nx Ny iz, 1 + Nx Ny iz + Nx (Ny - 1)},

{1 + Nx Ny iz + Nx - 1, 1 + Nx Ny iz + Nx (Ny - 1) + Nx - 1} ] -> 1/dyˆ2, {iz, 0, Nz - 1}]

], {NN, NN}];

(38)

(*Periodic boundary conditions -> we must choose a position that is equal to zero and implement in both the charge distribution and in the Laplace matrix.

We choose a point in the gap between the two layers of copper atoms.*) indx = Nx*Ny*150;

L[[indx, ;;]] = 0; L[[indx, indx]] = 1; rhoVektor[[indx]] = 0;

(*Solve Poisson's equation, phi = the electric potential*)

Electrostatic forces on a CO molecule