Higher modes and higher harmonics in the non-contact atomic force microscopy

Higher modes and higher harmonics in the non-contact atomic force microscopy

International Journal of Non-Linear Mechanics 110 (2019) 33–43 Contents lists available at ScienceDirect International Journal of Non-Linear Mechani...

806KB Sizes 0 Downloads 23 Views

International Journal of Non-Linear Mechanics 110 (2019) 33–43

Contents lists available at ScienceDirect

International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm

Higher modes and higher harmonics in the non-contact atomic force microscopy Mohammad Safikhani Mahmoudi, Arash Ebrahimian, Arash Bahrami ∗ School of Mechanical Engineering, College of Engineering, University of Tehran, Tehran, Iran

ARTICLE

INFO

Keywords: Non-contact atomic force microscopy Higher modes and higher harmonics Direct harmonic balance Nonlinear oscillations

ABSTRACT The present paper proposes a new mathematical approach to study the role of higher harmonics and higher mode shapes in the dynamics of an atomic force microscope microcantilever. The basic principle of the present method, called the direct harmonic balance, is similar to the well-known method of harmonic balance, but the key difference is that this method does not require discretizing the partial differential equation, governing the probe dynamics. In contrast to the traditional harmonic balance, the direct harmonic balance method, introduced here, attacks the partial differential equation directly, leading to a more efficient solution procedure. The accuracy of this method is demonstrated using a perturbation method. The direct harmonic balance method is then employed to investigate the nonlinear oscillations of the atomic force microscope probe in the non-contact mode. The significance of higher mode shapes and higher harmonics is illustrated through several frequency-response curves, presented for various cases. It is found out that while higher modes are more sensitive to the amplitude of excitation, they have a smaller frequency shift. It is also found out that the second harmonic exhibits high sensitivity to the initial tip/sample separation distance. In addition, it is shown that reducing the tip/sample initial distance results in a considerable growth in the magnitude of the second harmonic, which can be used as a reliable signal to explore topography of the sample.

1. Introduction Since the revolutionary idea of atomic force microscopy (also called AFM) was introduced by Binning et al. [1] in mid-eighties, it has evolved into a standard method for material characterization and one of the most powerful techniques in nanoscale imaging. The most popular operational mode of the AFM is the dynamic mode. In this mode (also called dAFM), a vibrating probe is used to scan topography of the surface [2]. High-resolution images of DNA and live samples [3–5], nanometer-scale topography of polymers [6,7] and also true-atomicresolution images of semiconductors [8,9] have been obtained using dynamic AFM, so far. A variety of signals may be used to connect the probe dynamic behavior to surface properties. Amplitude of oscillation, frequency shift and phase lag between the excitation and response are main signals through which data is gathered [10,11]. Dynamic atomic force microscopy has undergone significant developments during the past three decades. Image resolution and contrast, ability to gather new kind of information about the sample, and compatibility with various environments have been remarkably improved in recent years. However, there is still room for improvement. Garcia et al. [12,13] have proposed a method to enhance image contrast without tip entrance to the contact region. They excited the first two cantilever eigenmodes simultaneously. Lozano and Garcia have shown coupling

between the modes and sensitivity of the second mode to long-range attractive forces allow to probe compositional changes of the sample more accurately [14]. During the past decade, exciting a single higher vibrating mode [15,16] or simultaneous excitation of two or more mode shapes have been attracted much attention [17–20]. In addition, it is well established by now that analyzing higher harmonics and frequencies is useful even if only one of the cantilever frequencies is excited [21–25]. It has been shown that use of higher eigenmodes are beneficial to avoid the jump-to-contact phenomenon [26,27]. High-resolution images with improved contrast of living bacterium [28] and living cells [29] have been obtained using higher harmonics. As it is rather hard to measure higher harmonics in air, it has been suggested to use special cantilevers which allows to tune higher harmonics with a higher eigenmode [30]. Furthermore, some researchers have suggested to employ higher torsional harmonics since they are easier to detect [31,32]. As the dynamic AFM techniques have evolved in recent years, mathematical models to simulate AFM complex dynamics have been developed too. Mathematical models are employed to predict behavior exhibited by a real AFM and understand the underlying physics. However, they usually fail to provide a quantitative agreement with the experiment. The discrepancy between theory and experiment is

∗ Corresponding author. E-mail address: [email protected] (A. Bahrami).

https://doi.org/10.1016/j.ijnonlinmec.2019.01.006 Received 2 December 2018; Received in revised form 11 January 2019; Accepted 12 January 2019 Available online 21 January 2019 0020-7462/© 2019 Elsevier Ltd. All rights reserved.

M.S. Mahmoudi, A. Ebrahimian and A. Bahrami

International Journal of Non-Linear Mechanics 110 (2019) 33–43

Fig. 1. Schematic of the AFM microcantilever.

2. Mathematical modeling

caused by different factors such as difficulty of considering all forces and parameters and/or using inefficient mathematical methods. Early works documented in the literature, have used point-mass model to explore the dynamic behavior of AFM [33–36]. This method can describe some nonlinear aspects of the tip motion but ignores distributed nature of the AFM microcantilever. In fact, the lumped-parameter approach is incapable of modeling the microcantilever higher vibration modes and also is not accurate enough [37–39]. Considering continuous nature of the AFM structure in modeling and analysis can provide a deeper insight into physics of the AFM [40]. Such models can capture all the nonlinear dynamics of the tip motion such as chaos, presence of higher harmonics and eigenmode in the response, and nonlinear coupling between mode shapes. Point-mass models result in a single nonlinear ODE that describes the probe motion. On the contrary, in continuous models the governing equation is a nonlinear PDE. Since working with a nonlinear PDE is, generally, more difficult compared to a nonlinear ODE, little attention has been paid to these models. In the literature, main approach to solve continuous-model equations is discretizing the PDE to spatial and temporal parts through the singlemode Galerkin procedure and then applying a numerical integration scheme to solve temporal equations. Bahrami and Nayfeh [41,42] have shown even while an AFM is operating in the air a single-mode discretization may lead to quantitative and/or qualitative errors. They showed precise capturing of all nonlinear phenomena is just possible after involving four or even more eigenmodes in the Galerkin method. On the other hand, only a few articles have tried to solve the governing equations with semi-analytical methods without any discretization of the PDE. Turner [43] has approximated the response of a contact AFM to the primary resonance by applying the method of multiple scales directly to the partial differential equation. He did not discretize the PDE, governing the AFM motion. Abdel-Rahman and Nayfeh [44] have used the same method to obtain the response of a contact AFM to a subharmonic excitation.

An atomic force microscope is basically composed of a microcantilever with a nanoscale tip at its free end. The microcantilever is mounted on a vibratory base and is excited through a piezoelectric patch to oscillate over the sample. Approximating the three-dimensional AFM structure as an Euler–Bernoulli beam and assuming strains are small, the governing equation of the AFM microcantilever and associated boundary conditions are derived. The AFM tip is assumed to be of negligible mass. Fig. 1 illustrates geometry of the AFM microcantilever and the sample. l and EI are, respectively, microcantilever length and bending rigidity. 𝑍̂ indicates initial tip/sample separation distance ̂ 𝑥, in the reference configuration. 𝑤( ̂ 𝑡̂) demonstrates deflection of the microcantilever relative to a frame attached to the base and 𝑦( ̂ 𝑡̂) signifies harmonic base excitation 𝑦( ̂ 𝑡̂) = 𝑌̂ cos(𝛺̂ 𝑡̂). Also, 𝑧( ̂ 𝑡̂) denotes the instantaneous tip/sample separation distance, determined as 𝑧( ̂ 𝑡̂) = 𝑍̂ − 𝑢(𝑙, ̂ 𝑡̂) where 𝑢( ̂ 𝑥, ̂ 𝑡̂) is absolute deflection of the beam. Actually, 𝑢( ̂ 𝑥, ̂ 𝑡̂) can be written as: ̂ 𝑥, 𝑢( ̂ 𝑥, ̂ 𝑡̂) = 𝑦( ̂ 𝑡̂) + 𝑤( ̂ 𝑡̂)

(1)

Furthermore, 𝑥̂ and 𝑡̂ represent axial coordinate and time, respectively. A variety of formulations has been used so far to model the interaction force between the sample and the tip. Most of these formulations have been developed by assuming a spherical apex for the tip and a flat surface for the sample. For such a configuration, nonlinear attractive van der Waals force may be approximated by: 𝐹𝑣𝑑𝑊 (𝑡̂) = −

𝐻𝑅 6𝑧̂ 2 (𝑡̂)

(2)

In Eq. (2), 𝐹𝑣𝑑𝑊 (𝑡) represents the van der Waals non-contact force. Also, 𝐻 and 𝑅 are the Hamaker constant and the tip radius, respectively. 𝐹𝑣𝑑𝑊 (𝑡) is a conservative force, thus its potential energy can be defined as: 𝑢̂

𝑉𝑣𝑑𝑊 (𝑢) ̂ =

In the present work, we propose a new approach to solve the initial– boundary value problem, governing the AFM probe oscillations. In the proposed method, which is referred to as the direct harmonic balance method, the partial differential equation is attacked directly and without any reduction to ordinary differential equations. One of the main advantages of the present method is that it does not require any discretization. Moreover, it is not limited to the AFM problem and could be applied to a variety of nonlinear initial–boundary value problems, faced in vibrations as well as other fields of engineering and physics. Our focus is on the higher mode shapes and higher harmonics and the possible role that these signals can play to improve the AFM performance. In what follows, first the direct harmonic balance method is introduced and its corresponding numerical results are compared with those of the method of multiple scales. Then, the direct harmonic balance method is employed to investigate the behavior of the higher modes and higher harmonics in the non-contact atomic force microscopy.

∫0

𝑢̂

−𝐹𝑣𝑑𝑊 𝑑𝑢 =

∫0

−𝐻𝑅 6(𝑍̂ − 𝑢) ̂

2

𝑑 𝑢̂ =

−𝐻𝑅 6(𝑍̂ − 𝑢) ̂

(3)

2.1. Equations of motion To study influence of the higher eigenmodes and higher harmonics in the AFM nonlinear oscillations, partial differential equation of the microcantilever is considered. As aforementioned, the cantilever governing equation is derived based on Euler–Bernoulli beam theory. Total kinetic energy and total potential energy of the system are, respectively, given by Eqs. (4) and (5). Also, the work done by the air viscous damping is presented by Eq. (6): 𝑙

𝑇 =

1 2 𝜌𝐴(𝑢( ̂̇ 𝑥, ̂ 𝑡̂)) 𝑑 𝑥̂ 2 ∫0

𝑉 =

1 𝜕 2 𝑢̂ 𝐻𝑅 𝐸𝐼( ) 𝑑 𝑥− ̂ ∫ 2 0 𝜕 𝑥̂ 2 6(𝑍̂ − 𝑢( ̂ 𝑥, ̂ 𝑡̂))

(4)

2

𝑙

(5)

𝑙

𝛿𝑊 = − 34

∫0

̂̇ 𝑥, ̂ 𝑥, 𝑐( ̂ 𝑤( ̂ 𝑡̂) + 𝑦( ̂̇ 𝑡̂))𝛿 𝑤( ̂ 𝑡̂)𝑑 𝑥̂

(6)

M.S. Mahmoudi, A. Ebrahimian and A. Bahrami

International Journal of Non-Linear Mechanics 110 (2019) 33–43

Table 1 Values of the AFM parameters. Symbol

Properties

Value

Unit

𝑅 𝑙 𝐴 𝐼 𝜌 𝐸 𝐻 𝑐̂

Tip radius Length Cross-section area Second moment of area Density Young’s modulus Hamaker constant Viscous damping per unit length

10 225 7.2𝑒−11 3.57𝑒−23 2300 170 2.96𝑒−19 2𝑒−4

nm µm m2 m4 kg∕m3 GPa J kg∕ms

where 𝑐̂ is the damping coefficient per unit length of the microcantilever. Using the Hamilton extended principle, equation of motion and associated boundary conditions are determined: 𝐸𝐼 𝑢̂ 𝑖𝑣 (𝑥, ̂ 𝑡̂) + 𝜌𝐴𝑢( ̂̈ 𝑥, ̂ 𝑡̂) + 𝑐̂𝑢( ̂̇ 𝑥, ̂ 𝑡̂) = 0 𝑢(0, ̂ 𝑡̂) = 𝑦( ̂ 𝑡̂), 𝑢̂ ′ (0, 𝑡̂) = 0, 𝑢̂ ′′ (𝑙, 𝑡̂) = 0 𝐻𝑅 𝐸𝐼 𝑢̂ ′′′ (𝑙, 𝑡̂) = − 2 ̂ 6(𝑍 − 𝑢(𝑙, ̂ 𝑡̂))

(7a) (7b) (7c)

Fig. 2. Variation of the tip static deflection with the initial tip/sample separation.

A prime indicates derivative with respect to 𝑥̂ and a dot represents derivative with respect to 𝑡̂.

3. Direct harmonic balance method The method of harmonic balance is one of the most powerful techniques among the frequency-domain methods to solve nonlinear ordinary differential equations and obtain the steady-state periodic response of a system. The conventional method of harmonic balance is used to solve nonlinear ordinary differential equations. In the present study, we modify this method to solve partial differential equations directly and without discretization.

2.2. Static deflection Cantilever’s static deflection is a result of nonlinear force acting on the tip. Ratio of the van der Waals force to stiffness of the cantilever can be found by measuring deflection of the tip. Understanding the static behavior of the AFM probe provides us with useful information about the device operation in different regimes and also sheds light on the dynamic problem. Neglecting terms including time derivative and excitation in Eqs. (7a)–(7c), the static equation and associated boundary conditions are determined: 𝑖𝑣

𝐸𝐼 𝑢̂⋆ (𝑥) ̂ =0 ′ ′′ 𝑢̂ ⋆ (0) = 𝑦(𝑡), 𝑢̂⋆ (0) = 0, 𝑢̂⋆ (𝑙) = 0 ′′′ 𝐻𝑅 ⋆ ̂ 𝐸𝐼 𝑢 (𝑙) = − 2 6(𝑍̂ − 𝑢̂⋆ (𝑙) )

3.1. Solution procedure In the standard form of the method of harmonic balance, the response is assumed to be a linear sum of harmonic terms with frequencies that are integer multiples of the excitation frequency (for more details see Nayfeh’s nonlinear Oscillations book [45]). In the present research, however, we seek a solution for a partial differential equation. Therefore, we modify the method of harmonic balance and assume that the microcantilever deflection can be expanded in the following form:

(8a) (8b) (8c)

𝑢̂⋆ (𝑥) ̂ denotes static deflection of the cantilever and is a function of tip/sample initial distance 𝑍̂ as well as material and geometrical properties. ▵ represents tip/sample static separation distance, defined as: ▵= 𝑍̂ − 𝑢̂⋆ (𝑙)

𝑢(𝑥, 𝑡) = 𝜑0 (𝑥) + 𝜑1 (𝑥) cos(𝛺𝑡) + 𝜓1 (𝑥) sin(𝛺𝑡) + 𝜑2 (𝑥) cos(2𝛺𝑡) +𝜓2 (𝑥) sin(2𝛺𝑡) + ⋯

(9)

(12)

where 𝜑0 (𝑥), 𝜑𝑖 (𝑥) and 𝜓𝑖 (𝑥) are unknown functions which are found later by substituting 𝑢(𝑥, 𝑡) given by Eq. (12) into the PDE and associated boundary conditions. We note that 𝛺 is the nondimensional excitation frequency. We also note that the number of harmonic terms in Eq. (12) may be increased to satisfy the desired convergence criterion. To demonstrate how this method, which is called the direct harmonic balance method, works we consider a two-term harmonic expansion. Next, we substitute the two-term harmonic expansion into Eq. (11a) and set coefficients of sine, cosine and constant terms equal to zero to obtain: Time-independent term:

Numerical values of the AFM parameters are listed in Table 1 [14]. Fig. 2 represents variation of static separation distance with initial separation gap. Solving Eqs. (8a)–(8c), yields two branches of solution for ▵. As observed in the figure, the static response includes both stable and unstable solutions. The stability of each branch is determined by analyzing the potential function of the system. The static bifurcation diagram illustrated in Fig. 2 is in qualitative agreement with that of Bahrami and Nayfeh [41]. We note that the AFM tip does not touch the sample in the non-contact operation of the AFM. Results of static deflection are used to make parameters nondimensional. Dimensionless parameters are defined as: √ √ 𝑝𝐴𝑙4 𝑥̂ 𝑢̂ 𝐸𝐼 𝐻𝑅𝑙3 ̂ ̂ 𝑥 = ,u = ,t = t ,𝛺 = 𝛺 ,𝐻 = 𝑙 𝛥 𝐸𝐼 𝑝𝐴𝑙4 6𝐸𝐼𝛥3 (10) 𝑢̂ 𝑦̂ 𝑌̂ 𝑍̂ 𝑢̂ 𝑐𝑙 ̂2 , 𝑦 = , 𝑌 = , 𝑍 = , 𝑢 = , u⋆ = c= √ 𝛥 𝛥 𝛥 𝛥 𝛥 𝐸𝐼𝜌𝐴

𝜑0 𝑖𝑣 (𝑥) = 0

(13)

Coefficient of cos(𝛺𝑡) 𝜑1 𝑖𝑣 (𝑥) − 𝛺2 𝜑1 (𝑥) + 𝑐𝛺𝜓1 (𝑥) = 0

(14)

Coefficient of sin(𝛺𝑡)

Substituting Eq. (10) into Eqs. (8a) through (8c), yields nondimensional governing equation and associated boundary conditions as:

𝜓1 𝑖𝑣 (𝑥) − 𝛺2 𝜓1 (𝑥) − 𝑐𝛺𝜑1 (𝑥) = 0

(15)

Coefficient of cos(2𝛺𝑡) 𝑢𝑖𝑣 (𝑥, 𝑡) + 𝑢(𝑥, ̈ 𝑡) + 𝑐 𝑢(𝑥, ̇ 𝑡) = 0 𝑢(0, 𝑡) = 0, u′ (0, 𝑡) = 0, u′′ (1, 𝑡) = 0 𝐻 𝑢′′′ (1, 𝑡) = − (𝑍 − 𝑢(1, 𝑡))2

𝜑2 𝑖𝑣 (𝑥) − 4𝛺2 𝜑2 + 2𝛺𝑐𝜓2 = 0

(11a) (11b)

(16)

Coefficient of sin(2𝛺𝑡)

(11c)

𝜓2 𝑖𝑣 (𝑥) − 4𝛺2 𝜓2 − 2𝛺𝑐𝜑2 = 0 35

(17)

M.S. Mahmoudi, A. Ebrahimian and A. Bahrami

International Journal of Non-Linear Mechanics 110 (2019) 33–43

By substitution of the two-term harmonic expansion into the boundary conditions given by Eqs. (11b) and (11c), separating coefficients of linearly independent terms, and setting them equal to zero, the following equations are found (we note that the nonlinear boundary condition given by Eq. (11c) is expanded after the substitution): Time-independent terms: 𝜑0 (0) = 0, 𝜑′0 (0) = 0, 𝜑′′ (1) = 0 0 ( 𝐹0 𝜑0 (1), 𝜑1 (1), 𝜑2 (1), 𝜓1 (1), 𝜓2 (1), 𝜑′′′ 0 (1), 𝜑′′′ 1 (1), ) 𝜑′′′ 2 (1), 𝜓 ′′′ 1 (1), 𝜓 ′′′ 2 (1), 𝑌 , 𝐻, 𝑍 = 0

Overall, we have 36 unknown coefficients including 32 constants appearing in Eq. (27) and those showing up in Eq. (23). Substituting Eqs. (23) and (27) into the boundary conditions given by Eqs. (18)– (22) leads to a system of 20 algebraic equations. This system consists of 15 linear equations corresponding to 3 linear boundary conditions and 5 nonlinear equations associated with the nonlinear boundary condition. Since we have 36 unknowns and 20 equations, we need 16 more independent equations to find the unknown coefficients. To this end, we substitute values obtained for 𝜆1 and 𝜆2 back into Eq. (25) and obtain 16 independent equations. We note that the present method may be applied to solve similar nonlinear initial–boundary value problems, especially problems with linear partial differential equations and nonlinear boundary conditions. The proposed method is well-suited to treat such problems, arising frequently in dynamics of microelectromechanical systems.

(18)

Coefficients of cos(𝛺𝑡): 𝜑1 (0) = 0, 𝜑′1 (0) = 0, 𝜑′′ (1) = 0 1 ( 𝐹1 𝜑0 (1), 𝜑1 (1), 𝜑2 (1), 𝜓1 (1), 𝜓2 (1), 𝜑′′′ 0 (1), 𝜑′′′ 1 (1), ) 𝜑′′′ 2 (1), 𝜓 ′′′ 1 (1), 𝜓 ′′′ 2 (1), 𝑌 , 𝑍 = 0

(19)

3.2. Verification using perturbation theory

Coefficients of sin(𝛺𝑡): 𝜓1 (0) = 0, 𝜓1′ (0) = 0, 𝜓1′′ (1) = 0 ( 𝐺1 𝜑0 (1), 𝜑1 (1), 𝜑2 (1), 𝜓1 (1), 𝜓2 (1), 𝜑′′′ 0 (1), 𝜑′′′ 1 (1), 𝜑′′′ 2 (1), ) 𝜓 ′′′ 1 (1), 𝜓 ′′′ 2 (1), 𝑌 , 𝑍 = 0

In this section, perturbation theory is utilized to verify the direct harmonic balance method. To this end, the method of multiple scales is used to solve governing equations of the system presented in Eqs. (11a)– (11c). To implement the method of multiple scales, the equations are first written around cantilever equilibrium position 𝑢⋆ :

(20)

Coefficients of cos(2𝛺𝑡): 𝜑2 (0) = 0, 𝜑′2 (0) = 0, 𝜑′′ (1) = 0 2 ( 𝐹2 𝜑0 (1), 𝜑1 (1), 𝜑2 (1), 𝜓1 (1), 𝜓2 (1), 𝜑′′′ 0 (1), 𝜑′′′ 1 (1), 𝜑′′′ 2 (1), ) 𝜓 ′′′ 1 (1), 𝜓 ′′′ 2 (1), 𝑌 , 𝑍 = 0

𝑢(𝑥, 𝑡) = 𝑈 (𝑥, 𝑡) + 𝑢⋆ (𝑥)

with 𝑈 (𝑥, 𝑡) being the dynamic deflection around the static equilibrium position. To impose the nonlinear boundary condition given by Eq. (11c), the Taylor series expansion around the equilibrium position is used:

(21)

Coefficients of sin(2𝛺𝑡): 𝜓2 (0) = 0, 𝜓2′ (0) = 0, 𝜓2′′ (1) = 0 ( 𝐺2 𝜑0 (1), 𝜑1 (1), 𝜑2 (1), 𝜓1 (1), 𝜓2 (1), 𝜑′′′ 0 (1), 𝜑′′′ 1 (1), 𝜑′′′ 2 (1), ) 𝜓 ′′′ 1 (1), 𝜓 ′′′ 2 (1), 𝑌 , 𝑍 = 0

− (22)

𝐻 (𝑍 − 𝑈 (1, 𝑡) − 𝑢⋆ (1))2

= 𝐻[−1 − 2𝑈 (1, 𝑡) − 3𝑈 2 (1, 𝑡) − 4𝑈 3 (1, 𝑡)]

(29)

Substitution of Eq. (28) into Eqs. (11a)–(11c) and Eq. (29) yields the equation governing the motion of the microcantilever about its equilibrium position and associated boundary conditions:

where functions 𝐹0 , 𝐹1 , 𝐹2 , 𝐺1 and 𝐺2 appearing in Eqs. (18)–(22) come from expansion of the nonlinear boundary condition. These functions are presented in the Appendix. Eqs. (13) through (17) are linear ordinary differential equations with constant coefficients. Such equations are solved easily. Solution to Eq. (13) is given by: 𝜑0 (𝑥) = 𝐴𝑥3 + 𝐵𝑥2 + 𝐶𝑥 + 𝐷

(28)

𝑈 𝑖𝑣 (𝑥, 𝑡) + 𝑈̈ (𝑥, 𝑡) + 𝑐 𝑈̇ (𝑥, 𝑡) = 0 𝑌 𝑈 (0, 𝑡) = (𝑒𝑖𝛺𝑡 + 𝑐𝑐), U′ (0, t) = 0, U′′ (l, t) = 0 2 ′′′ U (1, t) = 𝛽1 U(1, 𝑡) + 𝛽2 U2 (1, 𝑡) + 𝛽3 U3 (1, 𝑡)

(30a) (30b) (30c)

where:

(23)

𝛽1 = −2𝐻, 𝛽2 = −3𝐻, 𝛽3 = −4𝐻

Four unknown coefficients appear in Eq. (23): 𝐴, 𝐵, 𝐶 and 𝐷. In addition, general solutions to Eqs. (14)–(17) may be assumed to have the following forms:

Substituting these expressions in Eqs. (14)–(17) yields:

It should be noted that the term 𝑐𝑐 in Eq. (30b) represents complex conjugate of the preceding term. Next, we use the method of multiple scales to solve the governing equations given by Eqs. (30a)–(30c). However, no discretization scheme is used here and the partial differential equation is solved directly by assuming a solution in the following form:

(𝜆41 − 𝛺2 )𝐸 + 𝑐𝐹 𝛺 = 0

𝑈 (𝑥, 𝑇0 , 𝑇1 , 𝑇2 ) = 𝜀U1 (𝑥, 𝑇0 , 𝑇1 , 𝑇2 ) + 𝜀2 U2 (𝑥, 𝑇0 , 𝑇1 , 𝑇2 )

𝜑1 = 𝐸𝑒𝜆1 𝑥 , 𝜓1 = 𝐹 𝑒𝜆1 𝑥 , 𝜑2 = 𝐺𝑒𝜆2 𝑥 , 𝜓2 = 𝐻𝑒𝜆2 𝑥

(𝜆41 − 𝛺2 )𝐹 − 𝑐𝐸𝛺 = 0 (𝜆42 (𝜆42

(24)

+𝜀3 U3 (𝑥, 𝑇0 , 𝑇1 , 𝑇2 ) + ⋯

(25)

2

− 4𝛺 )𝐺 + 2𝑐𝐻𝛺 = 0

where 𝑇0 , 𝑇1 and 𝑇2 are different time scales with 𝑇𝑛 = Since a true approximation requires the damping, forcing, and the cubic nonlinearity to occur at the same level, we have rescaled damping and force coefficients [46]. The method of multiple scales splits a problem into several samescaled problems. By substituting Eq. (31) and rescaled coefficients into Eqs. (30a)–(30c) and separating terms of different orders, we obtain: Order of 𝜖:

− 4𝛺2 )𝐻 − 2𝑐𝐺𝛺 = 0

Linear algebraic equations given by Eq. (25) have nontrivial solutions if: (𝜆41 − 𝛺2 )2 + 𝑐 2 𝛺2 = 0

(26)

(𝜆42 − 4𝛺2 )2 + 4𝑐 2 𝛺2 = 0

Solving equations presented in Eq. (26) results in 8 values for 𝜆1 and also 8 values for 𝜆2 . Therefore, Eq. (24) is rewritten as: 𝜑1 = 𝜓2 =

8 ∑ 𝑖=1 8 ∑

𝐸𝑖 𝑒𝜆1𝑖 𝑥 , 𝜓1 =

(31) 𝜀𝑛 𝑡.

8 ∑ 𝑖=1

𝐹𝑖 𝑒𝜆1𝑖 𝑥 , 𝜑2 =

8 ∑ 𝑖=1

𝐺𝑖 𝑒𝜆2𝑖 𝑥 , (27)

𝑈1𝑖𝑣 (𝑥, 𝑡) + 𝑈̈ 1 (𝑥, 𝑡) = 0

(32a)

𝑈1 (0, 𝑡) = 0, 𝑈1′ (0, t) = 0, 𝑈1′′ (l, t) = 0

(32b)

𝑈1′′′ (l, t) =𝛽1 U1 (1, 𝑡)

(32c)

Order of

𝐻𝑖 𝑒𝜆2𝑖 𝑥

𝜖2 :

𝑈2𝑖𝑣 (𝑥, 𝑡) + 𝑈̈ 2 (𝑥, 𝑡)

𝑖=1

36

= −2𝐷0 𝐷1 𝑈1 (𝑥, 𝑡)

(33a)

M.S. Mahmoudi, A. Ebrahimian and A. Bahrami

International Journal of Non-Linear Mechanics 110 (2019) 33–43

𝑈2 (0, 𝑡) = 0, 𝑈2′ (0, t) = 0, 𝑈2′′ (l, t) = 0

(33b)

𝑔 ′′ (1, 𝑇2 ) = 0 ,

𝑈2′′′ (l, t) =𝛽1 U2 (1, 𝑡) + 𝛽2 U21 (1, 𝑡)

(33c)

𝑔 ′′′ (1, 𝑇2 ) = 𝛽1 𝑔(1, 𝑇2 ) + (

Order of 𝜖 3 :

(𝜔 −𝛺)

associated with the 𝑗th mode, defined as 𝜎𝑗 = 𝑗𝜖2 . We note that Eqs. (34a)–(34c) are not solved here and they are only used to determine a relation for 𝐴𝑗 , introduced in Eq. (35). To this end, we write 𝐴𝑗 in the polar form 𝐴𝑗 = 12 𝑝𝑗 𝑒𝑖𝑞𝑗 𝑇2 (𝑝𝑗 and 𝑞𝑗 are unknown functions of 𝑇2 ), substitute it into Eq. (35) and impose the solvability condition to obtain the following equation: √ [ ]2 ⎡ √ √[ ]2 1 ⎢ √ 1 ′′′ (0)𝑌 2 (𝑥)𝑑𝑥 √ ± 𝜎𝑗 = 𝜓 − 𝜔 𝑐𝑝 𝜓 ⎢ 𝑗 𝑗∫ 𝑗 𝑗 1 0 2𝑝𝑗 𝜔𝑗 ∫0 𝜓𝑗2 (𝑥)𝑑𝑥 ⎢ ⎣ ( )⎤ √ −4𝛽22 ⎥ 1 + 𝑝𝑗 3 + 2𝛽22 𝐺( 2𝜔𝑗 , 1) + 3𝛽3 ⎥ (43) 4 3 + 𝛽1 ⎥ ⎦

(34a) 𝑌 𝑖𝛺𝑡 (𝑒 + 𝑐𝑐), 𝑈3′ (0, t) = 0, 𝑈3′′ (l, t) = 0 2 ′′′ 𝑈3 (1, t) =𝛽1 U3 (1, 𝑡) + 2𝛽2 U1 (1, 𝑡)𝑈2 (1, 𝑡) + 𝛽3 U33 (1, 𝑡)

(34b) (34c)

Now, we solve equations associated with different orders of 𝜖: Order of 𝜖 problem: solution to Eqs. (32a)–(32c) for the 𝑗th mode is given by: 𝑈1 (𝑥, 𝑇0 , 𝑇1 , 𝑇2 ) = 𝐴𝑗 𝑒𝑖𝜔𝑗 𝑇0 𝜓𝑗 (𝑥) + 𝑐𝑐

(35)

where 𝐴𝑗 are unknown coefficients which will be found later and 𝜔𝑗 is the 𝑗th nondimensional natural frequency of the system. The nondimensional excitation frequency 𝛺 is assumed here to be near one of the natural frequencies 𝜔𝑗 . In addition, 𝜓𝑗 (𝑥), mode shapes of the linearized system, are given by: 𝜓𝑗 (𝑥) = 𝜓0 (sin 𝛾𝑗 𝑥 − sinh 𝛾𝑗 𝑥 − 𝐻0 (cos 𝛾𝑗 𝑥 − cosh 𝛾𝑗 𝑥))

Furthermore, the phase shift corresponding to the 𝑗th mode 𝜙𝑗 is given by:

(36) −1

𝜙𝑗 = tan

where: 𝜓0 =

cos 𝛾𝑗 + cosh 𝛾𝑗 2(sin 𝛾𝑗 cosh 𝛾𝑗 − sinh 𝛾𝑗 cos 𝛾𝑗 )

, 𝐻0 =

sin 𝛾𝑗 + sinh 𝛾𝑗 cos 𝛾𝑗 + cosh 𝛾𝑗

⎞ ⎛ 1 ⎟ ⎜ −𝜔𝑗 𝑐𝑝𝑗 ∫0 𝜓𝑗2 (𝑥)𝑑𝑥 ⎟ (44) ⎜ ( ) 2 ⎟ ⎜ √ 1 2 1 3 −4𝛽2 2 ∫ 2 𝜎 𝑝 𝜔 𝜓 (𝑥)𝑑𝑥 − 𝑝 + 2𝛽 𝐺( 2𝜔 , 1) + 3𝛽 ⎜ 𝑗 𝑗 𝑗 0 𝑗 3 ⎟ 𝑗 2 8 𝑗 3+𝛽1 ⎠ ⎝

Finally, the total deflection of the AFM tip may be represented as:

Value of dimensionless wave number 𝛾𝑗 is determined by solving this nonlinear characteristic equation [43]: 𝛾 3 (cosh 𝛾 cos 𝛾 + 1) + 𝛽1 (sin 𝛾 cosh 𝛾 − sinh 𝛾 cos 𝛾) = 0

(3 + 𝛽1 )

√ (42c) + 2𝛽2 2 𝐺( 2𝜔𝑗 , 1) + 3𝛽3 )𝐴𝑗 2 𝐴̄𝑗

where a dot indicates derivative with respect to 𝑇2 and a prime denotes derivative with respect to 𝑥. Moreover, 𝜎𝑗 is a detuning parameter

( ) 𝑈3𝑖𝑣 (𝑥, 𝑡) + 𝑈̈ 3 (𝑥, 𝑡) = −2𝐷0 𝐷1 𝑈2 (𝑥, 𝑡) − 2𝐷0 𝐷2 + 𝐷1 2 + 𝐶𝐷0 𝑈1 (𝑥, 𝑡) 𝑈3 (0, 𝑡) =

−4𝛽22

𝜀2 𝛽2 𝑝𝑗 2 ( ) 𝑢(1, 𝑡) = 𝑈 (1, 𝑡) + 𝑢⋆ (1) = − ( ) + 𝑢⋆ + 𝜀𝑝𝑗 cos 𝛺𝑡 − 𝜙𝑗 2 3 + 𝛽1 (45) √ 𝜀2 𝛽2 𝑝𝑗 2 𝐺( 2𝜔𝑗 , 1) ( ) + cos 2𝛺𝑡 − 2𝜙𝑗 2 In addition, starting with Eq. (43) and differentiating with respect to 𝜎𝑗 , maximum value of the tip amplitude corresponding to the 𝑗th mode is determined:

(37)

Moreover, the natural frequencies are related to wave numbers through equation 𝜔𝑗 4 = 𝛾𝑗 2 . Order of 𝜖 2 problem: 𝑈1 given by Eq. (35) is substituted into Eqs. (33a)–(33c). Imposing the solvability condition for the non-homogeneous part of Eqs. (33a)–(33c) implies that 𝐷1 𝐴𝑗 = 0, which means 𝐴𝑗 is a function of 𝑇2 only. Solution to order of 𝜖 2 problem may be presented as: √ 𝛽2 𝐴𝑗 𝐴̄ 𝑗 2 (38) 𝑈2 = ( ) 𝑥 (𝑥 − 3) + 𝛽2 𝐴𝑗 2 𝐺( 2𝜔𝑗 , 𝑥)𝑒2𝑖𝜔𝑗 𝑇0 + 𝑐𝑐 2 3 + 𝛽1

𝑝max𝑗 =

2𝜓𝑗′′′ (0)𝑌 1

𝜔𝑗 𝑐 ∫0 𝜓𝑗2 (𝑥)𝑑𝑥

occurring at: ) ( )2 ( −4𝛽 2 √ 2 𝐺( 2𝜔 , 1) + 3𝛽 2 𝜓𝑗′′′ (0)𝑌 + 2𝛽 𝑗 3 2 3+𝛽1 𝜎𝑗 = ( )3 1 2𝜔𝑗 3 ∫0 𝜓𝑗2 (𝑥)𝑑𝑥 𝑐 2

where the function 𝐺(𝜆, 𝑥) is given by: (cos 𝜆 + cosh 𝜆)(sinh 𝜆𝑥 − sin 𝜆𝑥) + (cos 𝜆 + cosh 𝜆)(sinh 𝜆𝑥 − sin 𝜆𝑥) 𝐺(𝜆, 𝑥) = 2(𝜆3 (1 + cos 𝜆 cosh 𝜆) − 𝛽1 (sinh 𝜆 cos 𝜆 − sin 𝜆 cosh 𝜆))

(46)

(47)

For more details on the method of multiple scales, see [45] and [46].

(39) 4. Numerical results

Order of 𝜖 3 problem: substituting 𝑈1 and 𝑈2 given by Eqs. (35) and (38) into Eqs. (34a)–(34c) yields:

In this section, numerical results are presented for the AFM with geometrical and material properties listed in Table 1. First, the accuracy of our proposed method (i.e., direct harmonic balance) is verified by comparing the results with those obtained based on the perturbation method. Then, the direct harmonic balance method is used to study nonlinear oscillations of the AFM probe. It should be noted that all figures presented in this section are plotted in terms of nondimensional quantities. For all numerical calculations, performed in this section the initial tip/sample separation distance is assumed to be 𝑍̂ = 8 nm, unless otherwise mentioned.

𝑈3 𝑖𝑣 + 𝐷02 𝑈3 = −𝑖𝜔𝑗 (2𝐴̇ 𝑗 + 𝑐𝐴𝑗 )𝑒𝑖𝜔𝑗 𝑇0 𝜓𝑗 (𝑥) (40a) 𝑌 𝑖𝛺𝑇0 ′ ′′ 𝑈3 (0, 𝑡) = 𝑒 + 𝑐𝑐, U 3 (0, 𝑡)= 0, 𝑈3 (1, 𝑡) = 0 (40b) 2 √ −4𝛽22 2 2 ̄ 𝑖𝜔𝑗 𝑇0 ′′′ 𝑈3 (1, 𝑡) = 𝛽1 𝑈3 (1, 𝑡) + ( + 2𝛽2 𝐺( 2𝜔𝑗 , 1) + 3𝛽3 )𝐴𝑗 𝐴𝑗 𝑒 3 + 𝛽1 √ +(𝛽3 + 2𝛽2 𝐺( 2𝜔𝑗 , 1))𝐴𝑗 3 𝑒3𝑖𝜔𝑗 𝑇0 + 𝑐𝑐 (40c) Now, we assume that the solution to Eqs. (40a) to (40c) is of the form: 𝑈3 = 𝑔(𝑥, 𝑇2 )𝑒𝑖𝜔𝑗 𝑇0 + 𝑐𝑐

(41)

4.1. Method verification

Substituting Eq. (41) into Eqs. (40a)–(40c) results in the governing equation for 𝑔 which may be represented as: 𝑔 𝑖𝑣 − 𝜔𝑗 2 𝑔 = −𝑖𝜔𝑗 (2𝐴̇ 𝑗 + 𝑐𝐴𝑗 )𝜓𝑗 (𝑥) 𝑌 𝑔(0, 𝑇2 ) = 𝑒𝑖𝜎𝑗 𝑇2 , 𝑔 ′ (0, 𝑇2 ) = 0 2

To verify the direct harmonic balance method, introduced in previous sections, the method of multiple scales is employed. Nonlinear governing equations of the system are solved for different values of excitation amplitude using both method of multiple scales and direct

(42a) (42b) 37

M.S. Mahmoudi, A. Ebrahimian and A. Bahrami

International Journal of Non-Linear Mechanics 110 (2019) 33–43

method and the method of multiple scales are in a very close agreement for excitations near first three natural frequencies. The accuracy and effectiveness of the direct harmonic balance method have been demonstrated for different excitation frequencies, so far. In the next section, the direct harmonic balance method is employed to study nonlinear oscillations of the AFM microcantilever with focus on higher modes and higher harmonics. It should be mentioned here that both methods presented in this paper predict steady-state responses, only. 4.2. Nonlinear oscillations of the AFM probe In this section nonlinear oscillations of AFM is investigated through several numerical examples based on the direct harmonic balance method. Fig. 5 displays the frequency-response curve for the system near its first three natural frequencies. Because of numerical limits, this figure does not include antiresonance points [47], but the direct harmonic balance method determines resonance frequencies with a high accuracy. As observed in the figure, excitation of the higher mode shapes results in a larger tip amplitude for the same value of the excitation amplitude. In the other words, sensitivity to excitation amplitude increases for higher mode shapes. Detailed frequency-response curves for different values of excitation amplitude are depicted in Fig. 6 for 𝑌̂ = 0.001 nm, 𝑌̂ = 0.002 nm and 𝑌̂ = 0.003 nm. As seen in the figure, a similar trend is observed for excitations near first three natural frequencies. The higher the excitation amplitude, the larger the tip amplitude is. In addition, numerical studies conducted by the authors reveal that amplitude of the response varies almost linearly with amplitude of the excitation. To study effects of the initial tip/sample separation distance on the dynamics of AFM probe, frequency-response curves of the system are presented in Fig. 7 for 𝑌̂ = 0.003 nm and several separation distances. Depending on how close the tip is to the sample in the reference configuration, frequency-response curves change. It is seen all the frequency-response curves follow a similar trend. However, a frequency shift occurs when the initial tip/sample separation distance varies. Moreover, it is observed that higher mode shapes exhibit a smaller frequency shift. This frequency shift is not identical for all values of the initial tip/sample separation distance. For small values of 𝑍̂ change in 𝑍̂ results in a larger frequency shift compared to change in

Fig. 3. Time history for tip displacement for 𝑌̂ = 0.003 nm and 𝛺̂ = 48.1 kHz.

harmonic balance method. Fig. 3 illustrates a comparison between time histories for the tip displacement obtained based on direct harmonic balance and perturbation methods. As seen in the figure, results are in excellent agreement. It should be noted that the excitation frequency is close to the first natural frequency of the system. First three natural frequencies of the system are 48.9 kHz, 306.6 kHz and 858.2 kHz, respectively. It should be also noted that to implement the direct harmonic balance method, the device is not necessarily required to be excited near a specific frequency (e.g., natural frequencies). In fact, the method may be implemented with the same numerical and analytical procedures even if the device is excited far from its natural frequencies. However, nonlinear oscillations of the AFM microcantilever are studied here near its first three natural frequencies using both the direct harmonic balance and the method of multiple scales. Phase portraits of the system are plotted in Fig. 4 for three different values of the excitation amplitude; that are 𝑌̂ = 0.001 nm, 𝑌̂ = 0.002 nm and 𝑌̂ = 0.003 nm. As observed in the figure, results obtained based on the direct harmonic balance

Fig. 4. Phase portraits of the system for excitation of the first three modes with 𝑌̂ = 0.001 nm, 𝑌̂ = 0.002 nm, 𝑌̂ = 0.003 nm and 𝑍̂ = 8 nm.

38

M.S. Mahmoudi, A. Ebrahimian and A. Bahrami

International Journal of Non-Linear Mechanics 110 (2019) 33–43

Fig. 5. Frequency response for a wide range of frequency for 𝑌̂ = 0.003 nm and 𝑍̂ = 8 nm.

Fig. 6. Frequency-response curves near first three natural frequencies for 𝑌̂ = 0.001 nm, 𝑌̂ = 0.002 nm, 𝑌̂ = 0.003 nm and 𝑍̂ = 8 nm.

̂ In addition, by reducing frequency shift, occurring for larger values of 𝑍. 𝑍̂ the softening behavior of the system becomes more apparent as the frequency-response curve is bent slightly to the left.

the second harmonic are shown in Fig. 9 for several values of the ̂ the second initial tip/sample separation gap. By changing value of 𝑍, harmonic experiences both frequency shift and a drastic change in the amplitude peak. Numerical results imply that the second harmonic exhibits a higher sensitivity to variation of initial gap rather than the first harmonic. We note that, in contrast to the first harmonic, the amplitude of the second harmonic decreases with initial tip/sample gap ̂ This may be attributed to the fact that when the tip comes close 𝑍. to the sample, tip/sample nonlinear interaction forces increase which consequently results in a higher amplitude for the second harmonic. As seen in the figure, for a rather large initial separation distance (i.e., 𝑍̂ = 15 nm) the second harmonic part of the response is almost zero compared to smaller initial separation gaps. We also note that amplitude of the second harmonic is, in fact, related to square of the excitation amplitude. This behavior of the second harmonic is also different from that of the first harmonic. For instance, if the excitation amplitude is doubled, the second harmonic amplitude is multiplied by 4. For the first harmonic, however, the response amplitude is doubled when the

So far, we have investigated the nonlinear oscillations of the AFM tip subjected to harmonic base excitations near its first three natural frequencies and presented numerical results for each excitation. Now, to obtain more information about the sample, we explore the second harmonic of the response for the same excitations. Fig. 8 presents the frequency-response curves for second harmonic of the response when the system is excited near its first natural frequency. As seen in the figure, magnitude of the second harmonic is very small in comparison with that of the first harmonic. In fact, the second harmonic is, generally, at least 1000 times smaller than the first harmonic [48] and is only a small part of the response. Nevertheless, due to its desirable properties, this small part may be used as a signal to improve the information, obtained by atomic force microscopy. In addition, amplitude of the second harmonic is more sensitive to amplitude of the excitation compared to that of the first harmonic. Frequency-response curves for 39

M.S. Mahmoudi, A. Ebrahimian and A. Bahrami

International Journal of Non-Linear Mechanics 110 (2019) 33–43

Fig. 7. Frequency-response curves near first three natural frequencies for 𝑌̂ = 0.003 nm and different values of the initial tip/sample separation distance.

Fig. 8. Second harmonic of the system for excitations near its first three natural frequencies for 𝑌̂ = 0.001 nm, 𝑌̂ = 0.002 nm, 𝑌̂ = 0.003 nm and 𝑍̂ = 8 nm.

excitation amplitude is doubled. Furthermore, we have studied the effect of viscous damping on the response amplitude. Our numerical studies indicate that decreasing the damping coefficient has exactly the same effect as increasing the excitation amplitude.

approaches to solve such problems include reducing the PDE to a set of ordinary differential equations and then applying numerical or semi-analytical methods to solve ODEs. The process of discretization is complex itself and may introduce numerical errors and reduce accuracy of final results. In this paper, however, a simple and efficient method which can be applied to a wide range of nonlinear problems has been introduced. This method, referred to as the direct harmonic balance method, does not require any discretization of governing equations. We have successfully applied the direct harmonic balance method to solve the initial–boundary problem governing the motion of the AFM microcantilever. A major advantage of the direct harmonic balance

5. Conclusions In a vast number of nonlinear oscillation problems, response to a harmonic excitation involves higher harmonics in addition to the fundamental frequency. On the other hand, real engineering problems are usually described by partial differential equations. Conventional 40

M.S. Mahmoudi, A. Ebrahimian and A. Bahrami

International Journal of Non-Linear Mechanics 110 (2019) 33–43

Fig. 9. Second harmonic of the system for excitations near its first three natural frequencies for 𝑌̂ = 0.003 and different values of the initial tip/sample separation distance.

method is that its implementation is very simple. On the other hand, contributions of all modes and all harmonics can be detected with no limit or complexity. After verifying our proposed method with the method of multiple scales, it has been used to study dynamics of the AFM with a focus on higher harmonics and higher mode shapes. It is found out that at the at the same conditions, excitation of a higher mode leads to a larger oscillation amplitude but less sensitivity to the initial tip/sample separation distance and a smaller frequency shift. The second harmonic response has been also studied. The second harmonic of the response in the non-contact mode is much smaller than the first harmonic and this makes it hard to detect it in the air ambient. However, our numerical simulations reveal that it is very sensitive to change in the initial tip/sample separation distance. In addition, while the main harmonic shows only a frequency shift by decreasing the initial tip/sample separation distance, second harmonic experiences both a frequency shift and a considerable change in its amplitude. Magnitude of the amplitude for the second harmonic grows quickly by reducing the tip/sample separation distance. We note that existence of the higher harmonics in the response stems from the tip/sample nonlinear interactions. Since this nonlinear force grows as the initial tip/sample separation distance decreases, larger amplitudes are expected for the second harmonic at smaller tip/sample distances. The keen sensitivity of higher harmonics makes them desirable signals to explore topography of various samples.

− (𝜑′′′ (𝑥) 𝜓1 (𝑥)2 )∕4 − 𝑌 𝑍 𝜑′′′ (𝑥) + 𝑌 𝜑0 (𝑥) 𝜑′′′ (𝑥) + 𝑌 𝜑′′′ (𝑥) 𝜑1 (𝑥) 1 1 1 0 + (𝑌 𝜑1 (𝑥) 𝜑′′′ (𝑥))∕2 + (𝑌 𝜑′′′ (𝑥) 𝜑2 (𝑥))∕2 − 2 𝑍 𝜑0 (𝑥) 𝜑′′′ (𝑥) 1 1 0 − 𝑍 𝜑1 (𝑥) 𝜑′′′ (𝑥) − 𝑍 𝜑2 (𝑥) 𝜑′′′ (𝑥) 1 1 + 𝜑0 (𝑥) 𝜑1 (𝑥) 𝜑′′′ (𝑥) + 𝜑0 (𝑥) 𝜑2 (𝑥) 𝜑′′′ (𝑥) 1 1 + (𝜑1 (𝑥) 𝜑′′′ (𝑥) 𝜑2 (𝑥))∕2 + (𝑌 𝜓1 (𝑥) 𝜓2′′′ (𝑥))∕2 1 + (𝑌 𝜓1′′′ (𝑥) 𝜓2 (𝑥))∕2 − 𝑍 𝜓1 (𝑥) 𝜓1′′′ (𝑥) − 𝑍 𝜓2 (𝑥) 𝜓2′′′ (𝑥) + 𝜑0 (𝑥) 𝜓1 (𝑥) 𝜓1′′′ (𝑥) + 𝜑0 (𝑥) 𝜓2 (𝑥) 𝜓2′′′ (𝑥) + (𝜑1 (𝑥) 𝜓1 (𝑥) 𝜓2′′′ (𝑥))∕2 + (𝜑1 (𝑥) 𝜓1′′′ (𝑥) 𝜓2 (𝑥))∕2 + (𝜑′′′ (𝑥) 𝜓1 (𝑥) 𝜓2 (𝑥))∕2 1 − (𝜑2 (𝑥) 𝜓1 (𝑥) 𝜓1′′′ (𝑥))∕2 + 𝐻 𝐹1 (𝜑0 (𝑥), 𝜑1 (𝑥), 𝜑2 (𝑥), 𝜓1 (𝑥), 𝜓1 (2), 𝜑0 ′′′ (𝑥), 𝜑′′′ (𝑥), 𝜑2 ′′′ (𝑥), 𝜓1′′′ (𝑥), 𝜓2 ′′′ (2), 𝑌 , 𝑍) = 1 (3 𝑌 2 𝜑′′′ (𝑥))∕4 + 𝑍 2 𝜑′′′ (𝑥) + 𝜑0 (𝑥)2 𝜑′′′ (𝑥) 1 1 1 + (3 𝜑1 (𝑥)2 𝜑′′′ (𝑥))∕4 + (𝜑′′′ (𝑥) 𝜑2 (𝑥)2 )∕2 1 1 + (𝜑′′′ (𝑥) 𝜓1 (𝑥)2 )∕4 + (𝜑′′′ (𝑥) 𝜓2 (𝑥)2 )∕2 − 2 𝑌 𝑍 𝜑′′′ (𝑥) 1 1 0 − 𝑌 𝑍 𝜑′′′ (𝑥) + 2 𝑌 𝜑0 (𝑥) 𝜑′′′ (𝑥) 1 0 + 𝑌 𝜑0 (𝑥) 𝜑′′′ (𝑥) + 𝑌 𝜑′′′ (𝑥) 𝜑2 (𝑥) 1 0 + (3 𝑌 𝜑1 (𝑥) 𝜑′′′ (𝑥))∕2 + 𝑌 𝜑2 (𝑥) 𝜑′′′ (𝑥)− 1 1

Appendix

2 𝑍 𝜑0 (𝑥) 𝜑′′′ (𝑥) − 2 𝑍 𝜑′′′ (𝑥) 𝜑1 (𝑥) 1 0 − 𝑍 𝜑1 (𝑥) 𝜑′′′ (𝑥) − 𝑍 𝜑′′′ (𝑥) 𝜑2 (𝑥) 1 1

𝐹0 (𝜑0 (𝑥), 𝜑1 (𝑥), 𝜑2 (𝑥), 𝜓1 (𝑥), 𝜓1 (2), 𝜑0 ′′′ (𝑥), 𝜑′′′ (𝑥), 𝜑2 ′′′ (𝑥), 1

+ 2 𝜑0 (𝑥) 𝜑′′′ (𝑥) 𝜑1 (𝑥) + 𝜑0 (𝑥) 𝜑1 (𝑥) 𝜑′′′ (𝑥) 0 1

𝜓1′′′ (𝑥), 𝜓2 ′′′ (2), 𝑌 , 𝐻, 𝑍) = (𝑌 2 𝜑′′′ (𝑥))∕2 + (𝑌 2 𝜑′′′ (𝑥))∕4 + 𝑍 2 𝜑′′′ (𝑥) 0 1 0 2 ′′′ ′′′ 2 + 𝜑0 (𝑥) 𝜑0 (𝑥) + (𝜑0 (𝑥) 𝜑1 (𝑥) )∕2 + (𝜑′′′ (𝑥) 𝜑2 (𝑥)2 )∕2 + (𝜑1 (𝑥)2 𝜑′′′ (𝑥))∕4 0 1 ′′′ 2 + (𝜑0 (𝑥) 𝜓1 (𝑥)2 )∕2 + (𝜑′′′ (𝑥) 𝜓 (𝑥) )∕2 2 0

+ 𝜑0 (𝑥) 𝜑′′′ (𝑥) 𝜑2 (𝑥) + 𝜑′′′ (𝑥) 𝜑1 (𝑥) 𝜑2 (𝑥) 1 0 + 𝜑1 (𝑥) 𝜑2 (𝑥) 𝜑′′′ (𝑥) + (𝑌 𝜓1 (𝑥) 𝜓1′′′ (𝑥))∕2 1 + 𝑌 𝜓2 (𝑥) 𝜓2′′′ (𝑥) − 𝑍 𝜓1 (𝑥) 𝜓2′′′ (𝑥) − 𝑍 𝜓1′′′ (𝑥) 𝜓2 (𝑥) + 𝜑0 (𝑥) 𝜓1 (𝑥) 𝜓2′′′ (𝑥) + 𝜑0 (𝑥) 𝜓1′′′ (𝑥) 𝜓2 (𝑥) + 𝜑′′′ (𝑥) 𝜓1 (𝑥) 𝜓2 (𝑥) 0 + (𝜑1 (𝑥) 𝜓1 (𝑥) 𝜓1′′′ (𝑥))∕2 + 𝜑1 (𝑥) 𝜓2 (𝑥) 𝜓2′′′ (𝑥) 41

M.S. Mahmoudi, A. Ebrahimian and A. Bahrami

International Journal of Non-Linear Mechanics 110 (2019) 33–43 [6] R. Magerle, Nanotomography, Phys. Rev. Lett. 85 (13) (2000) 2749. [7] J.C. Fernandes, V. Mareau, L. Gonon, Afm-raman colocalization setup: Advanced characterization technique for polymers, Int. J. Polymer Analy. Charact. 23 (2) (2018) 113–119. [8] F.J. Giessibl, Atomic resolution of the silicon (111)-(7x7) surface by atomic force microscopy, Science 267 (5194) (1995) 68–71. [9] H. Labidi, M. Koleini, T. Huff, M. Salomons, M. Cloutier, J. Pitters, R.A. Wolkow, Indications of chemical bond contrast in afm images of a hydrogen-terminated silicon surface, Nature Commun. 8 (2017) 14222. [10] R. Garcıa, R. Perez, Dynamic atomic force microscopy methods, Surface Sci. Rep. 47 (6–8) (2002) 197–301. [11] G. Rega, U. Andreaus, L. Placidi, V. Settimi, Nonlinear dynamics of atomic force microscopy, in: International Conference on Nonlinear Dynamics in Engineering: Modeling, Analysis and Applications, Aberdeen, Scotland, August, 2013. [12] T.R. Rodrıguez, R. Garcıa, Compositional mapping of surfaces in atomic force microscopy by excitation of the second normal mode of the microcantilever, Appl. Phys. Lett. 84 (3) (2004) 449–451. [13] N. Martinez, S. Patil, J.R. Lozano, R. Garcia, Enhanced compositional sensitivity in atomic force microscopy by the excitation of the first two flexural modes, Appl. Phys. Lett. 89 (15) (2006) 153115. [14] J.R. Lozano, R. Garcia, Theory of phase spectroscopy in bimodal atomic force microscopy, Phys. Rev. B 79 (1) (2009) 014110. [15] J.R. Lozano, D. Kiracofe, J. Melcher, R. Garcia, A. Raman, Calibration of higher eigenmode spring constants of atomic force microscope cantilevers, Nanotechnology 21 (46) (2010) 465502. [16] D. Kiracofe, A. Raman, On eigenmodes, stiffness, and sensitivity of atomic force microscope cantilevers in air versus liquids, J. Appl. Phys. 107 (3) (2010) 033506. [17] S.D. Solares, G. Chawla, Frequency response of higher cantilever eigenmodes in bimodal and trimodal tapping mode atomic force microscopy, Meas. Sci. Technol. 21 (12) (2010) 125502. [18] N. Martinez, J.R. Lozano, E. Herruzo, F. Garcia, C. Richter, T. Sulzbach, R. Garcia, Bimodal atomic force microscopy imaging of isolated antibodies in air and liquids, Nanotechnology 19 (38) (2008) 384011. [19] D. Forchheimer, R. Forchheimer, D.B. Haviland, Improving image contrast and material discrimination with nonlinear response in bimodal atomic force microscopy, Nature Commun. 6 (2015) 6270. [20] R. Garcia, E.T. Herruzo, The emergence of multifrequency force microscopy, Nature Nanotechnol. 7 (4) (2012) 217. [21] O. Sahin, C.F. Quate, O. Solgaard, A. Atalar, Resonant harmonic response in tapping-mode atomic force microscopy, Phys. Rev. B 69 (16) (2004) 165416. [22] M. Stark, R.W. Stark, W.M. Heckl, R. Guckenberger, Inverting dynamic force microscopy: From signals to time-resolved interaction forces, Proc. Natl. Acad. Sci. 99 (13) (2002) 8473–8478. [23] R.W. Stark, W.M. Heckl, Higher harmonics imaging in tapping-mode atomic-force microscopy, Rev. Sci. Instrum. 74 (12) (2003) 5111–5114. [24] A. Belianinov, S.V. Kalinin, S. Jesse, Complete information acquisition in dynamic force microscopy, Nature Commun. 6 (2015) 6550. [25] U. Andreaus, L. Placidi, G. Rega, Microcantilever dynamics in tapping mode atomic force microscopy via higher eigenmodes analysis, J. Appl. Phys. 113 (22) (2013) 224302. [26] S. Kawai, S.-i. Kitamura, D. Kobayashi, S. Meguro, H. Kawakatsu, An ultrasmall amplitude operation of dynamic force microscopy with second flexural mode, Appl. Phys. Lett. 86 (19) (2005) 193107. [27] O. Pfeiffer, C. Loppacher, C. Wattinger, M. Bammerlin, U. Gysin, M. Guggisberg, S. Rast, R. Bennewitz, E. Meyer, H.-J. Güntherodt, Using higher flexural modes in non-contact force microscopy, Appl. Surf. Sci. 157 (4) (2000) 337–342. [28] R.D. Turner, J. Kirkham, D. Devine, N.H. Thomson, Second harmonic atomic force microscopy of living staphylococcus aureus bacteria, Appl. Phys. Lett. 94 (4) (2009) 043901. [29] A. Raman, S. Trigueros, A. Cartagena, A. Stevenson, M. Susilo, E. Nauman, S.A. Contera, Mapping nanomechanical properties of live cells using multi-harmonic atomic force microscopy, Nature Nanotechnol. 6 (12) (2011) 809. [30] O. Sahin, G. Yaralioglu, R. Grow, S. Zappe, A. Atalar, C. Quate, O. Solgaard, High-resolution imaging of elastic properties using harmonic cantilevers, Sensors Actuators A 114 (2–3) (2004) 183–190. [31] O. Sahin, S. Magonov, C. Su, C.F. Quate, O. Solgaard, An atomic force microscope tip designed to measure time-varying nanomechanical forces, Nature Nanotechnol. 2 (8) (2007) 507. [32] M. Dong, O. Sahin, A nanomechanical interface to rapid single-molecule interactions, Nature Commun. 2 (2011) 247. [33] W. van de Water, J. Molenaar, Dynamics of vibrating atomic force microscopy, Nanotechnology 11 (3) (2000) 192. [34] L. Nony, R. Boisgard, J.-P. Aimé, Nonlinear dynamical properties of an oscillating tip–cantilever system in the tapping mode, J. Chem. Phys. 111 (4) (1999) 1615– 1627. [35] X.-H. Long, G. Lin, B. Balachandran, Grazing bifurcations in an elastic structure excited by harmonic impactor motions, Physica D 237 (8) (2008) 1129–1138. [36] I. Chakraborty, B. Balachandran, Noise influenced elastic cantilever dynamics with nonlinear tip interaction forces, Nonlinear Dynam. 66 (3) (2011) 427–439.

𝐹2 (𝜑0 (𝑥), 𝜑1 (𝑥), 𝜑2 (𝑥), 𝜓1 (𝑥), 𝜓1 (2), 𝜑0 ′′′ (𝑥), 𝜑′′′ (𝑥), 𝜑2 ′′′ (𝑥), 𝜓1′′′ (𝑥), 𝜓2 ′′′ (2), 𝑌 , 𝑍) = 1 (𝑌 2 𝜑′′′ (𝑥))∕2 + (𝑌 2 𝜑′′′ (𝑥))∕2 + 𝑍 2 𝜑′′′ (𝑥) 0 1 1 ′′′ 2 2 ′′′ + (𝜑0 (𝑥) 𝜑1 (𝑥) )∕2 + 𝜑0 (𝑥) 𝜑1 (𝑥) + (𝜑1 (𝑥)2 𝜑′′′ (𝑥))∕2 + (3 𝜑2 (𝑥)2 𝜑′′′ (𝑥))∕4 1 1 ′′′ 2 ′′′ 2 − (𝜑0 (𝑥) 𝜓1 (𝑥) )∕2 + (𝜑1 (𝑥) 𝜓1 (𝑥) )∕2 + (𝜑′′′ (𝑥) 𝜓2 (𝑥)2 )∕4 − 𝑌 𝑍 𝜑′′′ (𝑥) 1 1 ′′′ ′′′ + 𝑌 𝜑0 (𝑥) 𝜑′′′ (𝑥) + 𝑌 𝜑 (𝑥) 𝜑 1 (𝑥) + 𝑌 𝜑1 (𝑥) 𝜑1 (𝑥) 1 0 ′′′ ′′′ + 𝑌 𝜑1 (𝑥) 𝜑2 (𝑥) − 2 𝑍 𝜑0 (𝑥) 𝜑1 (𝑥) − 2 𝑍 𝜑′′′ (𝑥) 𝜑2 (𝑥) − 𝑍 𝜑1 (𝑥) 𝜑′′′ (𝑥) 0 1 + 2 𝜑0 (𝑥) 𝜑′′′ (𝑥) 𝜑 (𝑥) + 𝜑 (𝑥) 𝜑 (𝑥) 𝜑′′′ (𝑥) + 𝜑1 (𝑥) 𝜑′′′ (𝑥) 𝜑2 (𝑥) 2 0 1 0 1 1 ′′′ ′′′ ′′′ + 𝑍 𝜓1 (𝑥) 𝜓1 (𝑥) − 𝜑0 (𝑥) 𝜓1 (𝑥) 𝜓1 (𝑥) + 𝜑2 (𝑥) 𝜓1 (𝑥) 𝜓1 (𝑥) + (𝜑2 (𝑥) 𝜓2 (𝑥) 𝜓2′′′ (𝑥))∕2 𝐺1 (𝜑0 (𝑥), 𝜑1 (𝑥), 𝜑2 (𝑥), 𝜓1 (𝑥), 𝜓1 (2), 𝜑0 ′′′ (𝑥), 𝜑′′′ (𝑥), 𝜑2 ′′′ (𝑥), 𝜓1′′′ (𝑥), 𝜓2 ′′′ (2), 𝑌 , 𝑍) = 1 2 ′′′ (𝑌 𝜓1 (𝑥))∕4 + 𝑍 2 𝜓1′′′ (𝑥) + 𝜑0 (𝑥)2 𝜓1′′′ (𝑥) + (𝜑1 (𝑥)2 𝜓1′′′ (𝑥))∕4 + (𝜑2 (𝑥)2 𝜓1′′′ (𝑥))∕2 + (3 𝜓1 (𝑥)2 𝜓1′′′ (𝑥))∕4 + (𝜓1′′′ (𝑥) 𝜓2 (𝑥)2 )∕2 − 𝑌 𝑍 𝜓2′′′ (𝑥) + 𝑌 𝜑0 (𝑥) 𝜓2′′′ (𝑥) + 𝑌 𝜑′′′ (𝑥) 𝜓2 (𝑥) + (𝑌 𝜑1 (𝑥) 𝜓1′′′ (𝑥))∕2 0 ′′′ + (𝑌 𝜑1 (𝑥) 𝜓1 (𝑥))∕2 − 2 𝑍 𝜑0 (𝑥) 𝜓1′′′ (𝑥) − 2 𝑍 𝜑′′′ (𝑥) 𝜓1 (𝑥) − 𝑍 𝜑1 (𝑥) 𝜓2′′′ (𝑥) 0 − 𝑍 𝜑′′′ (𝑥) 𝜓2 (𝑥) + 𝑍 𝜑2 (𝑥) 𝜓1′′′ (𝑥) 1 ′′′ + 𝑍 𝜑1 (𝑥) 𝜓1 (𝑥) + 2 𝜑0 (𝑥) 𝜑′′′ (𝑥) 𝜓1 (𝑥) + 𝜑0 (𝑥) 𝜑1 (𝑥) 𝜓2′′′ (𝑥) 0 ′′′ + 𝜑0 (𝑥) 𝜑1 (𝑥) 𝜓2 (𝑥) − 𝜑0 (𝑥) 𝜑2 (𝑥) 𝜓1′′′ (𝑥) − 𝜑0 (𝑥) 𝜑′′′ (𝑥) 𝜓1 (𝑥) 1 ′′′ + 𝜑′′′ (𝑥) 𝜑 (𝑥) 𝜓 (𝑥) − 𝜑 (𝑥) 𝜑 (𝑥) 𝜓1 (𝑥) 1 2 2 0 0 ′′′ + (𝜑1 (𝑥) 𝜑′′′ (𝑥) 𝜓 (𝑥))∕2 + 𝜑 (𝑥) 𝜑 (𝑥) 𝜓1 (𝑥) + 𝜓1 (𝑥) 𝜓2 (𝑥) 𝜓2′′′ (𝑥) 1 2 1 1 𝐺2 (𝜑0 (𝑥), 𝜑1 (𝑥), 𝜑2 (𝑥), 𝜓1 (𝑥), 𝜓1 (2), 𝜑0 ′′′ (𝑥), 𝜑′′′ (𝑥), 𝜑2 ′′′ (𝑥), 𝜓1′′′ (𝑥), 𝜓2 ′′′ (2), 𝑌 , 𝑍) = 1 ′′′ (𝜓2 (𝑥) 𝑌 2 )∕2 − 𝜓1′′′ (𝑥) 𝑌 𝑍 + 𝜓1′′′ (𝑥) 𝑌 𝜑0 (𝑥) + 𝜓2′′′ (𝑥) 𝑌 𝜑1 (𝑥) + 𝜑′′′ (𝑥) 𝑌 𝜓1 (𝑥) 0 ′′′ + 𝜑′′′ (𝑥) 𝑌 𝜓 (𝑥) + 𝜓 (𝑥) 𝑍 2 − 2 𝜓2′′′ (𝑥) 𝑍 𝜑0 (𝑥) 2 1 2 ′′′ ′′′ − 𝜓1 (𝑥) 𝑍 𝜑1 (𝑥) − 𝜑1 (𝑥) 𝑍 𝜓1 (𝑥) − 2 𝜑′′′ (𝑥) 𝑍 𝜓2 (𝑥) + 𝜓2′′′ (𝑥) 𝜑0 (𝑥)2 0 ′′′ + 𝜓1 (𝑥) 𝜑0 (𝑥) 𝜑1 (𝑥) + 𝜑′′′ (𝑥) 𝜑0 (𝑥) 𝜓1 (𝑥) 1 ′′′ + 2 𝜑′′′ (𝑥) 𝜑 (𝑥) 𝜓 (𝑥) + (𝜓 (𝑥) 𝜑1 (𝑥)2 )∕2 0 2 0 2 ′′′ ′′′ + 𝜑0 (𝑥) 𝜑1 (𝑥) 𝜓1 (𝑥) + 𝜑1 (𝑥) 𝜑1 (𝑥) 𝜓2 (𝑥) + (𝜓2′′′ (𝑥) 𝜑2 (𝑥)2 )∕4 + (𝜑′′′ (𝑥) 𝜑2 (𝑥) 𝜓2 (𝑥))∕2 + (𝜓2′′′ (𝑥) 𝜓1 (𝑥)2 )∕2 1 + 𝜓1′′′ (𝑥) 𝜓1 (𝑥) 𝜓2 (𝑥) + (3 𝜓2′′′ (𝑥) 𝜓2 (𝑥)2 )∕4 References [1] G. Binnig, C.F. Quate, C. Gerber, Atomic force microscope, Phys. Rev. Lett. 56 (9) (1986) 930. [2] G. Haugstad, Atomic Force Microscopy: Understanding Basic Modes and Advanced Applications, John Wiley & Sons, 2012. [3] C. Möller, M. Allen, V. Elings, A. Engel, D.J. Müller, Tapping-mode atomic force microscopy produces faithful high-resolution images of protein surfaces, Biophys. J. 77 (2) (1999) 1150–1158. [4] S. Dutta, C. Rivetti, N.R. Gassman, C.G. Young, B.T. Jones, K. Scarpinato, M. Guthold, Analysis.of. single, Analysis of single cisplatin-induced dna bends by atomic force microscopy and simulations, J. Mol. Recognit. (2018) e2731. [5] M. Li, L. Liu, X. Xu, X. Xing, D. Dang, N. Xi, Y. Wang, Nanoscale characterization of dynamic cellular viscoelasticity by atomic force microscopy with varying measurement parameters, J. Mech. Behav. Biomed. Mater. 82 (2018) 193–201. 42

M.S. Mahmoudi, A. Ebrahimian and A. Bahrami

International Journal of Non-Linear Mechanics 110 (2019) 33–43 [42] A. Bahrami, A.H. Nayfeh, Nonlinear dynamics of tapping mode atomic force microscopy in the bistable phase, Commun. Nonlinear Sci. Numer. Simul. 18 (3) (2013) 799–810. [43] J. Turner, Non-linear vibrations of a beam with cantilever-hertzian contact boundary conditions, J. Sound Vib. 275 (1–2) (2004) 177–191. [44] E.M. Abdel-Rahman, A.H. Nayfeh, Contact force identification using the subharmonic resonance of a contact-mode atomic force microscopy, Nanotechnology 16 (2) (2005) 199. [45] A.H. Nayfeh, D.T. Mook, Nonlinear Oscillations, John Wiley & Sons, 2008. [46] A.H. Nayfeh, Introduction to Perturbation Techniques, John Wiley & Sons, 2011. [47] R. García, Amplitude Modulation Atomic Force Microscopy, John Wiley & Sons, 2011. [48] T.R. Rodrıguez, R. Garcıa, Tip motion in amplitude modulation (tapping-mode) atomic-force microscopy: Comparison between continuous and point-mass models, Appl. Phys. Lett. 80 (9) (2002) 1646–1648.

[37] U. Andreaus, L. Placidi, G. Rega, Soft impact dynamics of a cantilever beam: equivalent sdof model versus infinite-dimensional system, Proc. Inst. Mech. Eng. C 225 (10) (2011) 2444–2456. [38] U. Andreaus, P. Baragatti, L. Placidi, Experimental and numerical investigations of the responses of a cantilever beam possibly contacting a deformable and dissipative obstacle under harmonic excitation, Int. J. Non-Linear Mech. 80 (2016) 96–106. [39] U. Andreaus, G. Rega, L. Placidi, Higher order eigenmodes in tapping mode atomic force microscopy, in: 4th Canadian Conference on Nonlinear Solid Mechanics, Montréal, Canada, July, 2013. [40] A.J. Dick, B. Balachandran, H. Yabuno, M. Numatsu, K. Hayashi, M. Kuroda, K. Ashida, Utilizing nonlinear phenomena to locate grazing in the constrained motion of a cantilever beam, Nonlinear Dynam. 57 (3) (2009) 335–349. [41] A. Bahrami, A.H. Nayfeh, On the dynamics of tapping mode atomic force microscope probes, Nonlinear Dynam. 70 (2) (2012) 1605–1617.

43