Simulations for Kelvin Probe Force Microscopy

(1)

Simulations for Kelvin Probe Force Microscopy

Julian van Velzen

Supervisor: Rinke J. Wijngaarden

Bachelor project

Photo-conversion materials group

VU University Amsterdam

(2)

Abstract

Photo catalysts are evolving as a new method in order to split water into H2and O2, where hydrogen

can be used as a fuel. The group photo-conversion materials at the Vrije University is doing research on the properties of the photo catalyst NaTaO3. In order to understand its characteristics, the local

electric elds of NaTaO3are studied by the method of Kelvin Probe Force Microscopy. In order to

understand the measured force, this research provides a numerical analysis of the tip sample force and other relevant characteristics. The numerical analysis is conducted in the environment of the commercial software Matlab. A script has been deployed that utilizes the Finite element method to solve the Poisson equation, subsequently the Maxwell stress tensor is integrated over both the tip and sample's surface in order to determine the electrostatic force. Results have been obtained that have a relative error of around 1% for the KPFM geometry, which is good enough for the purpose of this project. By calculating the force on the sample by acknowledging the second law of Newton, the accuracy can be improved.

(3)

Introduction

Recent studies have shown the opportunities of photo-catalyst systems to split water into its sub-stituent parts, hydrogen and oxygen. While hydrogen fuel becomes increasingly attractive [15] due of its zero emission and because of the depletion of conventional fuels, a method is required to produce hydrogen. One such method is to split water using a photo-catalyst under radiation of either direct or articial light:

H2O2H2+O2.

In the search for such photo-catalyst, the group photo conversion materials at the free university in Amsterdam, is investigating the properties of NaTaO3. Although the splitting capabilities of

this catalyst seem promising, yielding in a splitting rate of 9.7 mmol/h under UV radiation [7], the characteristics of this catalyst need further investigation. In particular in the Amsterdam group, the local electric elds on the surface of NaTaO3 under UV irradiation were investigated with

Kelvin Probe Force Microscopy (KPFM). KPFM measures the electrostatic force between the tip and sample. In order to interpret this experimentally observed force, and in particular its distance dependence, comparison models are required. Although analytical expressions for the electrostatic force exists for simple geometries, the more complicated experimental geometries require numerical methods to obtain the electrostatic force and other relevant properties. In this project therefore, numerical simulations of the KPFM electrostatic tip-sample force are performed.

The numerical simulations are conducted in the environment of the commercial software package Matlab. Matlab provides a library for solving partial dierential equations (PDE) that is based upon the widespread method of nite elements. The nite element method (FEM) is a tool that can be used to discretize continuous equations with the premise that the solution of a PDE can be divided over an array of sub domains with homogeneous scalar values.

Below, after presenting the theoretical background, the simulation environment is discussed for calculating the electrostatic force and other relevant properties. This involves the identication and construction of the relevant geometry and boundary conditions. Subsequently dierent solution methods are presented and compared. And nally output of the method is presented for the electrostatic force and (self) capacitance for various geometries. This project provides a framework that can be further used for numerical calculations on various KPFM geometries.

(6)

Chapter 2

Theory

2.1 Kelvin Probe Force Microscopy

The numerical simulations in this research are done to better understand kelvin probe force mi-croscopy (KPFM). KPFM is a tool that is used to measure the local contact potential dierence (CPD) between a conducting tip and a sample. By measuring the CPD nano-scale electronic prop-erties of the sample can be mapped. It is therefore, and because of its high spatial resolution, extensively used [11]. In this section the basic principles of the kelvin probe force microscopy are discussed. Furthermore, to test the validity of numerical data to simulations done by this research, and attempt is been made to obtain analytical expressions for the electrostatic force and capacitance for simple geometries.

2.1.1 Principles of Kelvin probe force microscopy

Kelvin probe force microscopy determines local contact potential dierences. The contact potential dierences are related to the work function, which is a fundamental property of materials. It corresponds to the energy dierence between the energy level of an electron in a solid, and the vacuum level directly outside the solid surface [12]. The expression for the work function is given by:

W =Evac− µ

e (2.1)

where W is the work function, e is the electron charge and µ is the electrochemical potential or Fermi level of the electrons. When two objects are characterized by dierent work functions, and an electrically connection is established, an electric current will ow briey from the object with lower work function (weak binding) to the object with higher work function (strong binding) until the Fermi levels become aligned. The objects thus get charged and a contact potential dierence (CPD) of

Vcpd=

Wtip− Wsample

e (2.2)

is established. The Fermi levels of the two objects, which we took as tip and sample, are now aligned, although the vacuum energy levels dier. The CPD is very sensitive to changes in surface conditions, such as contamination and temperature [6]. It is therefore very useful for understanding

(7)

physical and chemical properties of materials. Due to any potential dierence between tip and sample, an electrostatic force acts between the two given by:

Fes= 1 2 dC dzV 2_. _(2.3)

Due to Vcpd this force is generally non zero. In the KPFM measurements this force is nullied

by the introduction of an external applied bias voltage Vdc. This external voltage is set equal in

magnitude to Vcpd, but with opposite sign. Thus the total electrostatic force is nullied. To create

a very sensitive measurement, an extra AC potential Vacis applied between the tip and sample,

causing the tip to oscillate: the applied Vdcand Vaccosωtresults in the electrostatic force of:

F (t) = 1 2

dC

dz(Vdc+ Vcpd+ Vaccosωt)

2_. _(2.4)

This force can be divided into three spectral components:

F (t) = Fdc+ Fωcosωt + F2ωcos2ωt (2.5) where: Fdc= dC dz 1 2(Vdc− Vcpd) 2₊1 4V 2 ac (2.6) Fω= dC dz(Vdc− Vcpd)Vac (2.7) F2ω = − 1 4 dC dzV 2 ac. (2.8)

Fω (equation 2.7) corresponds to the oscillation of the cantilever and is directly proportional to

the dierence between Vdc and Vcpd. A lock-in amplier is used to detect the oscillation of the

cantilever. Vcpd is equal to the known value of the applied Vdcwhen the lock-in ampliers output

equals zero. Because Vcpd equals the dierence in work function between the tip and sample, the

work function of the sample can be arrived when the work function of the tip is known.

KPFM has two operational modes; Amplitude Modulation (AM) and Frequency Modulation (FM). In the AM mode, the amplitude of the oscillation ω of the cantilever is measured directly by the lock-in amplier as we just discussed. Alternatively, in FM mode the frequency shift, which is proportional to the force gradient and therefore to (Vdc− Vcpd)Vac, is measured. In FM mode, the

DC current is then applied to nullify the frequency shift, and therefore again determines Vcpd. The

KPFM - FM mode has a higher spatial resolution. Although the KPFM method achieves a high spatial resolution, potential artifacts may occur. These include:

• non-perfect zeroing of the lock in signal • electric cross-talk eects in the lock in signal • extra charge on the tip or sample

• dielectric layers on the tip or sample [1]

To be able to understand the eects of such artifacts and thus to better interpret experimental results, we want to be able to calculate the tip sample force, and in particular Fω and ∂C_∂z for the

situations just discussed. As a rst step we now continue with analytical approximations that are used to compare the numerical data.

(8)

2.2 Calculating the Force analytically for simple geometries

Since the simulations are conducted solely on static congurations that have no time dependency, all the relations can be described by the theory of electrostatics. Two of Maxwell's equations, that provide the framework of all the relevant physics, and therefrom derived quantities are discussed in this section. The quantities concerned are the electric potential, the capacitance, and the maxwell stress tensor.

2.2.1 Electric potential

In the simulations done by this research, an electric potential dierence is maintained between two electrodes, in particular the KPFM tip and the sample. From the rst and second Maxwell equations a useful expression for the potential is obtained as follows. Starting with the rst equation, for a charge distribution ρ, the electric eld E is given by:

∇ · E = ρ 0

. (2.9)

Mathematically it follows from the second Maxwell equation that, ∇× E = 0, so that we may write

E = −∇V. (2.10)

Combining equation 2.9 and 2.10 one obtains the Poisson equation: ∇2_{V =} −ρ

0

. (2.11)

This equation forms the central part in this research and provides the backbone of the Finite element method, which is described in section 2.3. When no charge is enclosed, the Poisson equations simplies to the Laplace equation: ∇2_{V = 0}_.

2.2.2 Capacitance

The capacitance C is the ratio of stored charge Q, and the potential (dierence) V: C = Q

V. (2.12)

For numerical calculations, Q can easily be derived by the use of Gauss law if the electric eld is known in all space. The expression for Q then reads:

Q = ˛

0E · ˆ~ ndS (2.13)

In both equations, Q is the charge that is transported due to the creation of the potential (dif-ference). In a system of multiple conducting charges, the charge on any conductor will aect the potential of all other conductors in the system [13]. This relation is only dependent on the geomet-rical conguration and the capacity. In a system of two charges, the charge on a conductor can be written in terms of the potentials of neighboring conductors. The solutions are:

(9)

Q2= c12V1+ c22V2.

Here the quantities c11 and c22 are called the self capacitance. They are dened as the charge to

potential ratio on the rst and second conductor respectively, when the other conductor is grounded. The quantity c12and c21are called the mutual capacitance and can be dened as the induced charge

to potential ratio on the rst and second conductor respectively due to the potential of the other conductor. From symmetry it follows that cij = cji. When equation 2.14 is written in the form of

2.12, it follows that the total capacitance equals: C = c11c22− c

2 12

c11+ c22+ 2c12

2.2.3 Maxwell stress tensor

In order to nd the force in numerical simulations, one can use the method of the Maxwell stress tensor. In the static case, which we are discussing here, the Maxwell stress tensor gives the force per unit area. It can be derived from the Lorentz force. We will not give its derivation here but simply give its denition as:

σij ≡ 0(EiEj− 1 2δijE 2_{) +} 1 µ0 (BiBj− 1 2δijB 2_). _(2.15)

Were E and B is the electric and magnetic eld respectively, and δij is the Kronecker delta. Since

we are only interested in static electric congurations, equation (2.15) can be simplied to: σij= 0(EiEj−

1 2δijE

2_). _(2.16)

Further more, on a conductor the electric eld is always perpendicular to its surface. Equation 2.16 is therefore only non zero if i = j. The tensor can then be simplied to σij = ₂0E2. The

electrostatic force on a surface is obtained by integrating the stress tensor over the surface D. The integral is given by:

~ F = 0 2 ˆ D E2nds.ˆ (2.17)

Using the maxwell stress tensor to calculate the force gives the advantage that the electric eld just has to be calculated around the desired surface, instead of in all space. This eases and speeds up the computation. However, because the gradient of the force is greatest near the surface, it is expected that this will easily result in an error in the numerical calculations. Furthermore, calculating the force by its local contributions requires only one one computation of the electric eld, this in contrast with taking the force as the derivative of the systems energy, which requires calculating the force at two distances, hence requires calculating the electrostatic problem at two distances.

2.2.4 Analytical approximations

In order to validate our numerical simulations, analytical expressions are highly desirable. For the sphere plane geometry the method of images can be used to nd exact expressions for the capacitances. These expressions however, include innite sums and are not very practical to use. For

(10)

more complicated geometries no exact expressions are known, and one must turn to approximations. One such approximation is proposed by Hudlet [5] and was shown to be very accurate.

In Hudlets approximation the geometry concerned can be summarized by gure (2.1). It con-cerns a axial symmetric conductor, that represents the tip in the KPFM measurements, and a innite plane surface, that represents the sample. A voltage dierence of V is maintained.

Figure 2.1: Tip sample characteristics [5]

The electrostatic force between the tip and sample can be calculated by the Maxwell stress tensor integral, which is described in equation 2.17. From this equation, the vertical component is taken: Fz= 0 2 ˆ D E2ds ˆn · ˆz (2.18)

In order to nd the electric eld everywhere on the surface of the tip, Hudlet postulates that the electric eld on each innitesimal surface of the tip, is what would be created by the parallel plate capacitance of two innite plates that are in the same orientation as the innitesimal surface on the tip and the surface of the sample. The electric eld is then given by

E = − V

(11)

where l(M) is the length of the eld line from the sample to the innitesimal surface at the tip. This method is applied on both the apex and the cone, as dened in gure (2.1). The derivation of the expressions for the force are not included in this paper, and can be found in the paper by Hudlet [5]. The derived expressions for the apex, and cone contribution to the force are respectively:

Fapex= π0R2 1 − sinθ0 z [z + R(1 − sinθ0)] V2 (2.20) Fcone= π0V2 [ln tan(θ0/2)]2 lnza zb + zc 1 za − 1 zb . (2.21)

However, when these two contributions are combined by Hudlet, there seems to be an error in this paper [5]. The correct expression for the force is:

Ftotal= π0V2 _R2 z [z + R]+ k 2 − lnz + R H − 1 + R/sinθ0 z + R (2.22) where k2₌_{1/ ln tan(}θ0 2)

2_{. The expression for the capacitance can be arrived from these equations,} starting with the derivative of the energy of the system

F (z) = ∂C ∂z

V2

2 (2.23)

which shows that 2.22 yields ∂C

∂z to compare directly with KPFM experiments and which can be

integrated to: C = 2 V2 ˆ F (z)dz = 2π0 R ln(1 +R z) − 1 k2 _−R sin θ0 + R + z lnH + z R + z − H + R (2.24) The accuracy of the approximation is important for this research, since it is used as a comparison, and test for the numerical results. Hudlet states a good agreement between his approximation of the sphere plane geometry, and the exact one. For small lift height (z/r < 0.01) the error is less than 1%, and for large distances (z/r > 4) the error is 5% at most. These errors are acceptable for the aim of this research, however, it does put limitations on the verication of the accuracy of the algorithm.

2.3 Calculating the force numerically

2.3.1 Finite element method

In order to simulate the force and capacitance in KPFM measurements, computer models are used. Because computers can only handle a nite number of numbers, it is not possible to solve the dierential equations directly. Instead, space must be discretized.

One method to do that is the nite element method (FEM). FEM is a method to numerically solve boundary value problems. The strategy that is used is to divide the solution space into a large number of sub-domains. The concerning quantities are then approximated to have uniform value over the smaller sub-domains. These smaller sub-domains have a discrete value and can be calculated by a computer. The complete solution is then formed as a superposition of all sub-domains. In order to discretize and solve the equations, the following steps have to be addressed:

(12)

• Division of the solution space into small sub-domains • Derivation of linear equations for the sub-domains • Solving the equations for the sub-domains

2.3.2 Division into sub-domains

The set of sub-domains that forms the total region of interest is called the computational mesh. In cylindrical geometry, that is used in this research, the solution properties depend on the r and z direction, and are independent of the θ-coordinate. The sub-domain volumes are toroidal objects with arbitrary projections in the r-z plane, that extents in a rotation about the θ axis. The shape of the projection on the r-z plane can be chosen freely, however triangles are used because of their simple shape and because every domain can be decomposed into triangles.

Because of the approximation that the electric potential is uniform over the sub-domain volume, it is expected that an increase in the amount of triangles will result in a better solution. Therefore a good mesh has a high density of triangles in regions where the potential is rapidly changing.

2.3.3 Poisson equation in cylindrical coordinates

The dierential equation in particular that is solved with the nite element method in this research, is the Poisson equation (2.11). This equation can be written in general form as

∇ · (c(r)∇u(r)) = f (r). (2.25)

The electrostatic poison equation is obtained when c(r) = (r), u(r) = V (r) and f(r) = −ρ(r) are substituted. In our case cylindrical symmetry is used in order to reduce the 3-D problem to a 2-D problem in the r-z plane. The Laplacian must therefore be rewritten in cylindrical coordinates as ∇2_{f =} 1 r ∂f ∂r+ ( ∂2f ∂r2 + ∂2f

∂z2)where f is independent of θ because of the symmetry. This Laplacian can be obtained from equation 2.25 by choosing c(r) = r.

f = ∇ · (r∇u) = (∇r · ∇u) + r∇2u = ∂ ∂ru + r∇ 2_u

this results in the poison equation in cylindrical coordinates, 1 r ∂u ∂r + ( ∂2 ∂r2u + ∂2 ∂z2u) = f r (2.26)

This derivation shows that the 3-D situation can be simplied by transforming the poison equation in the cylindrical equation. We will now continue with the discretization of this dierential equation.

(13)

2.3.4 Derivation of linear equations

The FEM will solve the partial dierential equation by converting it in a discrete form of KU = F that can be solved by conventional methods. In this form K is the stiness matrix, F is the load vector and U is the solution vector. This conversion is done by choosing a nite number of N trial functions φj(r, z)and assembling the total solution from these trial functions.

To discretise the Poisson equation, it must rst written in its weak form. Transforming the equation in the weak form has the advantage that it now only includes rst derivatives of u and test functions v, whereas the strong formulation contains second derivatives. Both sides of equation (2.25) are multiplied by a test function v and integrated,

¨

(∇ · (c∇u))vdA = ¨

f vdA. (2.27)

Equation 2.27 holds for all possible v. This equation is then integrated by parts using Greens theorem to result in the weak form of Poisson's equation:

¨ c∇u · ∇vdA − ˆ ˆ n · (c∇u)vdS = ¨ f vdA (2.28)

which includes only rst derivatives.

Before this equation is discretized, boundary conditions can be applied to simplify equation 2.28. In Matlab two types of boundary conditions can be applied. Dirichlet boundary conditions are of the form:

hu = cst, (2.29)

where the direct value of u is known. The value of u can be imposed by setting the value of test function v to zero at the boundary points, and requiring u in equation 2.28 to attain the desired value.

Neumann boundary conditions are conditions of the form:

c∇u + qu = g. (2.30)

In this case the derivative of u is dened. The next step is to derive a discrete approximation to the weak form. This is done by choosing the solution u(r, z) as a linear combination of N trial functions φj(r, z)with coecients Uj, u(r, z) = N X j=1 Ujφj(r, z). (2.31)

In order to solve equation 2.28, with N dierent unknowns for Uj, you need N dierent equations

with instances of v. Therefore, v in Equation 2.28 is replaced with v = Vi, which we choose

subsequently to equal φi. Substituting equation 2.31 in equation 2.28 then yields:

ˆ c(r, z)∇φi(r, z) · N X j=1 Uj∇φjdrdz = ˆ f (r, z)φi(r, z)drdz (2.32)

which is our nal discretized Poisson equation. Equation 2.32 contains N equations, in which Uj

are the unknowns which can be determined. By dening Kij =

´

c(r, z)∇φi(r, z) · ∇φj(r, z)drdz

and Fi=

´

φj(r, z)f (r, z)drdz equation 2.32 can be written as a linear system of equations:

(14)

Where U is the solution vector with the nodal values uj. The starting point of the Poisson equation

(2.25) is now dened by the test functions φj for all the the sub domains.

2.3.5 Solving the equations for the sub domains

The test functions φj are polynomials that are dened on each triangle in the mesh with value 1

at its node, and 0 at other nodes. For the use of this paper, only linear test functions are used, meaning that φj = aj+ bjx + cjywhere a, b and c are unknowns that have a value for each triangle.

Because a triangular mesh is used with three nodes per triangle, this results in three equations and three unknowns per triangle that can be solved. Matlab will iteratively determine the value of the test functions at each triangle and determine their contribution to the matrices K and F. The linear system 2.33 now only depends on the test functions φj that are dened by the nodal points

(15)

Chapter 3

Method

In this chapter the method to obtain and calculate the relevant properties of Kelvin probe force measurement, in particular the electrostatic force and capacitance, are discussed. The platform that has been used for programming is Matlab, although other programs such as COMSOL and Mathematica have been considered. Matlab has been chosen because of its wide use in the scientic world, and its easy applicability for further simulations. This chapter will discuss the deployment of the script, and in particular the construction and properties of the geometry, the applied boundary conditions, and the post processing and calculations of the electrostatic force and capacitance. The source code of the script can be found in the appendix.

3.1 Geometry

The geometry concerned is that of the Kelvin probe force microscope that includes a spherical apex, cone, cantilever and sample. In order to clearly visualize the properties of this geometry, the dimensions have been altered for use in this paper. The dimensions of the geometry, used in the simulation are summarized in table 3.1. Although in real KPFM the radius of the spherical tip is around 20 nano meters, for this experiment a larger scale is chosen because it speeds up the calculation. All dimensions are therefore multiplied by the factor 'scale'. Also, in older Matlab versions there seems to be an error in the mesher function when creating meshes for geometries with a scale smaller than 10−6_{, but this is solved in the 2013 version of Matlab by specifying}

Mesherversion to be0_R2013a0_.

For the KPFM calculations the same method has been applied to two geometries, gure 3.1a and gure 3.1b. In both cases the vertical axis corresponds to the z-axis and represents the symmetry axis. The horizontal axis corresponds to the radial axis. The whole geometry is rotated around this z-axis and gure 3.1 should therefore be interpreted as a cross section. The sample (indicated with number 2 in gure 3.1a and the x-axis in gure 3.1b ) is therefore a circular disk, and the tip has cylindrical symmetry. In order to test the validity of the method, the capacitance and force are tested also for the sphere plane geometry. This geometrys has the same characteristics as gure 3.1 but with just a sphere and sample.

The Matlab function which was custom build within this project, geomBuilder takes all the arguments specied in gure 3.1 and creates the decomposed geometry matrix. This matrix is used by Matlab to contain all the information about the geometry and can be constructed in several ways.

(16)

The easiest way would be to use the Matlab PDEtool, which provides a graphical user interface to construct the geometry, set the boundary conditions, and solve for the potential. However, since we want to vary the characteristics of both the geometry as boundary conditions in simulations, the decomposed geometry matrix is constructed manually. The value of all dimensions described in picture (3.1) can be given as arguments to the function geomBuilder winch then returns the decomposed geometry matrix.

property value*scale R 20 w1 30000 w2 0.5/scale h1 10000 h2 20000 h3 1000 scale 10^-6

Table 3.1: dimensions of the KPFM geometry

(a) (b)

Figure 3.1: KPFM geometry

3.2 Boundary conditions

Because the geometry describes a nite size area, boundary conditions need to be dened at the borders. Matlab allows Neumann and Dirichlet boundary conditions to be specied on the bor-der (see section 2.3). The sample and tip are maintained at constant potential dierence by the specication of the Dirichlet boundary condition. The value of these potentials can be freely given as arguments to the function scalarBoundaryBuilder, which was custom build within this project. The function scalarBoundaryBuilder, only accepts scalar values for Neumann or Dirichlet boundary

(17)

conditions. On the symmetry axis Neumann boundary conditions of zero ux are applied, which follows from the azimuthal symmetry. The boundary conditions on the outside border are harder to dene. The only condition that is certain, is that the potential will go to zero for r, z → ∞, however since we are dealing with nite size meshes, an approximation has been used. Ideally, no boundary condition would be applied on the outside border, since the system is dened suciently by the boundary conditions on the tip and sample. Matlab's PDE solver however, requires bound-ary conditions to be specied on all outside borders. On the outside borders dierent boundbound-ary conditions have been tested. It can seem reasonable to apply Neumann conditions of zero ux, suggesting that the potential is constant at the border. Another option would be to set Dirichlet boundary conditions on the border, and specify a custom value. This has been done by treating the tip as a point charge, and analytically determine the Dirichlet value at the border. The boundary condition is then specied by:

r0∗ Vt

pr2_{+ (z − p)}2 (3.1)

where Vtis the potential of the tip, r0is the radius of a constructed sphere that has the same volume

as the tip, and p is the z-coordinate of the center of mass of the tip in the z direction. This is reasonable since the tip behaves as a point charge at distances where r r0. Since in the simulations

we vary the lift height, the center of mass of the tip changes every iteration, and accordingly the boundary conditions have to be modied. Equation (3.1) can be given as an argument, as well as the value of the potential for the tip and the sample, to the function scalarBoundaryBuilder to specify the Dirichlet boundary condition at the border. The function scalarBoundaryBuilder will then return the boundary condition matrix.

3.3 Mesh and Solution

In order to compute the solution for the potential, the nite element method divides the domain into small sub domains called the computational mesh. The Matlab function initMesh is been used for creating the mesh. The function accepts the decomposed geometry matrix and returns the mesh data that consists of a list of points p, edges, e and triangles t.

To rene the mesh, there are two options. One possibility would be to repeatedly divide each triangle into four. This way triangle density would be homogeneously distributed over the domain. Although this method is very fast, expected is that a denser triangle distribution in the area of interest, that is around the tip where eld lines are relatively dense, results in a better solution. This can be obtained by adaptive mesh renement function from Matlab, adaptMesh. This function solves the PDE problem each generation of mesh renement, where it then calculates an error estimation, and renes the mesh in areas with the highest estimated error. This method is expected to achieve a more eective mesh that results in a closer approximation, however, it requires a lot more computing force because is has to go trough the cycle of solving the PDE problem and calculating the error estimation each time.

With the mesh data, decomposed geometry matrix, and the Boundary Condition matrix the PDE problem can be solved. The equation that is solved is:

−∇ · (c∇u) + au = f (3.2)

(18)

transforms to the cylindrical electrostatic Laplace equation: 1 r ∂ ∂r r∂V ∂r +∂ 2_V ∂z2 = 0. (3.3)

For solving this equation the Matlab provided function pdenonlin has been used, although custom functions have been tested. The function accepts the decomposed geometry matrix, the Boundary Condition matrix, mesh data and PDE coecients, and returns the solution u. The function pdenonlinsolves for the potential by the nite element method (section 2.3).

3.4 Further calculations

After the nite element solution u for the potential is obtained further calculations can be done. In particular we are interested in the electric eld E, (mutual) capacitance C and force F .

3.4.1 Capacitance

We are interested in particular in the mutual capacitance between the tip and sample, and not particularly in the self capacitance of an individual object.

The tip has a capacitance versus the grounded sample, but it also has a self capacitance versus the grounded innite sphere surrounding it. Both capacitances are in parallel and hence have

Q

V = Cself + Cts where Ctsis the tip sample capacitance we want to determine. The capacitance

that is calculated by C = Q

∆V, where Q is the charge on one object and ∆V is the potential dierence

between the tip and sample and corresponds to the sum of the self, and mutual capacitance. Hence in order to obtain the mutual capacitance, the self capacitances must be subtracted from the total capacitance. This is done by simulating an individual object with a potential dierence with respect to innity. The sample is xed at V = 0, hence the self capacitance of the sample does not contribute, and only the self capacitance of the tip has to be determined.

In order to reduce systematical errors, the same geometrical dimensions and mesh are used for the calculations for the self capacitance and the total capacitance. There is no easy way to generate a mesh that can be used for dierent geometries and hence the following method is used. First instead of considering the sample as a region outside the region of interest that is dened by its boundary conditions, the area of the sample is considered as a separate minimal region (area 2 in gure 3.1a) inside the solution space of the PDE problem (area 1 + 2). In this area, separate PDE coecients can be dened, and solved for equation 3.3. The electric permittivity of a perfect conductor corresponds to = ∞ so high values for are applied for this region in the PDE problem in order to simulate a conductor. To bring the sample to a certain potential the left side of the area 1, that lies on the symmetry axis and contains a boundary condition, is applied a Dirichlet boundary condition. Since the sample is a conductor with innite the potential will homogeneous distribute over the sample.

Since the sample is now part of the PDE solution space, a mesh will be created in both areas. We can now calculate the self capacitance by solving a dierent partial dierential equation that uses the same geometry and mesh. We set = 1 in the sample, and the left side of the sample that contains a boundary condition to a Neumann condition of zero ux, in order so simulate a vacuum. The tip is now in absence of the sample, and the PDE problem can be solved for an individual object with a potential in respect with innity. Calculating C = Q

(19)

For both calculations, of the total capacitance and self capacitance, the relation C = Q

V applies.

Since the potential is xed, we only have to determine the accumulated charge. As discussed in section 2.2.2, we can use the Gauss law,

Q = ˆ

~

E · ˆndS (3.4)

and integrate over a surface that encloses all the charge.

3.4.2 Force

The electrostatic force can be calculated with the Maxwell stress tensor integral (section 2.2.3). In order to do so, we must nd the electric eld, and integrate it over the desired surface. However, since the tip contains complex geometry where the solution is expected to deviate from the real solution, the force is calculated both over the surface of the tip and over the surface of the sample. Due to the second law of Newton, Fij = −Fji, the force on the tip is equal in magnitude, but

opposite in sign to the force acting on the sample. The integral that has to be evaluated is: ~

F = 2

ˆ

E2ndSˆ (3.5)

For both the capacitance and the force, the electric eld has to be determined rst.

3.4.3 Electric eld

The electric eld E is minus the gradient of the potential u. Since we only need the electric eld at the surface of the tip to calculate the Maxwell stress tensor integral, a method has been developed to determine the electric eld from the scalar values of u for specic triangles by the use of the Green Gauss theorem. The discrete form of the potential yields the electric eld as follows,

E = −(∇φ)t= 1 V X f φfA~f (3.6)

where φ is the potential at the centroid of each triangle. In 3.6 the sum runs over all the faces of the triangle, where φf is the value of the potential at the triangles face, Af is the vector describing

the face, that follows from the two nodes whereat it is connected, and V the area of the triangle. Notwithstanding that the build-in Matlab function for taking the derivative uses the same princi-ple, it is better optimized and therefore runs a bit faster. Our simulations will therefore use the Matlab function pdegrad. However, because both methods approximate the value of the potential as constant in the triangle, an additional error might occur.

3.4.4 Integration

In order to compute the electrostatic force and capacitance, the integral of equations 3.5 as 3.4 have to be evaluated. For the force, we want to integrate the electric eld over the surface of either the tip or the sample. For the charge we can use the fact that Gauss law (3.4) applies to any path enclosing the charge, so we can choose a path that lies away from the surface, where the gradient of the electric eld is lower and therefore the solution is expected to have greater accuracy. The dierent integration paths that we will therefore use are:

(20)

• at the surface of the tip • at the surface of the sample

for the calculations of both the force and charge, and additionally for the charge • at a square path around the tip

• at a square path around the sample.

Because the geometry involves cylindrical symmetry the area element ˆndS for an arbitrary surface that is symmetric in the θ direction is

ˆ

ndS = 2πr (drˆz + dz ˆr) (3.7)

where r is the r-coordinate of the triangles face midpoint, and dr and dz are the lengths of the face. Because of symmetry the radial contribution to the force is zero, and the force is on the z-direction only. Substituting equation 3.7 in equations 3.5 and 3.4 yields the summation for the charge and force respectively: Q = 2πX S r (drˆz + dz ˆr) · ~E (3.8) F = πX S E2rdr (3.9)

(21)

Chapter 4

Results

To investigate the electrostatic force and capacitance of complex geometries, rst the validity of our method has to be conrmed. Because no error estimation is available for the KPFM geometry, the method will be tested for dierent situations and geometries. Numerous simulations have been performed for simpler boundary value problems whereas numerical data was compared to data derived from known analytic expressions. Those simulations involve the following; rstly, the electrostatic force and capacitance were calculated for two concentric spheres and a sphere plane geometry. For these geometries analytical expression are available and can be compared. For the two concentric spheres a relative accuracy 4F

F of up to 10

−6 _{was reached, depending upon the}

mesh density. This determination uses the same method of integration as all other geometries. However, the geometry of two concentric spheres is far simpler than the KPFM geometry, and the high accuracy is therefore no proof that the method will be successful in other cases as well. Secondly, an adaptive mesh rener has been used in order to increase mesh eectiveness. In gure 4.1 both the capacitance of the sphere plane geometry (blue line) as the computational time is plotted against the number of adaptive mesh renement generations. Notwithstanding that the adaptive mesh rener is eective when the area of interest is small compared to the solution space, adaptive meshing will create an ineective mesh for the geometry 3.1a that is described in section 3.4.1. This is because a high mesh density is created over the whole sample, although we are in particular interested in the region around the tip. Because the adaptive mesh renement has to go trough the cycle of solving the equation and determine the error function every time, this runs a lot slower, as can be seen in gure 4.1. Therefore for geometry 3.1a a iterative mesh rener is used, and for 3.1b a adaptive mesh rener is used.

(22)

Figure 4.1: Mesh renement

4.1 Force

The force has been determined for both the sphere as the KPFM geometry by integrating the maxwell stress tensor (Equation 3.5) over a surface. By use of newtons second law Fij = −Fji the

integration can be evaluated both over the tip as the sample. In Figure 4.2 the electrostatic force is plotted as a function of the lift height for the KPFM-plane geometry (gure 3.1b) and the sphere geometry. In both cases an adaptive mesh renement of 25 generations has been deployed. Figure 4.2a shows the electrostatic force of the sphere-plane geometry. The analytic expression (dotted line) has an average deviation from the numerical obtained value (red line) of less then 1 percent. In this case, no signicant improvement was made by integrating over the surface of the sample, instead of the tip. For the KPFM geometry the electrostatic force is plotted in gure 4.2b as a function of lift height. Here the integration over the sample surface (red line) shows to be more smooth than the integration over the tip (blue line).

(23)

(a) sphere-plane geometry (b) KPFM-plane geometry

Figure 4.2: Electrostatic force for the KPFM-plane and sphere-plane geometry

4.2 Capacitance

The capacitance has been determined of various geometries by integrating the electric eld directly at the surface of the tip, as well as a path at a distance from around the tip. For the sphere-plane geometry both integration methods led to a curve very similar to Hudlets approximation, however with an additional constant factor that can be interpreted as the self capacitance. To account for this self capacitance, the capacitance of an individual sphere with a potential with respect to innity was calculated. With the same geometry and mesh, the total capacitance was calculated. For the KPFM geometry, gure 4.3a represents the situation where the sample was simulated by an dielectric constant of = 1000 and a potential of V = 0. The boundary conditions at the border were set to Neumann conditions of zero ux. In gure 4.3b the sample was held at Neumann boundary conditions of zero ux and a dielectric constant of = 1, simulating open space. The boundary conditions at the border were set to Dirichlet boundary conditions, with a xed value treating the tip as a point charge.

(24)

(a) (b)

Figure 4.3: Determination of self capacitance. In gure (a) the sample is represented by a large dielectric constant, in gure (b) the dielectric constant is set to one in the area of the sample, simulating open space.

The same method was applied for both the sphere plane geometry, and the KPFM geometry. For the sphere plane geometry, the mutual capacitance, ie the total capacitance minus the self capacitance, was plotted in gure 4.4a as a function of lift height. The analytic value (dotted line) has an average deviation of 0.7 percent from the numerical obtained value (red line) when integrated over the surface of the tip. Other integration paths, such as integrating over a surface around the tip, results in lower accuracy. For the KPFM geometry the mutual capacitance was plotted in gure 4.4b as a function of the lift height. The numerical determined self capacitance of the KPFM tip is very constant with a standard deviation of 2 · 10−6_{, suggesting accurate measurement. Because}

there was no sample, the integration path could be chosen as seen in gure A.3, that can be found in appendix A. The relative error for the KPFM simulation resulted in approximately the same error as the sphere-plane geometry and the noise in gure 4.4b can be explained by to the relatively short distance between the tip and sample. The determination of the capacitance with the method of integrating along the tip surface seems to have a systematic error.

(25)

(a) sphere plane geometry (b) KPFM plane geometry

(26)

Chapter 5

Conclusion and discussion

In this research numerical calculations have been performed to simulate certain parameters (F, C,∂C ∂z)

of Kelvin Probe Force Microscopy. In the software environment of Matlab, a script has been de-veloped that simulates the electrostatic tip-sample force and other characteristics of the KPFM geometry. The script has been organized and commented, and can be used for further research.

The continuous Poisson equation had been discretized by the use of the nite element method, resulting in discrete values for the electrostatic potential. By use of the Maxwell stress tensor and Gauss's law, we could determine the tip-sample force and capacitance respectively.

For two concentric spheres a high accuracy of up to 0.99999 was reached for ∆C

C and

∆F F ,

depending on the mesh size. Also in the simulations to KPFM geometries reliable data has been generated that does not deviates from the analytical value more than 1 percent. Within this project the electrostatic force was obtained by integrating the Maxwell stress tensor over a surface. Because the tip contains complex geometry including sharp corners, where the potential is not well dened, and innite forces may occur, integration over the tip resulted in an extra error[4]. Accuracy improvement was therefore acquired by integrating over the surface of the sample by making use of the second law of Newton. Figure 4.2b clearly shows that by integrating over the surface of the sample, inaccuracies due to the complex KPFM tip can be corrected for. Other options would be to round of the corners of the tip, at the expense of extra triangles to mesh the curve. Also, instead of using the Maxwell stress tensor method that was used within this project, other methods to determine the force not from the local forces contribution, but somewhere away from the surface would possibly provide accuracy improvement. One such method could be the egg shell method [8], in which also the stress tensor is used, but is integrated over a area surrounding the surface.

Calculations of the capacitance yielded an initially unexpected additional constant that can be interpreted as the self capacitance. This could however, neatly be solved a trick. By calculating the capacitance of the tip, while the sample was eectively grounded due to a large dielectric constant, the self capacitance was determined. This self capacitance represents the ratio of stored charge and potential dierence between the tip and innite sphere surrounding it. To simulate this in the nite sized area of the PDE solver, Dirichlet boundary conditions where applied on the outside borders that obey the potential eld of a point charge. The mutual capacitance can then be obtained by subtracting the self capacitance from the total capacitance. For the sphere plane geometry this resulted in a accurate mutual capacitance that agrees well with the analytical expressions, in which the error is partly due to the error in the approximation of Hudlet. For the more complex geometry

(27)

of the KPFM the self capacitance was relatively constant with a standard deviation in∆C

C of 2·10 −6

for dierent lift heights with a mesh of 3·105_{triangles. This suggests that the calculation for the self}

capacitance is valid and accurate. Calculations for the total capacitance, resulted in an accuracy of

∆C

C of around 1%. Although further improvement is desirable, it is sucient for this project. The

integration path now used follows closely by the sample, which is likely to result in an additional error. Improvement could possibly be made when the integration path is optimized.

Dierent other options are available to improve accuracy. Generally, a greater mesh density results in an improved accuracy. However, since the Matlab solver pdenonlin has an computational complexity of O(n3₎_{[9], large meshes are expensive calculations. Especially adaptive mesh}

rene-ment is a expensive computation, since it requires the PDE problem to be solved in every generation of renement. Figure 4.1 shows adaptive mesh renement quickly becomes unusefull when work-ing with larger meshes. An option would be to divide the solution space into multiple parts, and solve for them separately. This technique is called domain decomposition. This technique has the advantage that the runtime can be reduced, because the number of triangles n in each sub domain is less then the total number of triangles. Also, the mesh density could be increased specially in the area of interest. In our case where dimension dierences are relatively large (z H) this is particularly interesting.

Although dierent improvements are possible and desirable, this project has provided reliable measurements with accuracy of less then 1 percent in ∆C

C . This should be enough to investigate the

eects of artifacts in KPFM measurements, and can therefore be used in further KPFM research to complement experimental data.

(28)

Bibliography

[1] J Collins, L. Kilpatrick. Dual harmonic kelvin probe force microscopy for surface potential measurements of ferroelectrics. 2012.

[2] J. I.; Bhaskaran M. Sriram S. Weber Stefan A. L.; Jarvis Suzi; Rodriguez Brian J. Collins, Liam; Kilpatrick. Dual harmonic kelvin probe force microscopy for surface potential measure-ments of ferroelectrics. Institute of Electrical and Electronics Engineers, 2012.

[3] D. J. Griths. Introduction to Electrodynamics. Pearson, 2013.

[4] V Hannot. S. Rixen, D. J. Rochus. Rounding the corners in an electromechanical fem model. 2007.

[5] M. Guthmann C. Berger J. Hudlet, S. Saint Jean. Evaluation of the capacitive force between an atomic force microscopy tip and a metallic surface. The European Physical Journal B, 2:510, 1998.

[6] I.D. Szente R.N. Jankov, I.R. Goldman. Principles of the kelvin probe force microscopy. 2000. [7] K. Kudo A Kato, H. Asakura. Highly ecient water splitting into H2 and O2 over

lanthanum-doped NaTaO3 photocatalysts with high crystallinity and surface nanostructure. 2003. [8] P. Mach F. Dolezel I. Korous, L. Karban. Ecient eggshell algorithm based on fully adaptive

higher-order nite element method.

[9] Mathworks. Matlab r2014a documentation. 2014.

[10] V. G. Mazauric. From thermostatistics to maxwells equations: A variational approach of electromagnetism. IEEE TRANSACTIONS ON MAGNETICS, 2004.

[11] J. Kummel A. Lee S. Melitz, W. Shen. Kelvin probe force microscopy and its application. Surface Science Reports, 2011.

[12] A. Sadeghi. Multiscale approach for simulations of kelvin probe force microscopy with atomic resolution. 2011.

[13] W. Smythe. Static and dynamic electricity. Pearson, 1950.

[14] G Strang. Computational Science and Engineering. Wellesley-Cambridge Press, 2007.

(29)

Appendix A

Additional gures

(30)

Figure A.2: Integration path used for total capacitance

(31)

Appendix B

Script

This appendix will include the Matlab code that has been used in this research. First a brief overview will be given of the hierarchy of the deployed scripts that are implemented within this project.

script name Child functions Description

Simulation Overall script that denes and runs a simulation

see table B.2 see table B.2

(32)

Function / script

name Child functions Description

params - Generates all settings and start

values for a simulation Methods for analytical values Calculations for geometry

runSimulation - iteratively runs a simulation

GeomBuilder Determines geometry

ScalarBoundaryBuilder Determines boundary conditions

adaptmesh Creates and solves mesh

adaptively (option 1)

initmesh Creates initial mesh

(option 2)

renemesh Renes mesh iteratively

(option 2)

pdenonlin Solves the mesh

(option 2)

numericalCalculations See table B.4

analyticalCalculations All analytical calculations based on the solution of the PDE

self_capacitance See table B.3

Table B.2: All child functions for the script 'Simulation'

Function / script

self_capacitance Determines the self capacitance

geomBuilder Creates geometry

pdenonlin Solves the PDE problem (with the original mesh)

numericalCalculations See table B.4

(33)

Function / script

numericalCalculations Performs all numerical calculations

pdegrad Calculates the gradient of the eld integrateSquarePath Integrates the electric eld over a square path integrateSurfacePath Integrates the electric eld over a surface path mst2 Calculates the electric eld, and integrates it

over the surface Table B.4: All child functions for the script 'numericalCalculations'

1 %example of running an simulation

2

3 %this will create an object that includes all the parameter values. All

4 %values are now at default and can be changed, for example to mode = sphere

5 simulationParams = params();

6 simulationParams.mode = ’sphere’;

7

8 %this function will run the simulation. It requires an object op the type

9 %params. It returns Flist and Clist, two lists with all the information

10 %about the force and capacitance

11 [Flist, Clist] = runSimulation(simulationParams);

1 classdef params < handle

2 % This object provides all the parameters used for simulation

3 %default properties

4 properties

5 % Geometry values and constants

6 scale = 10^-5; %everything gets scaled by this

7 r_scale = 20; %radius sphere

8 hcone_scale = 20000; %hoogte cone

9 hcentilever_scale = 1000; %hoogte centilever

10 rcentilever_scale = 30000; %radius centilever1

11 theta0 = 10/180 ∗ pi; %aperture of sphere

12 xQ_scale = 0.4; %breedte van verschillende potentialen ... op sample

13 yQ_scale = 20; %hoogte van dielectrisch materiaal, ... MOET kleiner dan h zijn

14 yplane_scale = 10000; %height of sample

15 xplane = 0.7; %width of sample

16 % e0 = 8.85∗10^-12; %electric permitivy of vacuum

17 e0 = 1;

18

19 % boundary conditions

20 Vt = ’1’; %Potential of tip

21 Vs = ’0’; %Potential of sample

22 Vb = ’neumann’; % ’dirichlet’ for dirichlet as point charge or ... ’neumann’ for neumann of zero flux at border

23 self_cap_test = false; %Extra option to calculate self-capacitance

24

(34)

26 variation = ’height’; %Either change height every iteration or mesh density

27 h_scale = 5; %starting lift height (h)

28 iterations = 1; %amount of iterations (h changes)

29 step_scale = 5; %amount h increases

30 ngen = 25; %number of generations for adaptive solving

31 adaptive = true; %true of adaptive, false for lineair refinement

32 ntriangle= 3e5; %max triangles in case of lineair refinement

33 gridsize = 5∗10^(-4); %used for integrating over a square path

34 %options for different calculations ...

[’surfacePath’,’squarePath’,’mst’,’tip’,’sample’]

35 calculate = [’surfacePath’,’squarePath’,’tip’,’sample’];

36

37 % model options -> kpfm/sphere/kpfm with centilever/cone

38 mode = ’sphere2’; 39 40 % Plot options 41 %plotparams = [’surfacePath’,’squarePath’,’geometry’,... 42 % ’selfCapacitance’,’solution’]; 43 plotparams = ’’; 44 end 45 methods

46 % methods, including the constructor are defined in this block

47 function obj = params()

48 end

49

50 %methods for constants ON scale

51 function r = r(obj)

52 r = obj.r_scale∗obj.scale;

53 end

54 function hcone = hcone(obj)

55 hcone = obj.hcone_scale∗obj.scale;

56 end

57 function hcentilever = hcentilever(obj)

58 hcentilever = obj.hcentilever_scale∗obj.scale;

59 end

60 function rcentilever = rcentilever(obj)

61 rcentilever = obj.rcentilever_scale∗obj.scale;

62 end

63 function yplane = yplane(obj)

64 yplane = obj.yplane_scale∗obj.scale;

65 end

66 function h = h(obj)

67 h = obj.h_scale∗obj.scale;

68 end

69 function step = step(obj)

70 if obj.adaptive && strcmp(obj.variation,’mesh’)

71 step = obj.step_scale; 72 else 73 step = obj.step_scale∗obj.scale; 74 end 75 end 76 function xQ = xQ(obj) 77 xQ = obj.xQ_scale∗obj.scale; 78 end 79 function yQ = yQ(obj) 80 yQ = obj.yQ_scale∗obj.scale; 81 end 82

(35)

83 %calculations for geometry

84 function value = ybol(obj)

85 value = obj.r_scale∗obj.scale∗(1-sin(obj.theta0)); 86 end 87 function zc = zc(obj) 88 zc = obj.r_scale∗obj.scale∗cos(obj.theta0)/tan(obj.theta0) -... 89 obj.r_scale∗obj.scale∗(1-sin(obj.theta0)); 90 end

91 function xbol = xbol(obj)

92 xbol = obj.r_scale∗obj.scale∗cos(obj.theta0);

93 end

94 function rcone = rcone(obj)

95 rcone = tan(obj.theta0)∗(obj.hcone_scale∗obj.scale+obj.r_scale∗... 96 obj.scale∗(1-sin(obj.theta0)+(obj.r_scale∗obj.scale∗... 97 cos(obj.theta0)/tan(obj.theta0)-obj.r_scale∗obj.scale∗... 98 (1-sin(obj.theta0))))); 99 end 100 function H = H(obj) 101 H = obj.r_scale∗obj.scale∗(1-sin(obj.theta0))+obj.hcone_scale∗... 102 obj.scale; 103 end

104 function hplane = hplane(obj)

105 hplane = (1-(obj.r_scale∗obj.scale∗(1-sin(obj.theta0)) + ... 106 obj.hcone_scale∗obj.scale))/2; 107 end 108 function k = k(obj) 109 k = 1/(log(tan(obj.theta0/2)))^2; 110 end 111 112 end 113 114 end

1 function [Flist, Clist, compTime] = runSimulation(s)

2 %output values 3 Flist = zeros(8,s.iterations); 4 Clist = zeros(22,s.iterations); 5 compTime = zeros(1,s.iterations); 6 7 ngen = s.ngen; 8 h = s.h; 9 ntriangle = s.ntriangle;

10 setenv(’PLOT’,s.plotparams);

11 12 13 14 for i=1:s.iterations 15 time = cputime; 16

17 %set iteration options

18 if strcmp(s.variation,’height’)

19 h = (i-1)∗s.step + s.h;

20 elseif strcmp(s.variation,’mesh’)

21 if s.adaptive

22 ngen = (i-1)∗s.step + s.ngen;

(36)

24 ntriangle = (i-1)∗s.step + s.ntriangle;

25 end

26 end

27

28 %Build geometry and boundary conditions and calculate u

29 [g,b,c,a,f,geom] = geomBuilder2(s.r,s.xbol,s.ybol,h,s.rcone,s.hcone,... 30 s.rcentilever,s.hcentilever,s.H,s.xQ,s.yQ,s.xplane,s.hplane,... 31 s.yplane,s.mode,s.Vt,s.Vs,s.Vb,false); 32 33 %solve for u 34 if s.adaptive

35 [u,p,e,t]=adaptmesh(g,b,c,a,f,’Ngen’,ngen,’Maxt’,5∗10^5 );

36 else

37 [p,e,t]=initmesh(g,’MesherVersion’,’R2013a’);

38 while (size(t,2)<ntriangle)

39 [p,e,t]=refinemesh(g,p,e,t);

40 fprintf(’#of triangles: %d\n’,size(t,2))

41 end

42 [u, res]=pdenonlin(b,p,e,t,c,a,f);

43 if (strfind(getenv(’PLOT’),’solution’))

44 pdeplot(p,e,t,’xydata’,u,’contour’,’on’);

45 end

46 end

47

48 %set default values;

49 Fsurface = 0; Emst = 0; Esurface = 0; EsquarePath = 0; integrate_path = ...

’regular’; Fmst = 0;

50

51 %Calculate force and integral of E on TIP

52 if (strfind(s.calculate,’tip’)) 53 geom.surface = ’tip’; 54 run(’numericalCalculations’); 55 end 56 57 %set values 58 Fsurface_tip = Fsurface; 59 Fmst_tip = Fmst; 60 Qmst_tip = Emst; 61 Qsurface_tip = Esurface; 62 Qsquare_tip = EsquarePath; 63

64 %Calculate force and integral of E on tip on SAMPLE

65 if (strfind(s.calculate,’sample’)) 66 geom.surface = ’sample’; 67 run(’numericalCalculations’); 68 end 69 70 %set values 71 Fsurface_sample = Fsurface; 72 Fmst_sample = Fmst; 73 Qmst_sample = Emst; 74 Qsurface_sample = Esurface; 75 Qsquare_sample = EsquarePath/2; 76

77 %obtain analytic Calculations

78 run(’analyticalCalculations.m’);

79 80

(37)

81 %saving values

82 Flist(1:8,i) = [h; ngen; Fanalytic; Fhudlet2; Fsurface_sample;...

83 Fsurface_tip ; Fmst_sample; Fmst_tip;];

84 Clist(1:10,i) = [h; ngen; Canalytic; Csphere; Qsurface_sample;...

85 Qsurface_tip; Qmst_sample; Qmst_tip; Qsquare_sample; Qsquare_tip];

86 compTime(1,i) = cputime-time;

87

88 %test for self capacitance

89 if s.self_cap_test 90 run(’self_capacitance.m’); 91 end 92 93 end 94 95 Fdiff = -0.5∗gradient(Clist(10,:),s.step); 96

97 filename = [’savedVar\’ s.mode ’ - It ’ mat2str(s.iterations) ’ - Ngen ’ ...

mat2str(s.ngen)];

98 save(filename);

99 end

1 function [gd,b,c,a,f,geom] = geomBuilder2(r,xbol,ybol,h,rcone,hcone,...

2 rcentilever,hcentilever,H,xQ,yQ,xplane,...

3 hplane,yplane,mode,Vt,Vs,Vb,self_cap_test)

4 %This function creates the geometry and boundary condition matrix

5 %The decomposed geometry matrix is speficied manualy

6 %

7 % % h = 0.002; %hoogte sphere

8 % % r = 0.001; %radius sphere

9 % % hcone = 0.15; %hoogte kentilever

10 % % rcone = 0.1; %radius Kentilever

11 % % rcentilever = 0.3; 12 % % hcentilever = 0.01; 13 % 14 %|_rcentilever__ 15 %| _______| hcentilever 16 %|rcone/ 17 %| / hcone 18 %| / 19 %| | 20 %|r / 21 %|_/ 22 %| h 23 %---24 %| |yQ | 25 %---26 % xQ 27 28 if strcmp(mode,’sphere’)

29 % rectangle is code 3, 4 sides,

30 % followed by x-coordinates and then y-coordinates

31 r1 = [3,4,0,1,1,0,0,0,1,1]’;

32 % Circle is code 1, center (0,h), radius r

33 C1 = [1,0,h+r,r]’;

34 % Pad C1 with zeros to enable concatenation with r1

(38)

36 geom = [r1,C1]; 37 ns = (char(’r1’,’C1’))’; 38 sf = ’r1-C1’; 39 40 % Create geometry 41 [gd dt] = decsg(geom,sf,ns); 42 43 %pde coefficients 44 c = ’x’; a = ’0’; f = ’0’; 45

46 %segments and start points for tip

47 geom.tip = [6 7]; 48 geom.sample = 1; 49 geom.border = [2 3]; 50 geom.rotaxis = [4 5]; 51 geom.bsample_in =1; 52 geom.start.sample = 1; 53 geom.start.tip = 6; 54 55 Vborder = 0; 56 Vb = ’neumann’; 57

58 %build boundary condition matrix

59 b = scalarBoundaryBuilder(geom,Vt,Vs,Vb,Vborder,mode,self_cap_test); 60 end 61 if strcmp(mode,’kpfm’) 62 gd = ... 63 [2 2 2 2 2 2 2 1; 64 0 1 1 xbol rcone 0 0 0; 65 1 1 0 rcone 0 0 0 xbol;

66 0 0 1 h+ybol h+ybol+hcone h+ybol+hcone 0 h;

67 0 1 1 h+ybol+hcone h+ybol+hcone 1 h h+ybol;

68 1 1 1 0 0 0 0 0;

69 0 0 0 1 1 1 1 1;

70 0 0 0 0 0 0 0 0;

71 0 0 0 0 0 0 0 h+r;

72 0 0 0 0 0 0 0 r];

73 % onder, rechts, boven, cone, cone boven, boven de cone, onder bol, bol

74

75 %pde coefficients

76 c = ’x’; a = ’0’; f = ’0’;

77

79 geom.tip = [5 4 8]; 80 geom.sample = 1; 81 geom.border = [2 3]; 82 geom.rotaxis = [6 7]; 83 geom.bsample_in =1; 84 geom.start.tip = 8; 85 geom.start.sample = 1; 86 87 Vborder = 0; 88 Vb = ’neumann’; 89

91 b = scalarBoundaryBuilder(geom,Vt,Vs,Vb,Vborder,mode,self_cap_test);

92 93

(39)

94 end

95 if strcmp(mode,’kpfm2’)

96 gd = ...

97 [2 2 2 2 2 2 ...

2 2 ...

2 1; %2 for line, 1 for sphere

98 0 1 1 xbol rcone rcentilever ...

rcentilever 0 0 0; %x1

99 1 1 0 rcone rcentilever rcentilever ...

0 0 0 xbol; %x2

100 0 0 1 h+ybol h+ybol+hcone h+ybol+hcone ...

h+ybol+hcone+hcentilever h+ybol+hcone+hcentilever 0 h; %y1

101 0 1 1 h+ybol+hcone h+ybol+hcone h+ybol+hcone+hcentilever ...

h+ybol+hcone+hcentilever 1 h h+ybol;%y2

102 1 1 1 0 0 0 ... 0 0 0 0; %left side ... region 103 0 0 0 1 1 1 ... 1 1 1 1; %right ... side region 104 0 0 0 0 0 0 ... 0 0 0 0; %x center ... sphere 105 0 0 0 0 0 0 ... 0 0 0 h+r; %y center ... sphere 106 0 0 0 0 0 0 ... 0 0 0 r]; %radius sphere

107 %onder rechts boven cone centilever onder centilever rechts ... centilever_boven boven centilever onder bol bol

108

110 c = ’x’; a = ’0’; f = ’0’;

111

113 geom.tip = [4 5 6 7 10]; 114 geom.sample = 1; 115 geom.border = [2 3]; 116 geom.rotaxis = [8 9]; 117 geom.bsample_in = 1; 118 geom.start.tip = 10; 119 geom.start.sample = 1; 120 121 Vborder = 0; 122 Vb = ’neumann’;

124 b = scalarBoundaryBuilder(geom,Vt,Vs,Vb,Vborder,mode,self_cap_test); 125 126 end 127 if strcmp(mode,’bol’) 128 gd = ... 129 [2 2 4 4 4 4; 130 0 0 6.1232339957367663E-18 0.1 5.5109105961630896E-17 0.9; 131 0 0 0.1 6.1232339957367663E-18 0.9 5.5109105961630896E-17; 132 0.9 -0.10000000000000024 -0.1 0 -0.9 0; 133 0.10000000000000024 -0.9 0 0.1 0 0.9; 134 1 1 0 0 1 1; 135 0 0 1 1 0 0;

(40)

136 0 0 0 0 0 0; 137 0 0 0 0 0 0; 138 0 0 0.1 0.1 0.9 0.9; 139 0 0 0.1 0.1 0.9 0.9; 140 0 0 0 0 0 0]; 141 142 b = ... 143 [1 1 1 1 1 1; 144 0 0 1 1 1 1; 145 1 1 1 1 1 1; 146 1 1 1 1 1 1; 147 48 48 1 1 1 1; 148 48 48 1 1 1 1; 149 48 48 48 48 48 48; 150 48 48 48 48 48 48; 151 49 49 49 49 49 49; 152 48 48 49 49 48 48]; 153 154 %pde coefficients 155 c = ’x’; a = ’0’; f = ’0’; 156

158 geom.tip = [3 4]; 159 geom.start = 2; 160 end 161 162 if strcmp(mode,’cone’) 163 gd = ... 164 [2 2 2 2 2 2 2; 165 0 1 1 0 rcone 0 0; 166 1 1 0 rcone 0 0 0; 167 0 0 1 h h+hcone h+hcone 0; 168 0 1 1 h+hcone h+hcone 1 h; 169 1 1 1 0 0 0 0; 170 0 0 0 1 1 1 1; 171 0 0 0 0 0 0 0; 172 0 0 0 0 0 0 0; 173 0 0 0 0 0 0 0];

174 % onder, rechts, boven, cone, cone boven, boven de cone, onder cone,

175

177 c = ’x’; a = ’0’; f = ’0’;

178

180 geom.tip = [4 5]; 181 geom.start = 5; 182 geom.border = [2 3]; 183 geom.rotaxis = [6 7]; 184 geom.sample.tip = 1; 185 % geom.sample.sample = 1; 186 187 b = scalarBoundaryBuilder(geom,Vt,Vs,Vb,Vborder,mode,self_cap_test); 188 end 189 190 if (strcmp(mode,’diffV’)) 191 xv1 = 0.1; 192 gd = ... 193 [2 2 2 2 2 2 2 2 2 2 1;

(41)

194 0 xv1 1 1 xbol rcone rcentilever rcentilever 0 0 0;

195 xv1 1 1 0 rcone rcentilever rcentilever 0 0 0 xbol;

196 0 0 0 1 h+ybol h+ybol+hcone h+ybol+hcone ...

h+ybol+hcone+hcentilever h+ybol+hcone+hcentilever 0 h;

197 0 0 1 1 h+ybol+hcone h+ybol+hcone ...

h+ybol+hcone+hcentilever h+ybol+hcone+hcentilever 1 h h+ybol;

198 1 1 1 1 0 0 0 0 0 0 0;

199 0 0 0 0 1 1 1 1 1 1 1;

200 0 0 0 0 0 0 0 0 0 0 0;

201 0 0 0 0 0 0 0 0 0 0 h+r;

202 0 0 0 0 0 0 0 0 0 0 r];

203 %onder %onder2 %rechts %boven %cone %centilever_onder ...

%centilever_rechts %centilever_boven %boven_centilever %onder_bol %bol

204 end 205 206 if (strcmp(mode,’diffQ’)) 207 %geometry 208 gd = ... 209 [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1; %2 for ... line, 1 for sphere

210 1 xbol rcentilever rcentilever 0 xQ xQ 0 xQ 1 1 rcone ...

0 0 0 0; %x1

211 0 rcone rcentilever 0 xQ 1 xQ xQ 1 1 1 rcentilever 0 ...

0 0 xbol; %x2

212 1 h+ybol h+ybol+hcone h+ybol+hcone+hcentilever yQ yQ 0 0 ...

0 0 yQ h+ybol+hcone 0 h+ybol+hcone+hcentilever yQ h; %y1

213 1 h+ybol+hcone h+ybol+hcone+hcentilever ...

h+ybol+hcone+hcentilever yQ yQ yQ 0 0 yQ 1 ...

h+ybol+hcone yQ 1 h h+ybol; %y2

214 1 0 0 0 1 1 2 2 3 3 3 0 0 0 0 0; %left side ... region 215 0 1 1 1 2 3 3 0 0 0 0 1 2 1 1 1; %right ... side region 216 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; %x center ... sphere 217 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 h+r; %y ... center sphere 218 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 r]; %radius sphere

219 %boven %cone %centilever_rechts %centilever_boven %Q_x1 %Q_x2 ... %midden_q %onder %onder2 %Q_rechts %rechts %centilever_onder ... %Q_links %boven_centilever %onder_bol %bol

220 221 % boundary 222 b = ... 223 [1 1 1 1 0 0 0 1 1 1 1 1 1 1 1 1; 224 0 1 1 1 1 1 1 1 1 0 0 1 0 0 0 1; 225 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1; 226 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1; 227 48 1 1 1 1 1 1 1 1 48 48 1 48 48 48 1; 228 48 1 1 1 1 1 1 1 1 48 48 1 48 48 48 1; 229 48 48 48 48 48 48 48 48 48 48 48 48 48 48 48 48; 230 48 48 48 48 48 48 48 48 48 48 48 48 48 48 48 48; 231 49 49 49 49 49 49 49 49 49 49 49 49 49 49 49 49; 232 48 49 49 49 48 48 48 48 48 48 48 49 48 48 48 49]; 233 234 %PDE coefficients 235 a = ’0.0!0.0!0.0’; 236 c = ’x!x!x’;

(42)

237 f = ’0!1!0’; 238 c = ’x’; a = ’0’; f = ’0’; 239 240 end 241 if strcmp(’freeSample’,mode) 242 243 gd = ... 244 [2 2 1 2 2 2 2 2 2 2 2 2 2 2 2; %2 for ... line, 1 for sphere

245 rcentilever 0 0 0 xplane xplane 0 1 1 xbol rcone ...

246 0 0 xbol xplane xplane 0 0 1 0 rcone rcentilever ...

247 h+ybol+hcone+hcentilever+yplane+hplane ...

h+ybol+hcone+hcentilever+yplane+hplane h+yplane+hplane ...

yplane+hplane yplane+hplane hplane hplane 0 1 ...

h+ybol+yplane+hplane h+ybol+hcone+yplane+hplane ...

h+ybol+hcone+yplane+hplane yplane+hplane 0 0; %y1

248 h+ybol+hcone+hcentilever+yplane+hplane 1 h+ybol+yplane+hplane ...

yplane+hplane hplane hplane hplane+yplane 1 1 ...

h+ybol+hcone+yplane+hplane h+ybol+hcone+yplane+hplane ...

h+ybol+hcone+hcentilever+yplane+hplane h+yplane+hplane hplane ...

0; %y2 249 0 0 0 1 1 1 0 1 1 0 0 0 0 0 1; %left side ... region 250 1 1 1 2 2 2 2 0 0 1 1 1 1 1 0; %right ... side region 251 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; %x center ... sphere 252 0 0 h+r+yplane+hplane 0 0 0 0 0 0 0 0 0 0 ... 0 0; %y center sphere 253 0 0 r 0 0 0 0 0 0 0 0 0 0 0 0]; %radius sphere

254 %centilever_boven %boven_centilever %bol %sample_boven ...

%sample_rechts %sample_onder %sample_links %rechts %boven %cone ... %centilever_onder %centilever_rechts %onder_bol %links_onder ... %onder 255 256 c = ’x!1000∗x’; 257 a = ’0!0’; 258 f = ’0!0’; 259 260 geom.tip = [1 3 10 11 12]; 261 geom.border = [8 9 15]; 262 geom.rotaxis = [2 13 14]; 263 geom.sample = 7; 264 geom.bsample_in = [4 5 6]; 265 geom.start.tip = 2;

266 geom.start.sample = 6;%find(p(1,:)==0 &p(2,:)==(hplane+yplane));

267

268 r0 = sqrt((rcone∗H/2+hcentilever∗rcentilever)/pi);

269 Vborder = [mat2str(r0∗(Vt-48)) ’./sqrt(x.^2+(y-’ ...

mat2str(2∗H/3+r+h+yplane+hplane+hcentilever/2) ’).^2)’]; 270 271 b = scalarBoundaryBuilder(geom,Vt,Vs,Vb,Vborder,mode,self_cap_test); 272 273 end 274 if strcmp(mode,’sphere2’) 275

Simulations for Kelvin Probe Force Microscopy