Reconstruction of shear force in Atomic Force Microscopy from measured displacement of the cone-shaped cantilever tip

In this paper, a dynamic model of reconstruction of the shear force $g(t)$ in the Atomic Force Microscopy (AFM) cantilever tip-sample interaction is proposed. The interaction of the cone-shaped cantilever tip with the surface of the specimen (sample) is modeled by the damped Euler-Bernoulli beam equation $\rho_A(x)u_{tt}$ $+\mu(x)u_{t}+(r(x)u_{xx}+\kappa(x)u_{xxt})_{xx}=0$, $(x,t)\in (0,\ell)\times (0,T)$, subject to the following initial, $u(x,0)=0$, $u_t(x,0)=0$ and boundary, $u(0,t)=0$, $u_{x}(0,t)=0$, $\left (r(x)u_{xx}(x,t)+\kappa(x)u_{xxt} \right )_{x=\ell}=M(t)$, $\left (-(r(x)u_{xx}+\kappa(x)u_{xxt})_x\right )_{x=\ell}=g(t)$ conditions, where $M(t):=2h\cos \theta\,g(t)/\pi$ is the momentum generated by the transverse shear force $g(t)$. For the reconstruction of $g(t)$ the measured displacement $\nu(t):=u(\ell,t)$ is used as an additional data. The least square functional $J(F)=\frac{1}{2}\Vert u(\ell,\cdot)-\nu \Vert_{L^2(0,T)}^2$ is introduced and an explicit gradient formula for the Fr\'echet derivative through the solution of the adjoint problem is derived. This allows to construct a gradient based numerical algorithm for the reconstructions of the shear force from noise free as well as from random noisy measured output $\nu (t)$. Computational experiments show that the proposed algorithm is very fast and robust. This allows to develop a numerical"gadget"for computational experiments of generic AFMs.


Introduction
Micro-cantilever plays a key role in nanomachining process using an AFM which was originally developed to provide surface topography information [1].Nowadays, AFM can provide high resolution images in different settings including ambient, aqueous and vacuum environments.In standard AFMs, the micro-cantilever is mounted horizontally and the devices are operated in a contact or intermittent-contact mode (Fig. 1).The cantilever tip-sample interaction creates a transverse shear force and a bending moment on the tip of the cantilever [2].Estimation of the unknown shear force signal allows better interpretation and understanding of scan results.Since this force can only be measured indirectly, via a laser based sensor system, various models and inversion algorithms was developed for reconstruction of the transverse shear force in atomic and dynamic force microscopy, through the measured cantilever tip deflection (see [3,4,5] and references therein).
For the AFM cone-shaped cantilever tip-sample interaction, a simple mathematical model for the shear force reconstruction problem has first been proposed in [4], within the Euler-Bernoulli beam theory.Namely, the model considers the cutting system as the inverse problem of reconstructing the cutting force F y (t) in EI y tt = 0, x ∈ (0, L), y(x, 0) = y t (x, 0) = 0, x ∈ (0, L), y(0, t) = y x (0, t) = 0, from the measured displacement Y (L, t) := y(L, t), t ∈ [0, T ]. ( using available displacement measurement.Here, F x = (2h cos θ F y ) /π for a cone-shaped cantilever with the half-conic angle θ, and h > 0 is the cantilever tip length.This is an inverse problem with two Neumann inputs.Note that a similar inverse problem with one Neumann input was proposed in [6].A detailed analysis of inverse problems of identifying the unknown transverse shear force in the Euler-Bernoulli beam with Kelvin-Voigt damping was given in [7].
It is important to emphasize that in the AFM cone-shaped cantilever tipsample interaction model ( 1) is based on the simplified and constant coefficient Euler-Bernoulli beam equation, without the viscous external (µ(x)u t ) and the internal or Kelvin-Voigt ((κ(x)u xxt ) xx ) damping terms.Thus, in these models, not all physical properties of the cantilever are taken into account.However, the influence of these above mentioned properties on the dynamic behavior of the AFM cantilever is enormous, and needs to be studied carefully [9,10].
In this paper we propose a mathematical model of tip-sample processing in AFM with two Neumann inputs.This model is a generalization of existing mathematical models in the sense that; (a) the Euler-Bernoulli equation contains all the physical variable coefficients, including the both damping terms; (b) the time interval during which it is necessary to produce an experimental data, i.e. measured output, can be small enough; (c) the measured output contains random noise; (d) the inputs in the model may not be smooth enough.
Within the proposed model, we formulate the inverse problem of reconstructing the unknown shear force from measured displacement of the coneshaped cantilever tip.We provide a detailed mathematical and numerical analysis of the problem.Based on this analysis, we derive an explicit gradient formula for the least square functional.This allows us to construct an effective and fast reconstruction algorithm, as the presented results of computational experiments show.

Vibration model of tip-sample processing: the reconstruction problem
The sample processing with AFM cone-shaped cantilever, shown schematically in Fig. 1, is modeled as a damped Euler-Bernoulli beam.This cantilever, with length > 0 and cross-sectional area A s (x), is clamped at the left end x = 0.The tip-sample contact is modeled by a vertical reaction force, which is the transverse shear force with the negative sign, that is −g(t), and the moment M (t) := −2h cos θ g(t)/π, generated by this force, where h, θ > 0 are the tip length and half-conic angle, respectively.Then the sample processing vibration model is governed by the following initial boundary value problem for the damped Euler-Bernoulli equation: where Ω T := (0, ) × (0, T ), and the final time instance T > 0 may be small enough.Here and below, ρ A (x) := ρ(x)A s (x), while ρ(x) > 0 and A s (x) > 0 are the mass density and the cross-sectional area of the nonhomogeneous cantilever, r(x) := E(x)I(x) > 0 is the flexural rigidity (or bending stiffness) of the cantilever while E(x) > 0 is the elasticity modulus and I(x) > 0 is the moment of inertia.The coefficient κ(x) := c d (x)I(x) represents energy dissipated by friction internal to the beam, while c d > 0 is the strain-rate damping coefficient [9].The external and internal damping mechanisms are given by the terms µ(x)u t and (κ(x)u xxt ) xx , respectively.The coefficients µ(x) ≥ 0 and κ(x) > 0 are called the viscous (internal) damping and the strain-rate or Kelvin-Voigt damping coefficients, respectively.The transverse shear force g(t) in ( 3) is assumed to be unknown and needs to be determined from knowledge of the measured displacement ν(t) of the cone-shaped tip: Thus the inverse boundary value problem here is to reconstruct the unknown transverse shear force g(t) in (3) from knowledge of the measured displacement ν(t) defined in (4).
As noted above, the model governed by (3) and ( 4) describes an inverse problem with two inputs M (t) and g(t).As we shall see, this results in a number of differences and additional problems, unlike the single-input inverse problems considered in [6,7].Note also that due to the conical geometry of the rod, the relationship M (t) := (2h cos θ g(t)) /π is defined between the inputs.
We assume that the following basic conditions are satisfied: Introduce the set of admissible shear forces where C g > 0 is a constant independent on g(t).Denote by u(x, t; g) the solution of the forward problem (3) for a given g ∈ G, while u( , t; g) in defined as an output.Introduce the Neumann-to-Dirichlet operator: (Ψg)(t) := u( , t; g), t ∈ [0, T ], Ψ : G ⊂ H 1 (0, T ) → L 2 (0, T ), (7) defined on the set of admissible shear forces.In view of this operator, we can reformulate the inverse problem as the linear operator equation: Since the measured output ν(t) obtained as a result of measurement, it contains random noise.Hence the exact equality between the output u( , t; g) and the measured outputs ν(t) can never be achieved.As a consequence, there can never be an exact solution to the inverse problem (3)-( 4).We introduce the Tikhonov functional and look for the quasi-solution of the inverse problem (3)-( 4): Find g ∈ G such that

Necessary estimates for the weak solution of problem (3)
In the case when M (t) = 0, the existence and uniqueness of the weak solution u ∈ L 2 (0, T ; V 2 (0, )), with u t ∈ L 2 (0, T ; L 2 (0, )) and u tt ∈ L 2 (0, T ; H −2 (0, )) of the initial boundary value problem ( 3) is proved in [7], where For the direct problem (3) the same results can be proved in the same way.We derive here some a priori estimates for the weak solution which are necessary in the analysis of the inverse problem ( 3)-( 4).
Theorem 1. Assume that the inputs in (3) satisfy the basic conditions (5).Then the following estimates holds: where and r 0 , ρ 0 , κ 0 > 0 are the constants introduced in (5) Proof.Multiply both sides of equation ( 3) by 2u t (x, t), integrate it over Ω t := (0, ) × (0, t), t ∈ (0, T ], and employ the identities Applying the integration by parts formula multiple times, using the initial and boundary conditions in (3) we obtain the following energy identity: for all t ∈ (0, T ].Now, applying the integration by parts to the right hand side integrals and the use the ε-inequality, we have: Use also the auxiliary inequalities ( [11], Ch. 11): for all t ∈ [0, T ].Taking these inequalities with the inequality into account in the energy identity above, we get the following inequality: where , C θ > 0 are the constants introduced in (12).Choosing here the arbitrary parameter ε > 0 from the condition r 0 − ε > 0 as ε = r 0 /(2 ) we finally obtain the main integral inequality: where C 0 > 0 is the constant introduced in ( 12).The first consequence of ( 14) is the inequality With the Grönwall-Bellmann inequality this implies: Both of the first two estimates in (11) are easily derived from this inequality.The second consequence of ( 14) is the inequality With (15) this leads to the third estimate in (11).The fourth estimate in (11) is proved in the same way.
Remark 1.The results of Theorem 1 are valid, with slightly different from the constants introduced in ( 12), also for the case where the consistency condition g(0) = 0 in ( 5) is not met.
Corollary 1. Assume that conditions of Theorem 1 hold.Then for the H 1norm of the output u( , t; g) the following trace estimate holds: Proof follows from the trace inequalities which are the consequence of the first inequality in (13) and estimates in (11).

Analysis of the inverse problem
The compactness property is one of the main properties of the inputoutput operators corresponding to problems, since the ill-posedness of an inverse problem is the result of this property.For the simplified version, with one Neumann input (M (t) = 0) and with κ(t) = 0, the compactness of the Neumann-to-Dirichlet operator ( 7) is proven in [6] for the regular weak solution.For the model we are considering, the regularity condition is not necessary, as we shall see below.That is, this property is also preserved in the case of the weak solution, which shows the role of the Kelvin-Voigt damping coefficient κ(x) > 0.
Lemma 2. Assume that the basic conditions (5) hold.Then the Neumannto-Dirichlet operator is Lipschitz continuous, that is with here L 0 = 3 /3 C 1 > 0 is the Lipschitz constant and C 1 > 0 is the constant introduced in (12).
Proof.Let u k (x, t) := u(x, t; g k ), k = 1, 2, be two weak solutions of the direct problem (3) corresponding to the inputs g 1 , g 2 ∈ G. Then the function δu(x, t) = u 1 (x, t) − u 2 (x, t) solves the problem subject to the inputs δM (t) = C 2 θ δg(t) and δg(t) = g 1 (t) − g 2 (t).By the definition (7) of the input-output operator we have: In view of the first inequality in (13) and the second estimate in (11) applied to the weak solution δu(x, t) of problem (18) we deduce that δu( , This leads to (17).
The Lipschitz continuity of the Neumann-to-Dirichlet operator leads to the Lipschitz continuity of the Tikhonov functional introduced in (9), and this, in turn, leads to the existence of the quasi-solution of the inverse problem (3)-( 4), by Theorem 6.5.2 [11].
Theorem 2. Assume that the inputs in (3) satisfy the basic conditions (5).Suppose that the measured output ν(t) belongs to L 2 (0, T ).Then there exists a quasi-solution of the inverse problem (3)-( 4) in the set of admissible shear forces G.

Fréchet differentiability of the Tikhonov functional and gradient formula
For g, g + δg ∈ G we find the increment δJ(g) := J(g + δg) − J(g) of the Tikhonov functional introduced in (9) is where δu(x, t) is the solution of the sensitivity problem (18).Multiplying both sides of equation ( 18) by arbitrary function φ(x, t), integrating it over (0, T ) and applying the integration by parts formula multiple times, we obtain: We require now φ(x, t) solves the well-posed backward problem The control function p(t) here is the arbitrary Neumann input and is specified below.
In view of the initial, final and boundary conditions in (3) and (18) we deduce from (22) the following integral relationship: where C θ > 0 is the constant introduced in ( 12).Taking into account the increment formula (20) we choose the control function p(t) as follows: The backward problem with this input, i.e. the problem is called the adjoint problem corresponding to the inverse problem (3)-(4).Substituting (24) into (23) we obtain the input-output relationship: which contains the output u( , t; g) and the measured output ν(t).Comparing ( 20) and ( 26) we deduce that Theorem 3. Assume that the inputs in (3) satisfy the basic conditions (5).Suppose, in addition, the measured output ν(t) belongs to H 1 (0, T ).Then the Tikhonov functional introduced in ( 9) is Fréchet differentiable.Furthermore, for the Fréchet gradient of this functional the following gradient formula holds: Proof.Applying the first inequality in (13) with the second estimate in (11) to the weak solution δu(x, t) of problem (18) we conclude that the second right hand side integral in (27) of the order O( g L 2 (0,T ) ).This means that the Tikhonov functional is Fréchet differentiable.The gradient formula (28) expressed in terms of the weak solution φ(x, t) of the adjoint problem (25) forms the basis of the algorithm for numerical solving the inverse problem (3)-(4).

Numerical Algorithms and Computational Experiments
In this section, a detailed description of an efficient numerical method is presented to solve the inverse problem (3)-( 4).This process has several steps and each of them should be considered carefully due to the sensitivity of the identification process.First, measured data ν(x) := u( , t) is generated by solving the direct problem.It is critical to keep the error as low as possible in this step.This requires a successful algorithm for the solution of the direct problem (3).Due to the effectiveness of the method of lines approach used in our several published previous studies ( [6], [12], [13], [14]) on an optimized mesh, an improved version of this method is employed here.

The Method of Lines (MOL) Approach for the Numerical Solution of
Direct Problem Basically, the MOL is based on the principle of independent discretization of space and time variables.More specifically, a semi-analytical structure is obtained by expressing the variational formulation in finite dimensional space denoted by V h .The method here is a finite element approximation with cubic Hermite basis functions which ensures continuity of both deflection and slope throughout the beam.These shape functions is defined on uniformly discretizing spatial domain 0 = Formally, the solution U h (t) := u h (•, t) ≈ u(•, t) satisfies the following semidiscrete version of the variational formulation of (3).
Here the symmetric bilinear functional a ψ : The next discretization step is performed for temporal derivatives.At this level the second order system of ODE in (29) can be approximately solved by using any temporal finite difference method.It is crucial that the approach to be used here have to be practical, fast and stable.These requirements can be met through the following second order backward finite difference approximations of U h and U h with uniform temporal discretization 0 = t 0 < t The full-discrete algebraic systems of equations are obtained by substituting these difference expressions with U h (t) and U h (t) in (29).Solutions of the resulted equations are provided desired approximations U j h ≈ u(x, t j ) for j = 0 : N .Note that for j = 1, 2, the necessary a priori approximations can be obtained by combining the ghost point technique within the central difference scheme.
Finally, several numerical tests are compared to determine the effective values of the pair (M, N ) and optimized with the ratio h/τ 142.

Reconstruction with Conjugate Gradient Algorithm (CGA)
The explicit gradient formula in (28) is very important in determining the minimizer of the least square functional (7) for any unconstrained optimization techniques.Here we use CGA, one of the most suitable and stable one.It is known that this method is based on the conjugate directions and these directions are determined by the solution of the adjoint problem (25).This requires the MOL technique at each iteration step.Although CGA is a self-stabilized method, the quality of the reconstruction process also depends on the success of solving both adjoint and direct problem.The details of the CGA is as follows.
• From g (i) (t), calculate the decent direction • Define the next iteration g (i+1) (t) = g (i) (t)+α (i) * p (i) (t).Here α * i solution of the minimization problem and has the following explicit form, • If the following stopping condition based on Mozorov's discrepancy principle holds, for known parameter ε > 0, stop the iteration; otherwise, repeat the process by taking g (i) (t) := g (i+1) (t).
For the first iteration, an arbitrary choice of g (0) (t) can be made, but if there is no prior knowledge it is better to choose g (0) (t) = 0 and p (i) (t) := −∇J(g (0) )(t).As a note, the first iteration has no significant effect to the success of the algorithm.
Here a standard method is used for the derivation of synthetic noise with a given noise level γ > 0. In deed, the formula ν γ (t j ) := ν(t j )+γ ν L 2 (0,T ) R j for j = 1, • • • , N generates measured noisy data.The vector R has M random numbers array normally distributed with mean 0 and standard deviation σ = 1.
In CGA steps, both Fréchet derivative ∇J(g) and L 2 (0, T ) norms are computed by Simpson's numerical integration while in MOL algorithm, a three-point Gauss quadrature rule is employed for all computation on each element.

Computational Experiments
In the reconstruction process, we work on two different test problems.One of them is based on engineering applications (realistic parameters), while the other one is preferred to test the applicability of the method.
It is a general approach to use error analysis when comparing the quality of the methods.In the literature, two quantities are frequently used.These are Convergence and Accuracy Errors as follows.
As can be seen from their definitions, the Accuracy Error determines the success of the reconstruction.On the other hand, especially in the case of noisy data, the stop criterion is very crucial to prevent divergence of the approximation and it is completely related to Convergence Error.Therefore, these quantities should be evaluated together to analyze the process.
For the first test problem, the parameters are selected in accordance with real engineering applications and are based on those proposed in [15,16,17].We take a beam of length of 200 nm and observe it for a time interval of 10 −3 s.After a simple change of variables, to re-scale the problem so that the length of the beam and the time observation length interval become = 1 and T = 1 respectively, the numerical values adopted for this study become as follows: ρ A (x) = 1.864 × 10 −7 kg/nm, µ(x) = 8.16 × 10 −6 kg s −1 /nm, r(x) = 2.265 × 10 −3 kg nm 3 /s −2 , κ(x) = 3.5875 × 10 −5 kg nm 3 /s, and domain parameters are = 1 and T = 1, both non-dimensional.As for the tip length, it usually ranges from 5 nm to 50 nm [2].After the re-scaling for doing our numerical simulations, we take as a reasonable value h = 0.2 (non-dimensional).
We tested the performance of the algorithm for the unknown shear force g(t) = t sin(7πt/2) with θ π = (cos(π/36)) /(5π).The graph on the left in Fig. 2 shows noisy free as well as random noisy output data with the noise levels γ = 3% and 6%.Then unknown target g(t) is identified by using each of these data.Results can be seen on the right in Fig. 2  Fig. 3 reveals the general characteristics of an iteration.Especially the rapid deterioration in the Accuracy Error indicates that the sensitivity of the stopping which is directly determined by the Convergence Error.In case this balance is not determined appropriately, the success of the construction process can be adversely affected.The second computational experiment aims to test the accuracy of CGA regardless of the realizability of the parameters.For this goal, reconstruction of the following discontinues target source g(t) is studied under high noise levels.
) and θ π = 2 cos(π/4) π Here H(x) is the Heaviside step function.Moreover, all problem parameters are imposed as non-constant case as follows with unit domain parameters = 1 and T = 1.
Synthetic noise free and noisy data are plotted in Fig. 4 (left) with noise levels γ = 5% and 10%.Then CGA is applied for identification of the temporal function g(t) and results are illustrated in Fig. 4 (right).Here, due to the effect of high noise levels and discontinuity on g(t), non-physical distortions are naturally observed in the reconstruction.
Convergence and Accuracy Errors are plotted in Figure 5 on the left and on the right, respectively.Similar behavior of these error quantities examined in the first problem is also observed in this second experiment.The results of the two experiments presented here show that CGA is effective and successful for the solution of the inverse problem under considera-tion, provided that certain sensitivities are taken into account.Nevertheless, the algorithm may need to be improved for further applications.Especially in the realistic cases, it is required to choose a small final time such as T = 10 −3 for stable calculations using the Finite Element Method.Since the method suggested here is just a preliminary numerical study of the inverse problem related to Atomic Force Microscopy, we only aimed to present the general principles.

Conclusions
In this study, a novel mathematical model of tip-sample processing with AFM cone-shaped cantilever is proposed.Compared to the models known in the literature, this model is a fairly advanced model, and takes into account both viscous and internal damping parameters.A detailed mathematical analysis of the model has been carried out.An explicit gradient formula for the Fréchet derivative of Tikhonov functional is derived through the weak solution of the appropriate adjoint problem.This allows us to construct the fast Conjugate Gradient Algorithm for the numerical reconstruction of the shear force.Numerical experiments carried out with real physical and geometric parameters show the high accuracy of the algorithm.

Figure 1 :
Figure 1: Schematic diagram of AFM cone-shaped cantilever tip-sample interaction