A comprehensive model for transient behavior of tapping mode atomic force microscope

Many investigations have focused on steady-state nonlinear dynamics of cantilevers in tapping mode atomic force microscopy (TM-AFM). However, a transient dynamic model—which is essential for a model-based control design—is still missing. In this paper, we derive a mathematical model which covers both the transient and steady-state behavior. The steady-state response of the proposed model has been validated with existing theories. Its transient response, however, which is not covered with existing theories, has been successfully verified with experiments. Besides enabling model-based control design for TM-AFM, this model can explain the high-end aspects of AFM such as speed limitation, image quality, and eventual chaotic behavior.


Introduction
Atomic force microscopy (AFM) is a powerful tool for imaging the surface of samples with a sub-nanometer resolution which has many different applications in experimental sciences such as physics and biology [1][2][3][4]. Nowadays, AFM is also attracting more and more attention in the semiconductor industry as an inspection and metrology tool. The latter applications demand much stricter requirements regarding throughput, accuracy, reliability, and nondestructiveness [5][6][7][8].
To improve the imaging throughput and reliability of AFM, it is essential to understand its physics and optimize its mechanical parts and controllers.
A popular and promising AFM-based measurement technique, both for research and industrial applications, is the amplitude modulation (AM-AFM) or the tapping mode AFM (TM-AFM). TM-AFM is mainly popular because it applies very small forces on the sample surface and, consequently, has lower probability of damaging the sample. The nondestructiveness of TM-AFM can be so important that, e.g., researchers in the field of biology refer to it as a blessing [9].
The TM-AFM works as follows: A cantilever with a sharp tip attached to its free end is excited at a frequency around its fundamental resonance frequency, and brought close to the surface of the sample. At a Fig. 1 Schematic representation of relationship between amplitude and distance, depicting the reasoning for parachuting effect, coexistence of attractive-repulsive regimes, snap-in phenomenon and transient behavior certain distance from the sample surface, the tip starts to interact with the sample via different forces such as van der Waals (vdW) interaction, Pauli repulsion, squeezed film damping, hydro-capillarity, electrostatic forces, etc [10]. As a result of the tip-sample interaction (TSI) forces, the amplitude of the motion of the cantilever reduces [11,12]. A feedback control loop adjusts the distance between the cantilever and the sample surface using a piezoelectric (z-stage) actuator, such that the amplitude of the cantilever is kept constant at a userdefined fraction of its free air amplitude. While raster scanning the sample surface, and keeping the amplitude constant, the output signal of the controller is recorded and interpreted as the height profile of the sample.
In case the amplitude of the cantilever has a oneto-one relationship with its distance to the surface, keeping the amplitude constant would mean keeping the distance constant. Therefore, any fluctuation in the height profile of the sample is compensated with the z-stage control unit, and the height profile of the sample can be measured via the control signal. However, as schematically shown in Fig. 1, the relationship between the amplitude and the distance is not always one-toone, and consequently, the height measurement is not always trivial. For example, at the far right-hand side of Fig. 1, where the sample is far from the cantilever, obviously there is no interaction and the amplitude is independent of the distance. If this situation happens during imaging, the error signal saturates and the so-called parachuting-type artifacts appear on the image [4,13]. Similarly, at the far left-hand side of the curve in Fig. 1, where the attractive forces exceed a certain value, the tip snaps onto the surface and sticks to it, the so-called snap-in phenomena. Thus, independent of small varia-tions of the distance, the amplitude remains zero [14]. Another situation where amplitude and height do not have a one-to-one relationship is the coexistence of the attractive and repulsive regimes [15]. At two different distances from the surface, the TSI force can reduce the amplitude of the cantilever to the same value. In one regime, the average force is attractive, and in the other, it is repulsive. In Fig. 1, if the distance between the cantilever and the sample is h 1 , the vibration amplitude could be either A 1 or A 2 . Similarly, if the measured amplitude is equal to A 2 , the distance could be either h 1 or h 2 . In this case, a certain amplitude may correspond to two different height values and vice versa. Hence the controller fluctuates between the two height values, which causes artifacts on the image. Researchers have studied these nonlinear effects, their consequences, and the related problems, extensively [15][16][17]. However, all these studies consider a steady-state condition and the transient response of the cantilever is neglected.
The nonlinearities are not the only cause for distortion of the amplitude-distance relationship, also transient behavior of the cantilever can distort this one-toone relationship. The amplitude-distance relationship in Fig. 1 is only valid for steady-state situations, i.e., where the amplitude signal is settled to a constant value. For example, suppose that the system is initially settled to Point A in Fig. 1, then, suddenly, the surface is retracted up to a distance corresponding to Point B. Obviously, the amplitude of the cantilever cannot adapt suddenly. Therefore, it might follow a different trajectory than the linear steady-state one. How fast the amplitude can adjust itself to the variations in the height totally depends on its dynamic trajectory. This problem is crucial for high-speed AFM, where the changes of the distance happen in time intervals that are comparable or shorter than the cantilevers response time. Note that, the transient response of the cantilever is not taken into account in control design for AFM. Hence, such a transient behavior of the cantilever may cause a closed-loop chaotic behavior in presences of nonlinear tip-sample interactions [18]. It has been previously reported that this chaotic behavior strictly limits the speed of the TM-AFM and can only be avoided by reducing the control gains and imaging speed [18].
To understand and improve the speed limit of AFM, an in-depth investigation of the transient response of the cantilever is crucial. However, the transient dynamics of the cantilever, which governs the overall performance of AFM, is often ignored or purely dis-cussed from a tip-sample interactions point of view [19]. Therefore, here we try to answer the following research question: how does the amplitude (and phase) of the motion of the cantilever evolve in transient situations? For this, we derive a set of governing differential equations based on Fourier components of the force and displacement. Different aspects of the proposed model are verified with experiments and compared with existing steady-state models. The proposed model is generic in a sense that it does not contain a certain tip-sample interaction model and enables graphical interpretations. A graphical interpretation is helpful for understanding the effects of cantilever dynamics on high-end aspects such as image quality and speed limit. Moreover, the proposed model can be a first step toward designing model-based controllers for high-speed nondestructive TM-AFM.
This paper is further organized as follows: In Sect. 2, we present a detailed derivation of the proposed mathematical model based on an averaging approach. In Sect. 3, we present analytic, numerical, and experimental results based on the proposed model. This is divided into four different subsections. The first subsection shows that the linear steady-state response of the proposed model is exactly equal to a steady-state response of a one-degree-of-freedom (DOF) resonator. The second subsection compares the nonlinear steadystate response of the proposed method with the existing theories for TM-AFM. In the third subsection, we introduce the transient response in time domain and verify it with experiments. Finally, in the fourth subsection, we use a specific experiment to verify the transient response in frequency domain and try to explain some of the existing experimental results. Section 4 is devoted to practical implications of the transient behavior for high-speed TM-AFM. This section consists of two subsections. The first subsection explains the effect of cantilever dynamics on image quality based on closed-loop nonlinear simulation. The second subsection graphically explains the reason behind the chaotic behavior of TM-AFM as a consequence of high control gains as previously reported in [18].

Mathematical modeling
According to experiments and previously achieved numerical results, despite the strong nonlinearities of the tip-sample interactions, the motion of the cantilever in conventional TM-AFM is harmonic up to a large extent [20]. The reason behind this harmonic motion is the extreme contrast between the sensitivity of the cantilever to different forces. The cantilever is highly sensitive to the forces that correspond to its resonance frequency, and almost not responsive to any other forces. Therefore, the Fourier component of the forces which correspond to the fundamental resonance frequency of the cantilever generates a measurable displacement, while the effects of the other forces are likely to be obscured by the noise [20]. When this contrast is low [low quality (Q)-factor [21] or multi-harmonic cantilevers [22,23]] or if the TSI force is so strong that it compensates for the low sensitivity, the motion of the cantilever can contain some higher frequency components. However, for conventional TM-AFM, and in the presence of measurement noises, it is almost impossible to measure the higher frequency content of the motion.
Considering only a single harmonic motion, it is a legitimate choice to model the AFM cantilever as a one-DOF mass-spring-damper system. For brevity, we start with a non-dimensional form of a one-DOF model of the AFM cantilever as: where x is the non-dimensional displacement of the tip, with over-dot representing the time derivative. ξ = 1 Q is the damping ratio, with Q being the Q-factor of the cantilever. F d , f ts and ω represent the normalized driving force, TSI force, and excitation frequency, respectively. Although Eq. (1) represents the dynamics of the cantilever up to a reasonably high precision, it is not very useful for control design purposes. As mentioned before, the motion of the cantilever only contains one harmonic component and the controller observes its amplitude, which by definition varies much slower than the motion of the cantilever. In practice, a lock-in amplifier (LIA) demodulates the motion of the cantilever into its amplitude and phase and feeds them to the controller. The details of the LIA are outside the scope of this paper, however, it is essential to note that there does not exist any transfer function or linear approximation for the LIA. Hence, to incorporate the functionality of the LIA, one should either solve for the strongly nonlinear multiple-time-scale dynamics of the cantilever which is coupled with the LIA, or derive a demodulated model, directly. Obviously, the first would not be efficient for applications such as control design or long time-horizon simulations. Moreover, the TSI force ( f ts ) in Eq. (1) is normally a continuous, but nondifferentiable function. Theoretically, such functions do not have any Taylor series approximation, which makes it even more elusive for stability analysis in control design. To overcome these problems, we opt for deriving a demodulated model for the cantilever.
Defining s 1 ẋ and s 2 x, Eq. (1) reads in statespace format: Since the amplitude and the phase evolve with a much slower time scale, it is useful to define a new time coordinate τ , assuming that the signals in τ domain are constant during one cycle of vibration of the cantilever. Thus, τ shall be used to describe any slowly varying function. A function f (τ ) is called "slow" in contrast to a rapidly varying periodic function g(t) with the periodicity of T , if it has the following property: Note that, τ by itself does not differ from t, it is only a notation which is introduced to explicitly separate the slowly varying functions from quickly changing ones. Assuming that each of the state variables are harmonic functions with slowly varying amplitude and phase (hereafter referred to as semi-harmonic), one can write: where i = 1, 2, j = √ −1, and indicates the real operator ( will be used for the Imaginary operator). Note that the static component and all the higher harmonics of the s signal is filtered by LIA, and bandwidth of the A and ϕ is limited to that of the LIA. By defining X i (τ ) A i (τ )e jϕ i (τ ) ∈ C, the state variables can be written as the following explicit multiplication of a slowly varying complex function and a pure harmonic function: Hence, the complex variable X i (τ ) represents the amplitude and phase of s i (t) corresponding to frequency ω. Using the chain rule of differentiation, and since differentiation is distributive for the operator, we have: In the same manner, the TSI force is a semi-periodic signal. Thus, it can be decomposed into its semi-harmonic components as: where F (n) ts (τ ) ∈ C represents the amplitude and phase of the n th harmonic component of the tip-sample interaction force. Note that Eq. (6) differs from the standard complex Fourier transform only in a way that the amplitude values are not necessarily constant, but represent slowly varying functions in time domain τ .
Substituting Eqs. (4) to (6) in Eq. (2) yields: Multiplying both sides of Eq. (7) with e jωt and integrating through a vibration cycle 2π ω 0 (·)e jωt dt , one can project the equations onto the space of the first harmonic component as: Notice that we deliberately expanded the zeroth and first harmonics of the TSI force (Γ 0 , and F (1) ts in Γ 1 ) out of the (sum) operator. In this way, the equations are rearranged in frequency order where Γ 0 , Γ 1 , and Γ n represent the terms with zero frequency, first harmonic, and higher harmonics, respectively.
Moreover, from the orthogonality of harmonic functions, it is easy to check that, for ∀ n ∈ N, c ∈ C we have: Considering the first statement of the orthogonality above, and the definition in Eq. (3), the last term in Eq. (8a) vanishes (Γ n = 0). Applying the second statement of orthogonality to the remaining terms in Eq. (8), the operator drops and the following differential equations are obtained: Equation (9) is an ordinary differential equation with complex variables and coefficients. By defining new state parameters as q and separating the real and imaginary parts of Eq. (9), the governing differential equations for the modulated system can be written in the standard real valued state-space form as: where is the slowly varying first Fourier component of the TSI force, A 2 (τ ) represents the amplitude of the motion, and ϕ 2 (τ ) is the phase.
It might be counterintuitive to observe that the second-order system in Eq. (1) is converted to a fourthorder system in Eq. (10). However, it is straightforward to check that Eq. (10) is a minimal realization, meaning that, four is the minimum dimensionality of the state-space that can accurately represent the modulated system.
It is useful to summarize the physical meaning of the equations and the derivation steps. We started with a single-DOF mass-spring model of the cantilever in Eqs. (1) and (2). Then, assumed that the amplitude and phase of the cantilever change significantly slower than the motion of the cantilever in Eqs. (3) and (4). Furthermore, by projecting the equations onto the Fourier kernel in Eq. (8), the differential equations for amplitude and phase at the slow timescale Eq. (10a) have been achieved. According to Eq. (10a), the amplitude modulated cantilever is a fourth-order linear dynamic system, with a nonlinear sensing (amplitude and phase). The dynamic properties of this system depend on its Qfactor and excitation frequency. The input of the modulated system only consists of the dither force and the first harmonic of the tip-sample interaction force.
The first harmonic of the force (F ts ) has a nonlinear relationship with the state variables (q i ) and the distance of the cantilever from the sample surface. "Appendix A" presents the derivation of F (1) ts as a function of the distance and the state parameters for the well-known Derjaguin-Muller-Toporov (DMT) force model.

Numerical and experimental results
In this section, the proposed formulation is used to study the dynamic behavior of AFM cantilevers. First, the steady-state response of the proposed model is verified with existing models. Next, the transient behavior of the cantilever is studied in time and frequency domains and verified with experiments.

Linear steady-state response
The most simple test to check the proposed model is to evaluate it for linear steady-state case. For this, the steady solution of Eq. (10) can be calculated by putting the TSI force and the left-hand side of Eq. (10a) equal to zero (F (1) The subscript ss stands for steady state. This algebraic set of equations can be solved analytically as follows.
The last two equations in Eq. (11) can be solved separately as: q 1 = −ωq 4 and q 2 = ωq 3 . Substituting this results in the first two equations of the Eq. (11) gives: Solving the above equation analytically, the stationary state variables become: Substituting the static solution for state variables into Eq. (10b), the following well-known relations are obtained: .
This shows that the static response of the proposed model in linear case indeed leads to the steady-state response of a one-DOF linear resonator [24].

Nonlinear steady-state response
The proposed model eliminates the fast timescale and represents the dynamics of AFM cantilever in slow timescale. Therefore, its static response corresponds to the steady-state response of the AFM cantilever, whereas its dynamic response represents the transient behavior. In this subsection, the steady response of the proposed model is compared to the existing theories presented in [12,17,25]. Similar to the previous subsection, by putting the left-hand side of Eq. (10a) equal to zero, the steady-state response of the AFM cantilever can be obtained. However, in contrast to the previous section, the TSI force is not equal to zero, but its relationship with state parameters has to be calculated via a so-called force model. For this, we used the Derjaguin-Muller-Toporov (DMT) model which consists of the attractive van der Waals (vdW), repulsive Hertz, and dissipative viscoelastic forces. For details on this model please refer to "Appendix A" and the references [26][27][28].
Using this model, the first harmonic component of the TSI force can be written as: , and γ = ωη R A 3 0 , are the coefficients of the vdW, Hertz and viscoelastic forces, respectively. H, R, E Eff , h, A 0 , A 2 , and σ are Hammaker constant, tip radius, effective stiffness of tip-sample contact, separation of the sample surface and the cantilever in its undeflected configuration, free air amplitude, actual amplitude at any time [as described in Eq. (10b)], and the intermolecular distance, respectively. The integral functions (I 1 , I 2 and I 3 ) as a function of their arguments (ζ 1 , ζ 2 ) are defined as: where the discontinuity function (a) D b is defined to impose the discontinuity of the forces during the con- . A detailed derivation of the force model for slow timescale and the details of DMT model is presented in "Appendix A". Typically, an arc-length continuation method is used to calculate the nonlinear frequency response of a nonlinear system. However, the frequency response of the system represented by Eq. (10) can always be transformed to a quadratic algebraic problem in terms of excitation frequency squared (ω 2 ), irrespective of Fig. 2 Nonlinear frequency response curves of cantilever in TM-AFM for different cantilever-sample separations; red represents the smallest and blue represent the largest distance. The greendashed curve schematically represents the amplitude-distance relationship as shown in Fig. 1. Simulation parameters are the same as reference [25], and the results in this figure are in good agreement with literature, for example see reference [25]. (Color figure online) the type of the nonlinearity. Hence, its nonlinear frequency response can also be calculated analytically as well. The details of calculating the nonlinear frequency response of the cantilever using the proposed model are presented in "Appendix B". Figure 2 shows the steady-state nonlinear frequency response of the cantilever considering different cantilever-sample separations. These results were calculated with the same parameters as in [25] (spring constant 2 N/m, resonance frequency 52.4 kHz, quality factor 66.7, tip radius 20nm Hammaker constant 2.96 × 10 − 19 J, intermolecular distance 2.0 Å, and free air amplitude 100 nm). The nonlinear frequency responses in Fig. 2 are in good agreement with the results presented by Lee et al. [25], whereas the dashed green curve represents the amplitude-distance relationship depicted in Fig. 1.
These results show that the static response of the proposed model agrees with the existing models for steadystate response of the AFM cantilever. The next sections will study the transient behavior of the AFM cantilever which captures the dynamic transition between each of the lines in Fig. 2.

Transient response in the time domain
To investigate the evolution of the amplitude and phase in the time domain, here the step response of the modulated system is studied. The assumption is that the cantilever is initially at rest (zero amplitude). Suddenly, a harmonic force with a constant amplitude and frequency is applied (i.e., dither piezo turned on). Fig-ure 3 shows the dynamic trajectory of the amplitude and phase from rest to its steady-state.
Note that the dynamic trajectory of the cantilever is a spiral of which the direction is determined by the ratio between the excitation frequency and the resonance frequency of the cantilever. If the excitation frequency is lower than the resonance frequency, the cantilever follows a counterclockwise spiral, and vice versa. This is illustrated via numerical and experimental results in Fig. 4. Here, the numerical results show the step response of linear part of Eq. (10) and experimental results are captured using a commercially available AFM and LIA , i.e., Bruker Fast-Scan and Zurich Instruments UHFLI 1.8 GSa/s 600 MHz. In the experiments, first the cantilever was settled to a certain amplitude and phase, then suddenly, we increase the excitation power to twice its steady value. Note that the initial point in experiments start from different phase values, with approximately identical amplitudes, because the excitation frequency was different for the three cases. Moreover, we did not start the experiment from rest condition, because the phase signal is defined based on the dither signal and cannot be defined while the cantilever is at rest.
The spiral shape of the trajectory is more relevant when high-speed AFM controllers are concerned. The spiral-shaped trajectory in Fig. 4 shows that there can be an initial response in the wrong direction, similar to the non-minimum-phase (MNP) behavior in linear systems [29]. Meaning that, if for any reason one expects the amplitude to drop, it might first increase and then drop, or vice versa. This wrong direction of the initial response makes it very challenging for high-speed con- Fig. 3 Step response of the amplitude and phase in polar and Cartesian coordinates. In the phasor plot, the distance from the center of coordinates shows the amplitude and the angle shows the phase delay. The reasoning behind overshoot of amplitude and the trajectory of phase is stem from the spiral trajectory shown in the phasor plot

Fig. 4
Dynamic trajectories of amplitude-phase pair with different excitation frequency. a Experimental: a cantilever with the Q factor ≈ 300, stiffness ≈ 1 N/m and fundamental resonance frequency ≈ 45 kHz) has been used in a commercial AFM system. While the cantilever was far from any sample surface, the excitation power has been suddenly increased from 3 to 6 V. b Numerical: step responses of Eq. (10). If the excitation frequency is less than the resonance frequency of the cantilever, the dynamic trajectory is a clockwise spiral and vice versa trollers. If the time constant of the controller is shorter than the settling time of the cantilever, the controller might take wrong actions based on the wrong directional response of the cantilever.

Nonlinear transient response in the frequency domain
This subsection investigates the frequency domain response of the demodulated system. In order to understand the model represented by Eq. (10) in frequency domain, it is assumed that the excitation force is modulated with another harmonic signal which has a much lower frequency than the resonance frequency of the cantilever. Considering the linear dynamics of the modulated system, i.e., Eq. (10a), one observes that the system has two pairs of complex conjugate poles [eigenvalues of matrix Λ in Eq. (10)] at: Normally, an AFM cantilever in ambient conditions is a highly underdamped system. Therefore one can simplify Eq. (16) further considering that ω ≈ 1 and ξ ≈ 0, and ξ 2 (1 − ω). 1 In that case, one observes that the system has one pair of dominant and one pair of non-dominant poles at: The non-dominant pair of poles, with complex part 2 (P 1,2 ), represent the up-modulation of the resonator and are not relevant for the TM-AFM problem because, besides being non-dominant, their response would be filtered by the LIA. The other two poles (P 3,4 ), which are dominant, actually determine the transient behavior of the AFM cantilever. Eq. (17) suggests that the imaginary part of the dominant poles of the modulated system (P 3,4 ) are approximately equal to the difference between the excitation frequency and the resonance frequency of the cantilever, i.e., |ω e − ω r | (or in non-dimensional form |ω − 1|), and their real part ( ξ 2 = 1 2Q ) represent the settling time of the system. In fact, the real part of the dominant poles which determine the relaxation time of the system is already known by AFM experts as it is a common knowledge that lower the Q-factor higher the speed of AFM. However, the imaginary part is not well understood.
Recently, the relationship between the real part of the poles (quality factor of the cantilever) has been experimentally confirmed by Adams et al. [30] They used different custom-made cantilevers with various Q-factors and showed the dependency of the settling time to the Q-factor. In their experiments, they used a piezoelectric actuator to modulate the surface height, artificially. In this way they measured the frequency response of the amplitude due to height variations. A very interesting and counterintuitive detail in their experiments is that the magnitude of the amplitude signal could go beyond zero decibel at certain frequencies. This effect which was not explained, is an evidence of the resonance in the modulated system (see Figure 4.b in [30]). Note that the presence of a resonance can only be explained with complex conjugate poles.
To systematically investigate the resonance of the amplitude signal, we repeated the experiment presented by Adams et al. [30]. However, with the difference that, instead of modulating the surface, we modulated the amplitude of the excitation signal. In this way the nonlinearities regarding the TSI force are out of the equation, and the pure dynamics of the cantilever is measured. For this experiment, two lock-in amplifiers were used as shown in Fig. 5. One of the LIAs provides the modulation/demodulation of the cantilever with the socalled carrier frequency, i.e., the excitation frequency of the dither (ω e ), whereas the second LIA modulates the amplitude of the carrier signal and extracts the amplitude and phase of the amplitude of the cantilever. In this way, effects of the fluctuation of harmonic force are measured depending on the frequency content of the fluctuation.  Figure 6 shows both the measured and calculated frequency response of the amplitude to these fluctuations. The numerical results were achieved as follows: The transfer function matrix for the linear part of the system [from F(ω M ) to q(ω M )] is calculated as the first column of (s I − Λ) −1 where s is the Laplace variable and I is a (4 × 4) unit matrix. Therefore, substituting a probe function as F d = sin(ω M τ ), the amplitude and phase signals (A 2 and ϕ 2 ) can be calculated analytically. Then, the functionality of the second LIA can be implemented by multiplying the amplitude signal ( A 2 ) by the probe signal and integrating as: In this way amplitude and phase (A A , ϕ A ) of the amplitude signal ( A 2 ) due to disturbance in total harmonic force (F d ) were calculated.
In Fig. 6, five different carrier frequencies were chosen around the resonance frequency with an interval of 500 Hz. As it can be seen, there is a good agreement between the experimental and numerical results. Both experimental and numerical results demonstrate that the peak of the amplitude occurs at a frequency equal to the difference between the carrier and the resonance frequencies (|ω e −ω r |), which confirms Eq. (16). The experiments were done using a commercial AFM and LIAs (Bruker Fast-Scan and Zurich Instruments UHFLI 1.8 GSa/s 600 MHz) with a standard cantilever which has a resonance frequency of ω r = 319.015 kHz and Q-factor Q ≈ 520. Simulations consider the same cantilever.
An other important observation from the frequency domain analysis is the slope of the decay line after resonance. The linear part of Eq. (10) suggests that the system has one pair of dominant complex conjugate poles. Hence, the system should behave like a secondorder system, and have a decay line with a slope of −40 dB per decade after the resonance. However, due to the output nonlinearity [see Eq. (10b)], the slope of the decay line can differ from −20 till −40 dB per decade, depending on the carrier frequency. This observation suggests that it would be very challenging for any system identification method to find any reliable integer-order fit for the system. This hinders the design of model-based controllers. Obviously, this problem would not appear without the output nonlinearity related to A and ϕ. Therefore, to design a modelbased controller, one has to either design fractionalorder controllers [31,32] or use the q 3 and q 4 as the control input instead of the amplitude.

Implications of transient behavior of cantilevers in TM-AFM experiments
In this section, the nonlinear closed-loop behavior of the system is studied. The results show the practical implications of the transient behavior of the TM-AFM on the final image quality. Also the origin for previously reported chaotic behavior of the TM-AFM [18] can be explained.

Nonlinear closed-loop behavior of the TM-AFM
To investigate the nonlinear behavior of AFM in a closed-loop setting, a force model and a model for the controller are needed, besides the model of the cantilever. For the force, we use the slow time domain DMT model as presented in previous subsections (Eq. 14). A detailed derivation of the slow time model is presented in "Appendix A". As the controller, an ideal proportional-integral (PI) controller is considered which assumes that the z-stage actuator is fast enough not to have any effect on the closed-loop dynamics of the system. In this way the distance between the cantilever and the sample (h) is defined as: where k p , k i , and A set are proportional gain, integral gain and set-point amplitude, respectively. Equations (10), (14) and (19) representing the cantilever (and LIA), the tip-sample interactions, and the controller are coupled with each other through shared signals A, h, and F (1) ts . Figure 7 shows the simulated amplitude, phase, and the magnitude of the harmonic component of the TSI forces during an approach process.
As it can be seen from Fig. 7, every time that the cantilever engages the surface (i.e., the distance becomes equal to the amplitude of the cantilever) the TSI force emerges and affects the cantilever and changes its amplitude and phase. It is important to notice that dur-  Height image of a calibration sample (UMG02B from Anfatec Instruments) with three different excitation frequencies, corresponding to Fig. 8. All three images are taken from the same spot on the sample, using a standard tapping mode cantilever with a resonance frequency of 72 kHz, Q factor of 200, and spring constant of 2.5 N/m. The free air amplitude, set-point amplitude, and scan rate for all three cases are 100 nm, 75 nm, and 1 Hz, respectively. The color masks are applied to exaggerate the imaging artifacts. Imaging with a lower excitation frequency adds fewer artifacts ing the engagement time, the TSI force causes always a considerable phase lead, while the amplitude mildly changes. However, after disengaging the surface, the amplitude value keeps reducing with a larger rate and then increases again. To further investigate the transient behavior of the amplitude and the phase, Fig. 8 shows the dynamic trajectory of the system in phasor plane for three different excitation frequencies. Using these graphs, we can study the following three questions: (i) Why do the tip-sample interactions always induce a phase lead? (ii) Why does the amplitude keep reducing and then increasing when disengaged from the surface? (iii) Why is the image quality better when we choose an excitation frequency slightly lower than the resonance frequency of the cantilever?
Before the cantilever engages to the surface, the motion of the cantilever has a phase delay with respect to the dither force. Since the TSI force is opposite to the displacement, the direction of the total harmonic force at the first moment of the contact is upwards (a positive phase). Also because the TSI force in first engagement is much stronger than the dither force, the total harmonic force is mainly dominated by the TSI force. In this situation, the phase increases and amplitude tries to increase by indenting the tip more and more into the sample surface. Hence, the cantilever follows  Fig. 9. The profiles measured with an excitation frequency less than the resonance frequency of the cantilever show a less noisy, more stable profile the circular arc in Fig. 8 in counterclockwise. At the moment that the phase lead does not need to increase anymore, the cantilever looses contact with the surface and only the dither force is acting on the cantilever. After losing the connection, the cantilever follows its free trajectory which first reduces the amplitude and then increases. As explained in Fig. 4, the free trajectory of the motion can be either a straight line or a spiral trajectory, depending on the excitation frequency. Repeating this engaging-disengaging process for few times, the cantilever reaches its steady-state situation. The damping needed for this process is provided by all three elements: the cantilever, the non-conservative part of TSI force, and the control action.
Comparing the three different cases in Fig. 8, for the lower excitation frequency, the cantilever has a lower phase delay before the contact. Thus, it needs to follow a shorter trajectory on the circular arc to lose contact. On the other hand, after losing contact, it follows a counterclockwise spiral to reach the surface again which is shorter than a clockwise trajectory to the surface. Therefore, in total the cantilever has a shorter path to the stable steady-state situation, when it is excited with a lower frequency. That is why the images captured with a lower excitation frequency have higher image quality, which is intuitively known by experienced AFM users and is experimentally demonstrated in Fig. 9. For this example, a calibration sample is measured with three different excitation frequencies where all other parameters are exactly the same for all the three cases. The difference in image quality can be observed more clearly in Fig. 10 which illustrates the added imaging artifacts on the cross section of one of the holes in Fig. 9.

Chaotic behavior
It has been reported that if the controller is tuned to be faster than a certain threshold, the closed-loop system shows a chaotic behavior [18]. Although the presence of chaos was confirmed by studying the Poincaré sections and Lyapunov exponents, the origin of the chaos was not explained. Figure 11 shows the difference between a stable and a chaotic trajectory, and Fig. 12 shows the images captured with conditions corresponding to Fig. 11. In the scenario as explained in the previous subsection, if the controller acts more aggressive than a certain limit, the cantilever experiences more than 2π radian phase lead before losing the surface and re-engaging. Therefore, the next attempt to re-engage does not start from a better initial point. Thus, the new initial condition is not closer to the stable point comparing the previous engaging point. Hence, Fig. 11 Closed-loop behavior of TM-AFM in polar coordinates. a Chaotic. b Stable. The aggressive control action causes more than 2π radians phase lead in each engage-disengage cycle and prevents the system from reaching a stable center point which generates a chaotic motion Fig. 12 Height image of the calibration sample (UMG02B from Anfatec Instruments) with parameters corresponding to Fig. 11. When the cantilever is in chaotic regime (the control gains are too high), the images do not provide any useful information the "engage/disengage" process repeats forever in a periodic or non-periodic manner, depending on the control gains. While the periodic "engage/disengage" generates a quasi-periodic regime, the non-periodic one represents the chaotic trajectory. Both these regimes have been reported in [18]. All in all, the chaotic behavior can be attributed to the wrong direction response of the amplitude signal.

Conclusions
In this paper, we presented a dynamic model for amplitude modulation (tapping mode) AFM as the first step toward a model-based control design. The model graphically explains the behavior of the AFM cantilever in a slow timescale, i.e., the changes in the amplitude, phase and the control signals. The proposed model has been verified with experiments and shows that the behavior of the AFM cantilever in slow timescale is profoundly affected by the excitation frequency and Q-factor of the cantilever. According to the presented model the amplitude per sé is not the best indication of the distance, and should not be used as the error signal in the control loop. Instead, to design highperformance controllers and avoid chaos, one should consider a modulated transient model of the cantilever as a multi-input-multi-output system. tion (a) D b is defined to impose the discontinuity of the forces during the contact as: Normalizing the model in Eq. (20) according to the same scales as in Eq. (1), the non-dimensional form of the TSI force f ts = f ph ts k A 0 can be written as: , and γ = ωη R A 3 0 , are the coefficients of the vdW, Hertz and viscoelastic forces, respectively. A 0 , and k represent the free air amplitude and the spring constant of the cantilever. Applying the Fourier operation as explained in Sect. 2 , we obtain: where the first harmonic of the van der Waals force (F  Considering the definition of X i = A i e ϕ i , i = 1, 2, and defining variables θ i = ωt + ϕ i , i = 1, 2, ζ 1 = h A 2 , and ζ 2 = σ A 2 , the integral in Eq. (23a) can be simplified as: Repeating the same procedure for Eqs. (23b) and (23c), we obtain: where K vdW , K H , and C vis are defined as: and can be considered as added negative stiffness due to the vdW force, added stiffness by the Hertz force, and the added damping of the viscoelastic force, respectively. The integral functions in Eq. (25) are defined as: and are illustrated in Fig. 13. Note that, the conservative part of the force ((K H − K vdW )X 2 ) is in the opposite Fig. 13 Integral functions I 1 till I 3 as a function their arguments direction of the modulated displacement (X 2 ), and the dissipative force (C vis X 1 ) is in the opposite direction of the modulated tip velocity (X 1 ). Therefore,(K H − K vdW ), and C vis can be considered as the added stiffness and the damping of the TSI force, respectively.

Appendix B: Calculation of the nonlinear frequency response curve
The steady-state frequency response of the cantilever can be calculated by substituting Eq. (24) into Eq. (10) and puttingq i = 0 yields: ⎡ As it can be seen from Eq. (25), the added nonlinear stiffness and damping values do not depend on the phase or frequency, instead they are only a function of amplitude. It is easy to check that this is the case for any time-invariant nonlinearity.
To avoid solving the nonlinear equation (27) for amplitude and phase, we suggest instead to scan the amplitude, and solve for frequency. Therefore, we defined Ξ and Υ in Eq. (27) such that Ξ X 2 = Υ . Considering Eq. (10b) (q 2 3 + q 2 4 = A 2 ), one can eliminate the phase by multiplying Eq. (27) by Ξ −1 , and left-multiply the transpose of the resultant by itself: Independent of the type of the nonlinearity (only if the nonlinearity is time independent), Ξ can always be written as a summation of a skew symmetric matrix and a scaled unit matrix. Therefore, ΞΞ T will always be diagonal, with both of the elements equal to each other which leads to the single algebraic equation: For one-DOF resonators, Eq. (29) is always quadratic in terms of ω 2 and has an analytic solution for any value of the amplitude. The square root of the real and positive solutions of the Eq. (29) show the frequency for each amplitude. After calculating the amplitudefrequency pair, the phase can be calculated simply by solving the Eq. (27) for X 2 , which is linear. The main advantage of this method as compared to arc-length methods is that in this way the computational time can be significantly reduced without compromising the accuracy. However, this method, in its present form, can only deal with one-DOF resonators with nonlinearities that are not explicitly time dependent.