Analysis of the Frequency Shift versus Force Gradient of a Dynamic AFM Quartz Tuning Fork Subject to Lennard-Jones Potential Force

Abstract

A self-sensing and self-actuating quartz tuning fork (QTF) can be used to obtain its frequency shift as function of the tip-sample distance. Once the function of the frequency shift versus force gradient is acquired, the combination of these two functions results in the relationship between the force gradient and the tip-sample distance. Integrating the force gradient once and twice elucidates the values of the interaction force and the interatomic potential, respectively. However, getting the frequency shift as a function of the force gradient requires a physical model which can describe the equations of motion properly. Most papers have adopted the single harmonic oscillator model, but encountered the problem of determining the spring constant. Their methods of finding the spring constant are very controversial in the research community and full of discrepancies. By circumventing the determination of the spring constant, we propose a method which models the prongs and proof mass as elastic bodies. Through the use of Hamilton’s principle, we can obtain the equations of motion of the QTF, which is subject to Lennard-Jones potential force. Solving these equations of motion analytically, we get the relationship between the frequency shift and force gradient.

1. Introduction

Quartz tuning forks (QTF) are widely used as a scanning probe for atomic force microscopy (AFM). Due to the inherent precise oscillation frequency of single-crystal quartz, the frequency shift is measured as a signal for the atomic and subatomic resolution images of surface of metals, dielectrics, and semiconductors [1,2]. Pawlak et al. [3] used QTF-based non-contact AFM (nc-AFM) to observe the island growth process and the hexagonal structure of fullerene assembly on a KBr(001). They found that Cu(111) substrate has a stronger molecule-substrate interaction than KBr(001), so that Cu(111) substrate can hold the fullerene and not be disrupted by the nc-AFM when its tip approaches the short-range repulsive force regime. The first time that individual atoms within the aromatic molecules were resolved was done by means of nc-AFM [4]. Both techniques of nc-AFM and STM (scanning tunneling microscope) can be used complementarily to study atomic bonding within and between molecules, molecular orbitals, molecular conformational changes, and chemical reactions at the single-molecule level [4,5]. Polesel-Marris et al. [5] used nc-AFM to observe the unfolding signatures during protein stretching and they also functionalized the tip, so as to graft proteins onto the substrate. The force-distance curve of AFM can be used to detect the stiffness change of diseased cells compared with the surrounding cell wall [6].

Certain researchers have enhanced the functions of AFM by modifying the structure of a QTF probe [7,8,9,10,11]. Bayat [7] coupled a cantilever probe called the Akiyama probe to a QTF. The Akiyama probe oscillating at its second eigenmode can be used to measure soft samples in tapping mode, while the QTF oscillating at its first eigenmode is used to operate in either tapping mode or non-contact mode. Hussain [9] modified the structure of the tip and its orientation, such that the QTF can perform sidewall roughness measurements.

Although the output of nc-AFM can provide the frequency shift as a function of the distance between the tip and sample [12,13,14,15,16], the tip-sample interaction force cannot be obtained directly from this kind of experiment. First, we need the equations of motion, which can describe exactly the motion behavior of the system composed of a QTF and the tip-sample interaction force. Tip-sample interaction force is equal to the negative gradient of the tip-sample potential and the gradient of the force (or force gradient) is the spring constant of the physically-modeled spring connecting the tip and sample. By solving the equations of motion, the frequency shift as a function of the force gradient can be obtained. The combination of these two functions yields the force gradient as a function of the tip-sample distance. Then, the integration of the force gradient provides the tip-sample interaction force and the second integration provides the interaction potential.

References [13,14,15] have modeled a QTF as a simple one-dimensional harmonic oscillator and they use their own method for estimating the spring constant of this oscillator. Castellanos-Gomez [17] modeled the QTF as two collinear coupling oscillators and he evaluated the spring constant of a single prong, as well as the spring constant, of the coupling spring connecting the two prongs by theoretical and experimental analysis. Falter et al. [18] calibrated the beam formula for determining the static spring constant of a qPlus sensor by shifting the origin of the prong to the zero stress zone inside the basis and including the torsion effect. Melcher et al. [19] analyzed a qPlus QTF sensor, in which the tip is of finite length and has axial offset. They improved the accuracy of the spring constant estimation by correct derivation of the angular displacement of the finite tip during the deformation of the slender prong, and also by considering the effect of the peaked ridges appearing in the side wall of the prongs. In this paper, we model the two prongs and the proof mass as continuum bodies instead of the spring-mass oscillator model, and establish the equations of motion through the use of Hamilton’s principle, including the Lennard-Jones potential force between the tip and sample. Although our equations of motion are more complicated than those of harmonic oscillators, there is no need for us to determine the equivalent spring constant.

The free prong, the prong loaded with interaction force, and the proof mass are taken apart. With the interface surfaces acted on by shear forces and moments of opposite direction, their own equations of motion are set up independently, and their general solutions are solved analytically. Finally, putting all these three sets of general solutions together with the boundary and interface conditions imposed results in the characteristic equation, which can be solved to get the natural frequency of the in-plane anti-phase resonant mode of the QTF. Based on this, the frequency shift, which is the difference between the loaded QTF and free QTF, can be obtained.

2. The Equations of Motion and the General Solution of the Prong Loaded with Tip-Sample Interaction Force

2.1. The Linearization of the Lennard-Jones Potential Force

The potential function of Lennard-Jones type is considered here as [20,21,22]

$$W(z)=-\frac{2}{3}{\pi }^2ϵ{\rho }_1{\rho }_2{\sigma }^{5}R[\frac{\sigma }{s}-\frac{{\sigma }^7}{210s^7}]$$

(1)

where $ϵ$ is the dielectric constant, $\sigma$ is the radius of the sample molecule, $R$ is the radius of the atom of the tip, ${\rho }_1\mathrm{and}{\rho }_2$ are the molecule density of the tip and sample, respectively, and $s$ is the distance between the tip and sample. The interaction force is the gradient of the potential

$$F=-\frac{\partial W}{\partial s}=-\frac{A_{1}R}{180s^8}+\frac{A_{2}R}{6s^2}$$

(2)

Referring to Figure 1, the parameter $z$ is the initial height between the tip and sample, with the interaction force not yet considered. $w^{\ast }$ is the displacement of the tip in static equilibrium, $\overline{w}\left(x,t\right)$ is the displacement of vibration around the equilibrium state. So the actual distance between the tip and sample in the deformed state of the QTF is $\left(z-w\left(x,t\right)\right)$ where

$$w\left(x,t\right)=w^{\ast }\left(x\right)+\overline{w}\left(x,t\right)$$

(3)

is the total displacement of the tip. Therefore, Equation (2) becomes

$$F=-\frac{A_{1}R}{180{(z-w)}^8}+\frac{A_{2}R}{6{(z-w)}^2}$$

(4)

where $A_1={\pi }^2{\rho }_1{\rho }_2{c}_1$ and $A_2={\pi }^2{\rho }_1{\rho }_2{c}_2$ are called Hamaker constants. In the practical application of nc-AFM, the amplitude $\overline{w}$ is kept small, thus, $w\ll z$, so we can rewrite (4) as

$$F(w)=-\frac{1}{180}{A}_{1}R{z}^{-8}{(1-\frac{w}{z})}^{-8}+\frac{1}{6}{A}_{2}R{z}^{-2}{(1-\frac{w}{z})}^{-2}$$

(5)

Taking a binomial expansion of those terms ${\left(1-w/z\right)}^{-8}$ and ${\left(1-w/z\right)}^{-2}$, Equation (5) reduces to

$$F\left(w;z\right)=k\left(z\right)w\left(x,t;z\right)+F_0\left(z\right)$$

(6)

where

$$k=(-\frac{2}{45}{A}_{1}R{z}^{-9}+\frac{1}{3}{A}_{2}R{z}^{-3}), F_0=-\frac{1}{180}{A}_{1}R{z}^{-8}+\frac{1}{6}{A}_{2}R{z}^{-2}$$

(7)

Notice that given a value of $z$ in (7), $k$, and $F_0$ becomes constant, and Equation (6) reveals that $k$ plays the role of the spring constant of a linear spring and the Lennard-Jones force can be represented by the linear spring plus a constant force, as shown in Figure 2.

Taking a derivative of Equation (6) with respect to displacement $w$ gives

$$k=\frac{\partial F}{\partial w}$$

(8)

Since $s=z-w$, for a given value of the initial height $z$ $ds=-dw$, then Equation (8) becomes

$$k=\frac{\partial F}{\partial w}=\frac{\partial F}{\partial s}\frac{ds}{dw}=-\frac{\partial F}{\partial s}.$$

(9)

Thus the spring constant has the physical meaning of negative force gradient.

2.2. The Deformation and Equations of Motion

The QTF is made of a $Z-cut2°$ ingle-crystal wafer. The second-order stress ${\sigma }_{ij}$ and strain tensors ${\epsilon }_{ij}$ are expressed in the vector form as $\left({\sigma }_1,{\sigma }_2,{\sigma }_3,{\sigma }_4,{\sigma }_5,{\sigma }_6\right)=$$\left({\sigma }_{11},{\sigma }_{22},{\sigma }_{33},{\sigma }_{23},{\sigma }_{13},{\sigma }_{12}\right)$, $\left(S_1,S_2,S_3,S_4,S_5,S_6\right)= \left({\epsilon }_{11},{\epsilon }_{22},{\epsilon }_{33},{\epsilon }_{23},{\epsilon }_{13},{\epsilon }_{12}\right)$. The constitutive law for stress and strain is $\sigma =CS$, the stiffness matrix $C$ of $Z-cut{2}^{\circ }$ quartz is [23]

$$c^E=\left[\begin{array}{cccccc}86.74& 6.99& 11.91& -17.91& 0& 0\\ 6.99& 86.74& 11.91& 17.91& 0& 0\\ 11.91& 11.91& 107.2& 0& 0& 0\\ -17.91& 17.91& 0& 57.94& 0& 0\\ 0& 0& 0& 0& 57.94& -17.91\\ 0& 0& 0& 0& -17.91& 39.88\end{array}\right]\times {10}^9 \mathrm{N}/{\mathrm{m}}^2$$

(10)

Since the lower prong loaded with a tip-sample potential force is slender, its deformation can be modeled by the well-known Euler beam theory, that is, the cross section perpendicular to the neutral axis remains perpendicular during deformation. Its free body diagram is shown in Figure 3. Here, the conjunctive surface of the proof mass and the prong is assumed to be rotatable and has no $x_1$ component of displacement vector. The external shear force $Q_{ext}$ and moment $M_{ext}$, acting at the left boundary, are exerted by the proof mass. The parameters $L_b,b_b,\mathrm{and}h$ denote the length, width, and height of the rectangular cross-section of the prong. For the QTF of $\left(ZYw\right)+2^{\circ }$ layout [24], the $x_2-\mathrm{axis}$ is designated to be along the longitudinal direction of the prong and the $x_3-\mathrm{axis}$ is perpendicular to the QTF surface.

Let $w\left(x_2,t\right)$ denote the transverse displacement of the central line of the prong and the displacement field $\left(u_1,u_2,u_3\right)$ of any point $\left(x_1,x_2,x_3\right)$ of the prong based on the Euler beam theory is

$$u_1=w\left(x_2,t\right),u_2=-x_1{w}_{,2}\left(x_2,t\right),u_3=0$$

(11)

where $w_{,2}=\partial w/\partial x_2$. Then the strain fields are

$$S_1=0, S_2=-x_1{w}_{,22}, S_3=S_4=S_5=S_6=0$$

(12)

The velocity of point $\left(x_1,x_2,x_3\right)$is

$$v={\dot{u}}_{1}i+{\dot{u}}_{2}j+{\dot{u}}_{3}k=\frac{\partial w\left(x_2,t\right)}{\partial t}i-x_1\frac{{\partial }^{2}w\left(x_2,t\right)}{\partial x_2\partial t}j$$

(13)

The Hamilton’ principle for an elastic body is

$${\displaystyle {\int }_{t_i}^{t_f}\left[\delta T-\delta U+\delta W^{\left(m\right)}\right]dt}=0$$

(14)

where $T$ is the kinetic energy, $U$ the elastic energy, and $W^{\left(m\right)}$ the work done by the applied mechanical force. At the interface, the moment $M_{ext}$ has the contribution $M_{ext}{{\delta }\left(w_{,2}\right)|}_{x_2=0}$ to ${\delta }{M}^{\left(m\right)}$ and no contribution for the shear force $Q_{ext}$. Because there is no $x_1-\mathrm{axis}$ component of displacement at $x_2=0$, the spring force and constant force $F_0$ have done virtual work, due to the virtual displacement ${{\delta }w|}_{x_2=L_b}$. So

$${\delta }{W}^{\left(m\right)}=M_{ext}{{\delta }\left(w_{,2}\right)|}_{x_2=0}+{\left(kw+F_0\right){\delta }w|}_{x_2=L_b}$$

(15)

The first term of (14), that is, the variation of kinetic energy, can be derived as (see Appendix A)

$${\displaystyle {\int }_{t_i}^{t_f}\left[{\delta }T\right]dt}={\rho }{\displaystyle {\int }_{t_i}^{t_f}\left\{{\displaystyle {\int }_{\partial \tau }{-\left[\left(x_1^2{\ddot{w}}_{,2}\right){\delta }w\right]|}_{x_2=0}^{x_2=L_b}ds+{\displaystyle {\int }_{\tau }\left(-\ddot{w}+x_1^2{\ddot{w}}_{,22}\right){\delta }wd\tau }}\right\}}dt$$

(16)

where $\tau$ and $\partial \tau$ are the volume and boundary surface of the prong. The strain energy per unit volume is

$$\frac{1}{2}{c}_{pq}^E{S}_p{S}_q=-\frac{1}{2}{\left[\begin{array}{l}0\\ S_2\\ 0\\ 0\\ 0\\ 0\end{array}\right]}^T\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& 0& 0\\ c_{12}& c_{22}& c_{23}& c_{24}& 0& 0\\ c_{13}& c_{23}& c_{33}& c_{34}& 0& 0\\ c_{14}& c_{24}& c_{34}& c_{44}& 0& 0\\ 0& 0& 0& 0& c_{55}& c_{56}\\ 0& 0& 0& 0& c_{56}& c_{66}\end{array}\right]\left\{\begin{array}{l}0\\ S_2\\ 0\\ 0\\ 0\\ 0\end{array}\right\}=\frac{1}{2}{c}_{22}{S}_2{}^2$$

(17)

Substituting Equation (12) into the second term of Hamilton’s principle (14) gives

$$\begin{array}{l}{\displaystyle {\int }_{t_i}^{t_f}{\displaystyle {\int }_{\tau }-\delta U}d\tau }dt={\displaystyle {\int }_{t_i}^{t_f}\left\{{\displaystyle {\int }_{\tau }\delta \left[-\frac{1}{2}{c}_{pq}^E{S}_p{S}_q\right]}d\tau \right\}}dt={\displaystyle {\int }_{t_i}^{t_f}\left\{{\displaystyle {\int }_{\tau }\delta \left[-\frac{1}{2}{c}_{22}{S}_2{}^2\right]}d\tau \right\}}dt\\ ={\displaystyle {\int }_{t_i}^{t_f}\left\{{\displaystyle {\int }_{\tau }\delta \left[-\frac{1}{2}{c}_{22}{\left(-x_1{w}_{,22}\right)}^2\right]}d\tau \right\}}dt={\displaystyle {\int }_{t_i}^{t_f}\left\{{\displaystyle {\int }_{\tau }\left[-\left(c_{22}{x}_1{}^2{w}_{,22}\right)\delta w_{,22}\right]}d\tau \right\}}dt\end{array}$$

(18)

The integration by parts of the right-hand side of (18) is given by Equations (A5) and (A6). Finally we have

$$\begin{array}{l}{\displaystyle {\int }_{t_i}^{t_f}\left\{{\displaystyle {\int }_{\tau }\delta \left[-\frac{1}{2}{c}_{pq}^E{S}_p{S}_q\right]}d\tau \right\}}dt\\ ={\displaystyle {\int }_{t_i}^{t_f}{c}_{22}\left\{{{\displaystyle {\int }_{\partial \tau }\left[-\left(x_1{}^2{w}_{,22}\right)\delta w_{,2}+\left(x_1{}^2{w}_{,222}\right)\delta w\right]|}}_{x_2=0}^{x_2=L_b}dS+{\displaystyle {\int }_{\tau }\left[-\left(x_1{}^2{w}_{,2222}\right)\delta w\right]}d\tau \right\}}dt\end{array}$$

(19)

The substitution of Equations (15), (16), and (19) into the Hamilton’s principle (14) results in

$$\begin{array}{l}{\displaystyle {\int }_{t_i}^{t_f}\{{\displaystyle {\int }_{t_i}^{t_f}\text{{}{\displaystyle {\int }_{\tau }\left(e{q}_1\delta w\right)}d\tau }+{\displaystyle {\int }_{\partial \tau }\left[{\left(B_1\delta w+B_2\delta w_{,2}\right)|}_{x_2=0}^{x_2=L_b}\right]}dS}\\ +{\left(M_{ext}\delta w_{,2}\right)|}_{x_2=0}+{\left(kw+F_0\right)\delta w|}_{x_2=L_b}\}dt=0\end{array}$$

(20)

The variation $\delta w$ in the volume integration of the first term of Equation (20) is arbitrary and its coefficient $e{q}_1$ must be zero, this leads to the equation of motion of the lower loaded prong as

$$\rho \left(-\ddot{w}+x_1{}^2{\ddot{w}}_{,22}\right)-c_{22}\left(x_1{}^2{w}_{,2222}\right)=0$$

(21)

Based on (20), at the left boundary $x_2=0$, ${\delta w|}_{x_2}=0$ leads to $B_1\ne 0$ and ${\delta w_{,2}|}_{x_2=0}\ne 0$ leads to ${\left(B_2+M_{ext}\right)|}_{x_2=0}=0$, which is the left boundary condition

$$-c_{22}\left(x_1{}^2{w}_{,22}\right)+M_{ext}=0,\mathrm{at}{x}_2=0$$

(22a)

At the right boundary $x_2=L_b$, both $\delta w$ and $\delta w_{,2}$ are not zero, therefore, their coefficients must be zero, which leads to the right boundary conditions as

$$\{\begin{array}{c}-\rho \left(x_1{}^2{\ddot{w}}_{,2}\right)+c_{22}\left(x_1{}^2{w}_{,222}\right)+\left(kw+F_0\right)=0,\\ -c_{22}\left(x_1{}^2{w}_{,22}\right)=0.\end{array}\mathrm{at}{x}_2=L_b$$

(22b)

For static equilibrium, all the velocities and accelerations are zero. By setting $\ddot{w}=0$ and ${\ddot{w}}_{,22}=0$, Equations (21) and (22a,b) become the static equation and boundary conditions for static displacement $w^{\ast }\left(x_2\right)$ as

$$c_{22}\left(x_1{}^2{w}^{\ast }{}_{,2222}\right)=0.$$

(23)

$$M^{\ast }{}_{ext}=c_{22}\left(x_1{}^2{w}^{\ast }{}_{,22}\right),\mathrm{at}{x}_2=0.$$

(24)

$$\{\begin{array}{c}{c}_{22}\left(x_1{}^2{w}^{\ast }{}_{,222}\right)+k{w}^{\ast }+F_0=0,\\ -c_{22}\left(x_1{}^2{w}^{\ast }{}_{,22}\right)=0.\end{array}\mathrm{at}{x}_2=L_b$$

(25)

Substituting (3) into (22a,b) and subtracting (23) from it gives

$$\rho \left(-\ddot{\overline{w}}+x_1{}^2{\ddot{\overline{w}}}_{,22}\right)-c_{22}\left(x_1{}^2{\overline{w}}_{,2222}\right)=0$$

(26)

and

$$\{\begin{array}{}\mathrm{(27a)}& {\overline{M}}_{ext}=c_{22}\left(x_1{}^2{\overline{w}}_{,22}\right),& \mathrm{at}{x}_2=0,\mathrm{(27b)}& -\rho \left(x_1{}^2{\ddot{\overline{w}}}_{,2}\right)+c_{22}\left(x_1{}^2{\overline{w}}_{,222}\right)+k\overline{w}=0,& \mathrm{at}{x}_2=L_b\mathrm{(27c)}& -c_{22}\left(x_1{}^2{\overline{w}}_{,22}\right)=0.& \mathrm{at}{x}_2=L_b\end{array}$$

where $M_{ext}^{}=M_{ext}^{\ast }\left(x_2\right)+{\overline{M}}_{ext}\left(x_2,t\right)$. The boundary condition (27c) is independent of the constant force $F_0$. It is clear that the natural frequencies of vibration around the static equilibrium state are not affected by the constant $F_0$, but influenced by the spring constant $k$.

2.3. The General Solution of the Lower Loaded Prong

The resonant solution of system (24) is of the form

$$\overline{w}\left(x_2,t\right)=W\left(x_2\right)e^{i\omega t}$$

(28)

where $\omega$ is the natural frequency and $W\left(x_2\right)$ is the mode shape. Substituting (28) into (26) and integrating over the cross-section area in the $x_1-x_3$ plane yields

$$c_{22}{I}_2{W}_{,2222}+\rho {\omega }^2{I}_2{W}_{,22}-\rho {\omega }^2{I}_{1}W=0$$

(29)

where

$$\begin{array}{c}{I}_1={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\displaystyle {\int }_{-\frac{b_b}{2}}^{\frac{b_b}{2}}}}d{x}_{1}d{x}_3=b_{b}h,\end{array}$$

(30a)

$$\begin{array}{c}{I}_2={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\displaystyle {\int }_{-\frac{b_b}{2}}^{\frac{b_b}{2}}{x}_1{}^2}}d{x}_{1}d{x}_3=\frac{b_b^{3}h}{12}.\end{array}$$

(30b)

The mode shape is assumed to be of the form $W\left(x_2\right)={\mathrm{Re}}^{\lambda x_2}$, and substituting this into (29) gives the characteristic equation

$$c_{22}{I}_2{\lambda }^4+\rho {\omega }^2{I}_2{\lambda }^2-\rho {\omega }^2{I}_1=0$$

(31)

The solution of (31) is

$${\lambda }^2=\frac{-\rho {\omega }^2{I}_2\pm \sqrt{{\left(\rho {\omega }^2{I}_2\right)}^2+4\left(c_{22}{I}_2\right)\left(\rho {\omega }^2{I}_1\right)}}{2\left(c_{22}{I}_2\right)}$$

(32)

The general solution of (29) is

$$W\left(x_2\right)=R_1\cosh {\xi }_1{x}_2+R_2\sinh {\xi }_1{x}_2+R_3\cos {\xi }_2{x}_2+R_4\sin {\xi }_2{x}_2$$

(33)

where ${\xi }_1$ and ${\xi }_2$ are related to the natural frequency $\omega$ by

$$\begin{array}{c}{\xi }_1=\sqrt{\frac{-\rho {\omega }^2{I}_2+\sqrt{\mathsf{\Gamma }}}{2\left(c_{22}{I}_2\right)}},\end{array}$$

(34a)

$$\begin{array}{c}{\xi }_2=\sqrt{\frac{\rho {\omega }^2{I}_2+\sqrt{\mathsf{\Gamma }}}{2\left(c_{22}{I}_2\right)}}\end{array}$$

(34b)

and

$$\mathsf{\Gamma }={\left(\rho {\omega }^2{I}_2\right)}^2+4\left(c_{22}{I}_2\right)\left(\rho {\omega }^2{I}_1\right)$$

(35)

3. The Deformation and Equations of Motion of the Proof Mass

3.1. The Deformation of the Proof Mass

When the QTF is in the state of free vibration at its in-plane anti-phase natural frequency, the two prongs vibrate in opposite directions with the same magnitude. These two prongs will exert shear force and bending moment of opposite direction and equal magnitude on the proof mass, as shown in Figure 4. When the QTF is put closely above the tested sample, the tip-sample interaction force will exert additional moment and shear to the proof mass at their conjunctive area. The net shear force and moment in the interface is shown in Figure 5, and they satisfy

$$Q_2=Q_1+F,M_2=M_1+M_F,$$

(36)

where $M_F=F{L}_b$, $F$ is the spring force acting on the tip of the lower prong.

3.1.1. Deformation Due to the Bending Moment of the Tip-Sample Interaction

The deformation of the proof mass subject to the bending moment $M_F$ could be of bending type or of pure shearing type. We have adopted the type of pure shear deformation in this paper, as shown in Figure 6. The cross-section is inclined by an angle $\beta \left(x_2,t\right)$.

The displacement of this pure shearing is

$$u_1^{\prime }=0,u_2^{\prime }=-x_1\beta ,u_3^{\prime }=0.$$

(37)

3.1.2. The Wrapping Deformation of the Proof Mass

The cross-section experiences wrapped deformation in the anti-phase mode, as shown in Figure 7. We assume that the rotation angle $\theta$ of any point $\left(x_1,x_2\right)$ in the cross-section is linearly proportional to the distance $x_1$, that is, $\theta =x_1\psi \left(x_2,t\right)$. We assume that there is no axial displacement at the conjunctive point $x_1=d/2$, i.e., $u_2^{″}\left(d/2,x_2,t\right)=0.$ Then we have

$$u_2^{″}\left(x_1,x_2,t\right)-u_2^{″}\left(d/2,x_2,t\right)={\displaystyle {\int }_{\frac{d}{2}}^{x_1}{x}_1\psi \left(x_2,t\right)d{x}_1}=\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right).$$

So the displacement field of the wrapping deformation is

$$u_1^{″}=0,u_2^{″}=\left(\frac{1}{2}{x}_1{}^2-\frac{d^2}{8}\right)\psi \left(x_2,t\right),u_3^{″}=0.$$

(38)

By rotating the wrapping displacement field around the $x_3$ axis by an angle $-\beta \left(x_2,t\right)$, the components resolved along the $x_1-x_2-x_3$ coordinates are obtained as

$$\left\{\begin{array}{c}{u}_1^{\left(w\right)}\\ \begin{array}{c}{u}_2^{\left(w\right)}\\ u_3^{\left(w\right)}\end{array}\end{array}\right\}=\left[\begin{array}{cc}\cos \beta & \begin{array}{cc}\sin \beta & 0\end{array}\\ \begin{array}{c}-\sin \beta \\ 0\end{array}& \begin{array}{cc}\begin{array}{c}\cos \beta \\ 0\end{array}& \begin{array}{c}0\\ 1\end{array}\end{array}\end{array}\right]\left\{\begin{array}{c}{u}_1^{″}\\ \begin{array}{c}{u}_2^{″}\\ u_3^{″}\end{array}\end{array}\right\}=\left\{\begin{array}{c}\psi \left(x_2,t\right)\sin \beta \left(x_2,t\right)\left[\frac{x_1^2}{2}-\frac{d^2}{8}\right]\\ \cos \beta \left(x_2,t\right)\psi \left(x_2,t\right)\left[\frac{x_1^2}{2}-\frac{d^2}{8}\right]\\ 0\end{array}\right\}$$

(39)

If we assume a small deformation, thus $\sin \beta \approx \beta$ and $\cos \beta \approx 1$, and drop the nonlinear term, we get

$$u_1^{\left(w\right)}=0,u_2^{\left(w\right)}=\psi \left(x_2,t\right)\left(\frac{x_1^2}{2}-\frac{d^2}{8}\right),u_3^{\left(w\right)}=0.$$

(40)

Equation (38) plus equation (40) gives the total displacement of the proof mass

$$u_1=0,u_2=-x_1\beta \left(x_2,t\right)+\left(\frac{x_1^2}{2}-\frac{d^2}{8}\right)\psi \left(x_2,t\right),u_3=0$$

(41)

The strain field of the proof mass is

$$\begin{array}{ccc}{S}_1=0,& S_2=-x_1{\beta }_{,2}+\left(\frac{x_1^2}{2}-\frac{d^2}{8}\right){\psi }_{,2},& S_3=0,\\ S_4=0,& S_5=0,& S_6=-\beta +x_1\psi .\end{array}$$

(42)

The strain energy per unit volume is

$$\frac{1}{2}{c}_{pq}^E{S}_p{S}_q=\frac{1}{2}{\left[\begin{array}{l}0\\ S_2\\ 0\\ 0\\ 0\\ S_6\end{array}\right]}^T\left[\begin{array}{cccccc}{c}_{11}& c_{12}& c_{13}& c_{14}& 0& 0\\ c_{12}& c_{22}& c_{23}& c_{24}& 0& 0\\ c_{13}& c_{23}& c_{33}& c_{34}& 0& 0\\ c_{14}& c_{24}& c_{34}& c_{44}& 0& 0\\ 0& 0& 0& 0& c_{55}& c_{56}\\ 0& 0& 0& 0& c_{56}& c_{66}\end{array}\right]\left\{\begin{array}{l}0\\ S_2\\ 0\\ 0\\ 0\\ S_6\end{array}\right\}=\frac{1}{2}\left(c_{22}{S}_2{}^2+c_{66}{S}_6{}^2\right)$$

(43)

3.2. Equations of Motion of the Proof Mass

The first term of the Hamilton’s principle (14) with the use of (41) is carried out as

$$\begin{array}{l}{\displaystyle {\int }_{t_i}^{t_f}\left\{\delta T\right\}dt={\displaystyle {\int }_{\tau }\rho \left\{-{\displaystyle {\int }_{t_i}^{t_f}{\ddot{U}}_i\delta U_{i}dt}\right\}d\tau }}\\ ={\displaystyle {\int }_{\tau }\rho \left\{-{\displaystyle {\int }_{t_i}^{t_f}\left[{\ddot{U}}_1\delta U_1+{\ddot{U}}_2\delta U_2+{\ddot{U}}_3\delta U_3\right]dt}\right\}d\tau }\\ ={\displaystyle {\int }_{\tau }\rho \{-{\displaystyle {\int }_{t_i}^{t_f}\left[-x_1\ddot{\beta }\left(x_2,t\right)+\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)\ddot{\psi }\left(x_2,t\right)\right]}}\\ \times \left[-x_1\delta \beta \left(x_2,t\right)+\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)\delta \psi \left(x_2,t\right)\right]dt\}d\tau \end{array}$$

(44)

The second term of Hamilton’s principle with the use of (43) is developed out as

$$\begin{gathered} \begin{array}{l}{\displaystyle {\int }_{t_i}^{t_f}{\displaystyle {\int }_{\tau }-\delta U}d\tau }dt=={\displaystyle {\int }_{t_i}^{t_f}\left\{{\displaystyle {\int }_{\tau }\delta \left[-\frac{1}{2}{c}_{pq}^E{S}_p{S}_q\right]}d\tau \right\}}dt=\\ {\displaystyle {\int }_{t_i}^{t_f}\{{\displaystyle {\int }_{\partial \tau }[\left(-c_{22}{\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)}^2{\psi }_{,2}+c_{22}\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)x_1{\beta }_{,2}\right)\delta \psi }}\\ {+\left(c_{22}\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)x_1{\psi }_{,2}-c_{22}{x}_1{}^2{\beta }_{,2}\right)\delta \beta ]}_{x_2=0}^{x_2=L_m}dS\end{array} \\ \begin{array}{l}+{\displaystyle {\int }_{\tau }[\left(c_{22}{\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)}^2{\psi }_{,22}-c_{22}\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)x_1{\beta }_{,22}-c_{66}\left(x_1\psi -\beta \right)x_1\right)\delta \psi }\\ +\left(-\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)c_{22}{x}_1{\psi }_{,22}+c_{22}{x}_1{}^2{\beta }_{,22}+c_{66}\left(x_1\psi -\beta \right)\right)\delta \beta ]d\tau \}dt\end{array} \end{gathered}$$

(45)

The virtual work done on the proof mass has two parts. One comes from the product of the moment $M_1$ as shown in Figure 4 and the virtual rotation $\delta \theta =x_1\delta \psi$ of the cross-section as shown in Figure 7 at $x_1=\pm d/2$. The other comes from the product of the moment $F{L}_b$ caused by the interaction force and the virtual rotation $\delta \beta$ of the cross-section due to the pure shearing deformation, that is,

$$\begin{array}{l}\delta W^{(m)}=\delta W_{\mathrm{int}erface}^{\left(m\right)}=\delta W_{warping}^{\left(m\right)}+\delta W_{\mathrm{Inclination}}^{\left(m\right)}\\ =\left(\left(-M_1\right)\delta \left(-\frac{d}{2}\psi \right)|{}_{x_2=L_m,,x_1=-\frac{d}{2}}+M_1\delta \left(\frac{d}{2}\psi \right)|{}_{x_2=L_m,x_1=\frac{d}{2}}\right)+F{L}_b\delta \beta |{}_{x_2=L_m}\\ =\left[2M_1\delta \left(\frac{d}{2}\psi \right)+F{L}_b\delta \beta \right]|{}_{x_2=L_m}\end{array}$$

(46)

where $L_m$ and $b_m$ are the length and width of the proof mass, respectively. Substituting Equations (44)–(46) into the Hamilton’s principle (14) yields

$$\begin{array}{l}{\displaystyle {\int }_{t_i}^{t_f}\{{\displaystyle {\int }_{\tau }\left(\overline{e}{\overline{q}}_1\cdot \delta \psi +\overline{e}{\overline{q}}_2\cdot \delta \beta \right)d\tau }}\\ +{\displaystyle {\int }_{\partial \tau }\left[{\left({\overline{B}}_1\cdot \delta \psi +{\overline{B}}_2\cdot \delta \beta \right)|}_{x_2=0}^{x_2=L_m}\right]}dS+\left[M_{1}d\delta \psi +F{L}_b\delta \beta \right]|{}_{x_2=L_m}\}dt=0\end{array}$$

(47)

In the volume integration of (47) the variations $\delta \psi$ and $\delta \beta$ are arbitrary, that their coefficients $\overline{e{q}_1}$ and $\overline{e{q}_2}$ must be zero gives the two equations of motion as

$$\begin{array}{l}{c}_{22}{\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)}^2{\psi }_{,22}-c_{22}\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)x_1{\beta }_{,22}-c_{66}\left(x_1\psi -\beta \right)x_1\\ +\rho x_1\left(\frac{1}{2}{x}_1{}^2-\frac{d^2}{8}\right)\ddot{\beta }\left(x_2,t\right)-\rho {\left(\frac{1}{2}{x}_1{}^2-\frac{d^2}{8}\right)}^2\ddot{\psi }\left(x_2,t\right)=0,\end{array}$$

(48)

$$\begin{array}{l}-\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)c_{22}{x}_1{\psi }_{,22}+c_{22}{x}_1{}^2{\beta }_{,22}+c_{66}\left(x_1\psi -\beta \right)\\ -\rho x_1{}^2\ddot{\beta }\left(x_2,t\right)+\rho x_1\left(\frac{1}{2}{x}_1{}^2-\frac{d^2}{8}\right)\ddot{\psi }\left(x_2,t\right)=0\end{array}$$

(49)

The left boundary conditions of the proof mass are

$$\psi =0,\beta =0,\mathrm{at}{x}_2=0.$$

(50)

The natural boundary at the right end of the proof mass can be derived from (47).

The variations $\delta \psi$ and $\delta \beta$ at $x_2=L_m$ are arbitrary, so their coefficients $\left({\overline{B}}_1+M_{1}d\right)$ and $\left({\overline{B}}_2+F{L}_b\right)$ must be zero, this leads to right boundary conditions

$$\{\begin{array}{}\mathrm{(51a)}& -c_{22}{\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)}^2{\psi }_{,2}+c_{22}\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)x_1{\beta }_{,2}+M_{1}d=0,\mathrm{(51b)}& c_{22}\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)x_1{\psi }_{,2}-c_{22}{x}_1{}^2{\beta }_{,2}+F{L}_b=0.\end{array}$$

The integration of Equation (51a,b) over the cross-section area in $x_1-x_3^{}$ plane gives

$$\{\begin{array}{c}-c_{22}{\overline{I}}_1{\psi }_{,2}+{\overline{I}}_5{M}_{1}d=0,\\ -c_{22}{\overline{I}}_3{\beta }_{,2}+{\overline{I}}_{5}F{L}_b=0.\end{array}\mathrm{at}{x}_2=L_m$$

(52)

3.3. General Solution of the Proof Mass

Integrating the Equations (48) and (49) of motion over the cross-section area gives

$$\{\begin{array}{c}{c}_{22}{\overline{I}}_1{\psi }_{,22}-c_{66}{\overline{I}}_3\psi -\rho {\overline{I}}_1\ddot{\psi }=0,\\ c_{22}{\overline{I}}_3{\beta }_{,22}-c_{66}{\overline{I}}_5\beta -\rho {\overline{I}}_3\ddot{\beta }=0,\end{array}$$

(53)

where

$${\overline{I}}_1={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\displaystyle {\int }_{-\frac{b_m}{2}}^{\frac{b_m}{2}}{\left(\frac{1}{2}{x}_1^2-\frac{d^2}{8}\right)}^2}}d{x}_{1}d{x}_3=\frac{h}{960}\left(3b_m{}^5-10b_m{}^3{d}^2+15b_m{d}^4\right),$$

(54a)

$$\begin{array}{c}{\overline{I}}_3={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\displaystyle {\int }_{-\frac{b_m}{2}}^{\frac{b_m}{2}}{x}_1^2}}d{x}_{1}d{x}_3=\frac{b_m^{3}h}{12},\end{array}$$

(54b)

$$\begin{array}{c}{\overline{I}}_5={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\displaystyle {\int }_{-\frac{b_m}{2}}^{\frac{b_m}{2}}}}d{x}_{1}d{x}_3=b_{m}h.\end{array}$$

(54c)

Equations (52) and (53) reveal that variables $\psi$ and $\beta$ are uncoupled, due to the fact that the nonlinear coupling term in the displacement field (39) is neglected for small amplitude vibration. The resonant solution is of the form

$$\begin{array}{c}\psi \left(x_2,t\right)=\mathsf{\Psi }\left(x_2\right)e^{i\omega t},\end{array}$$

(55a)

$$\begin{array}{l}\beta \left(x_2,t\right)=\mathsf{\Theta }\left(x_2\right)e^{i\omega t}.\end{array}$$

(55b)

Substituting (55) into (53) gives

$$\{\begin{array}{}\mathrm{(56a)}& c_{22}{\overline{I}}_1{\mathsf{\Psi }}_{,22}-c_{66}{\overline{I}}_3\mathsf{\Psi }+\rho {\omega }^2{\overline{I}}_1\mathsf{\Psi }=0,\mathrm{(56b)}& c_{22}{\overline{I}}_3{\mathsf{\Theta }}_{,22}-c_{66}{\overline{I}}_5\mathsf{\Theta }+\rho {\omega }^2{\overline{I}}_3\mathsf{\Theta }=0.\end{array}$$

The solution of the mode shape for wrapping vibration is of the form $\mathsf{\Psi }\left(x_2\right)=\overline{R}{e}^{\lambda x_2}$. Substituting this into (56a) gives

$$c_{22}{\overline{I}}_1{\lambda }_1{}^2-\left(c_{66}{\overline{I}}_3-\rho {\omega }^2{\overline{I}}_1\right)=0$$

(57)

The solution is ${\lambda }_1{}^2=\left(c_{66}{\overline{I}}_3-\rho {\omega }^2{\overline{I}}_1\right)/c_{22}{\overline{I}}_1.$ The mode shape is

$$\mathsf{\Psi }\left(x_2\right)={\overline{R}}_1\cosh {\xi }_3{x}_2+{\overline{R}}_2\sinh {\xi }_3{x}_2$$

(58)

where

$${\xi }_3=\sqrt{\frac{-\left(c_{66}{\overline{I}}_3-\rho {\omega }^2{\overline{I}}_1\right)}{c_{22}{\overline{I}}_1}}$$

(59)

Similarly, the mode shape of Equation (56b) is

$$\mathsf{\Theta }\left(x_2\right)={\overline{R}}_3\cosh {\xi }_4{x}_2+{\overline{R}}_4\sinh {\xi }_4{x}_2,$$

(60)

where ${\xi }_4$ is related to the natural frequency $\omega$ by

$${\xi }_4=\sqrt{\frac{-\left(c_{66}{\overline{I}}_5-\rho {\omega }^2{\overline{I}}_3\right)}{c_{22}{\overline{I}}_3}}$$

(61)

4. Anti-Phase Resonant Frequency of QTF

In this section, the three components of the QTF, one free prong, one loaded prong, and the proof mass, are considered as a whole. Their individual general solutions were obtained in previous sections, with the twelve coefficients of mode shapes and the resonant frequency left unknown, to be determined. When the twelve boundary and interface conditions are imposed, the anti-phase resonant frequency can be evaluated. Also, the frequency shift corresponding to the spring constant $k$ (or the negative force gradient) can be found. Three different sets of coordinates need to be used as shown in Figure 8.

4.1. Mode Shapes with Undetermined Coefficients

The mode shape of the upper free prong is almost similar to Equation (33)

$$W^{(2)}\left(x_2^{\left(2\right)}\right)=R_1{}^{(2)}\cosh {\xi }_1{x}_2^{\left(2\right)}+R_2{}^{(2)}\sinh {\xi }_1{x}_2^{\left(2\right)}+R_3{}^{(2)}\cos {\xi }_2{x}_2^{\left(2\right)}+R_4{}^{(2)}\sin {\xi }_2{x}_2^{\left(2\right)}$$

(62)

where ${\xi }_1$ and ${\xi }_2$ are functions of resonant frequency implicitly as

$${\xi }_1=\sqrt{\frac{-\rho {\omega }^2{I}_2+\sqrt{\mathsf{\Gamma }}}{2\left(c_{22}{I}_2\right)}}, {\xi }_2=\pm \sqrt{\frac{\rho {\omega }^2{I}_2+\sqrt{\mathsf{\Gamma }}}{2\left(c_{22}{I}_2\right)}}$$

(63)

where $\mathsf{\Gamma }$ is given by (35), coefficients $\left(R_1^{\left(2\right)},R_2^{\left(2\right)},R_3^{\left(2\right)},R_4^{\left(2\right)}\right)$ are to be determined.

The boundary conditions at the right free end are the same as Equation (27b,c), but with the interaction force $k\overline{w}$ and external moment $M_{ext}$ removed.

$$M_{ext}=c_{22}{\left(x_1^{\left(2\right)}\right)}^2{w}_{,22}^{\left(2\right)},\mathrm{at}{x}_2^{\left(2\right)}=0$$

(64a)

$$\{\begin{array}{}\mathrm{(64b)}& -\rho \left[{\left(x_1^{\left(2\right)}\right)}^2{\ddot{w}}_{,2}^{\left(2\right)}\right]+c_{22}\left[{\left(x_1^{\left(2\right)}\right)}^2{w}_{,222}^{\left(2\right)}\right]=0,\mathrm{at}{x}_2^{\left(2\right)}=L_b\mathrm{(64c)}& -c_{22}\left[{\left(x_1^{\left(2\right)}\right)}^2{w}_{,22}^{\left(2\right)}\right]=0,\mathrm{at}{x}_2^{\left(2\right)}=L_b\end{array}$$

Equation (64b) is the shear-force free condition and (64c) is the moment free condition.

Integrating (64b,c) over the cross-section area and using the mode shape (62) gives

$$\{\begin{array}{c}\rho {\omega }^2{I}_2{W}^{\left(2\right)}{}_{,2}+c_{22}{I}_2{W}^{\left(2\right)}{}_{,222}=0,\\ -c_{22}{I}_2{W}^{\left(2\right)}{}_{,22}=0.\end{array}\mathrm{at}{x}_2^{\left(2\right)}=L_b$$

(65)

The mode shapes of the proof mass for cross-section wrapping and shearing in terms of the newly-defined coordinates are

$$\{\begin{array}{c}{\mathsf{\Psi }}^{(1)}\left(x_2^{\left(1\right)}\right)={\overline{R}}_1{}^{(1)}\cosh {\xi }_3{x}_2^{\left(1\right)}+{\overline{R}}_2{}^{(1)}\sinh {\xi }_3{x}_2^{\left(1\right)},\\ {\mathsf{\Theta }}^{\left(1\right)}\left(x_2^{\left(1\right)}\right)={\overline{R}}_3{}^{(1)}\cosh {\xi }_4{x}_2^{\left(1\right)}+{\overline{R}}_4{}^{(1)}\sinh {\xi }_4{x}_2^{\left(1\right)}.\end{array}$$

(66)

where ${\xi }_3$ and ${\xi }_4$, being functions of resonant frequency, are given by (59) and (61).

The boundary and interface conditions are given by (50) and (52).

The mode shape of the lower loaded prong in terms of the new coordinates is

$${\overline{W}}^{(3)}(x_2^{\left(3\right)})=R_1{}^{(3)}\cosh {\xi }_1{x}_2^{\left(3\right)}+R_2{}^{(3)}\sinh {\xi }_1{x}_2^{\left(3\right)}+R_3{}^{(3)}\cos {\xi }_2{x}_2^{\left(3\right)}+R_4{}^{(3)}\sin {\xi }_2{x}_2^{\left(3\right)}$$

(67)

where the frequency dependent parameters ${\xi }_1$ and ${\xi }_2$ are given by Equation (34a,b). The boundary conditions are

$${\overline{M}}_{ext}=c_{22}\left(x_1^{\left(3\right)}\right){\overline{w}}_{,22}^{\left(3\right)},\mathrm{at}{x}_2^{\left(3\right)}=0.$$

(68a)

Integrating Equations (27b,c) over the cross-section area gives

$$\{\begin{array}{c}\rho {\omega }^2{I}_2{\overline{W}}_{,2}^{\left(3\right)}+c_{22}{I}_2{W}_{,222}^{\left(3\right)}+k{b}_{b}h{\overline{W}}^{\left(3\right)}=0,\\ -c_{22}{I}_2{W}_{,22}^{\left(3\right)}=0.\end{array}\mathrm{at}{x}_2^{\left(3\right)}=L_b$$

(68b)

4.2. The Twelve Boundary and Interface Conditions

The proof mass has no wrapping rotation and shear inclination at the clamped boundary, so

$${\mathsf{\Psi }}^{\left(1\right)}\left(x_2^{\left(1\right)}=0\right)=0,$$

(69a)

$${\mathsf{\Theta }}^{\left(1\right)}\left(x_2^{\left(1\right)}=0\right)=0$$

(69b)

At the interface, the free prong has no vertical displacement,

$$W^{\left(2\right)}\left(x_2^{\left(2\right)}=0\right)=0.$$

(69c)

The slope of the proof mass at the interface is equal to that of the free prong, i.e.,

$${\left({\mathsf{\Theta }}^{\left(1\right)}-\frac{d}{2}{\mathsf{\Psi }}^{\left(1\right)}\right)|}_{x_2^{\left(1\right)}=L_m,x_1^{\left(1\right)}=-d/2}=W_{,2}^{\left(2\right)}\left(x_2^{\left(2\right)}=0\right)$$

(69d)

The condition that the moment acting on the proof mass is equal in magnitude to the moment acting on the free prong at the upper interface is not used, because the end surface of the prong is flat while the proof mass is curvilinear. Instead, the condition that the work, $\delta W_{{}_{m2}}^{\left(1\right)}$, done by the moment $M_1$ in (51a) on the virtual angular displacement of the proof mass be equal to that, $\delta W_b^{(2)}$, done by the moment $M_{ext}\left(x_2^{\left(2\right)}\right)$ in (64a) on the virtual rotation of the prong end face is used here. Note that

$$\delta W_b^{(2)}={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\displaystyle {\int }_{-\frac{b_b}{2}}^{\frac{b_b}{2}}\left\{\left[-c_{22}\left({\left(x_1{}^{\left(2\right)}\right)}^2{w}_{,22}\right)\right]\delta w_{,2}\right\}}}d{x}_1^{\left(2\right)}d{x}_3^{\left(2\right)}$$

and

$$\begin{array}{l}\delta W_{{}_{m2}}^{\left(1\right)}={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\displaystyle {\int }_{-\frac{d}{2}-\frac{b_b}{2}}^{-\frac{d}{2}+\frac{b_b}{2}}\{\frac{1}{d}[c_{22}{\left(\frac{1}{2}{\left(x_1^{(1)}\right)}^2-\frac{d^2}{8}\right)}^2{\psi }_{,2}}}\\ -c_{22}\left(\frac{1}{2}{\left(x_1^{(1)}\right)}^2-\frac{d^2}{8}\right)x_1^{(1)}{\beta }_{,2}]\delta \left(-\frac{d}{2}\psi \right)\}d{x}_1^{\left(1\right)}d{x}_2^{\left(1\right)}\end{array}$$

The condition $\delta W_b^{(2)}=\delta W_{{}_{m2}}^{\left(1\right)}$ with the use of $\psi =w_{,2}$ and integration over the contact area leads to

$$d{I}_2{W}_{,22}\left(x_2{}^{(2)}=0\right)=-{\overline{I}}_1{\mathsf{\Psi }}_{,2}\left(x_2{}^{(1)}=L_m\right)$$

(69e)

where

• $I_2={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\displaystyle {\int }_{-\frac{b_b}{2}}^{\frac{b_b}{2}}\left(x_1{}^2\right)}}d{x}_{1}d{x}_3=\frac{1}{12}{b}_b{}^{3}h$,
• ${\overline{I}}_1={\displaystyle {\int }_{-\frac{h}{2}}^{\frac{h}{2}}{\displaystyle {\int }_{-\frac{d}{2}-\frac{b_b}{2}}^{-\frac{d}{2}+\frac{b_b}{2}}\left[{\left(\frac{1}{2}{\left(x_1^{(2)}\right)}^2-\frac{d^2}{8}\right)}^2\right]}}d{x}_1^{(2)}d{x}_3^{(2)}=\frac{h}{960}\left(3b_b{}^5+20b_b{}^3{d}^2\right)$

The shear force at zero at the free end of the upper prong is

$$\rho {\omega }^2{I}_2{W}_{,2}^{\left(2\right)}+c_{22}{I}_2{W}_{,222}^{(2)}=0,\mathrm{at}{x}_2^{\left(2\right)}=L_b.$$

(69f)

The bending moment at zero at the free end of the upper prong is

$$-c_{22}{I}_2{W}_{,22}^{\left(2\right)}=0.\mathrm{at}{x}_2^{\left(2\right)}=L_b.$$

(69g)

The lower prong has no displacement in the $x_1^{\left(3\right)}$ direction at the interface

$$W^{\left(3\right)}\left(x_2^{\left(3\right)}=0\right)=0.$$

(69h)

The slope of the lower prong is equal to that of the proof mass at the lower interface

$$W_{,2}^{\left(3\right)}\left(x_2{}^{(3)}=0\right)=\mathsf{\Theta }+\frac{d}{2}{\mathsf{\Psi }}^{(1)}\left(x_2{}^{(1)}=L_m,x_1{}^{(3)}=\frac{d}{2}\right)$$

(69i)

The virtual work $\delta W_b^{\left(3\right)}$ done by the moment in (24) on the virtual rotation of the prong is equal to the work $\delta W_{m3}^{\left(1\right)}$ done by the moment in (51a) on the virtual angular displacement of the proof mass, i.e., the condition $\delta W_b^{\left(3\right)}=\delta W_{m3}^{\left(1\right)}$ gives

$$d{I}_2{W}_{,22}\left(x_2{}^{(3)}=0\right)={\overline{I}}_1{\mathsf{\Psi }}_{,2}\left(x_2{}^{(1)}=L_m\right)$$

(69j)

The balance of the shear force at the right end of the lower loaded prong is

$$\rho {\omega }^2{I}_2{W}_{,2}^{(3)}+c_{22}{I}_2{W}_{,222}^{(3)}+k{b}_{b}h{W}_{}^{\left(3\right)}=0,\mathrm{at}{x}_2^{\left(3\right)}=L_b$$

(69k)

The moment-free of the right loaded prong is

$$-c_{22}{I}_2{W}_{{}_{,22}}^{(3)}=0.$$

(69l)

Substituting the general solutions (62), (66), and (67) into the twelve boundary and interface conditions given by the Equations (69a–l) yields the $12\times 12$ matrix-vector equation

$${\left[M_{ij}\right]}_{12\times 12}\left\{R_j\right\}=0$$

(70)

where $\left\{R\right\}=\left[\begin{array}{cccc}{\overline{R}}_1{}^{(1)}& \begin{array}{ccc}{\overline{R}}_2{}^{(1)}& {\overline{R}}_3{}^{(1)}& {\overline{R}}_4{}^{(1)}\end{array}& R_1{}^{(2)}& R_2{}^{(2)}\end{array}\begin{array}{cccc}{R}_3{}^{(2)}& R_4{}^{(2)}& R_1{}^{(3)}& \begin{array}{ccc}{R}_2{}^{(3)}& R_3{}^{(3)}& R_3{}^{(3)}\end{array}\end{array}\right]$.

The non-zero elements of the matrix $\left[M_{ij}\right]$ are given in Appendix B.

For the nontrivial solution of the parametric vector $\left\{R\right\}$ the determinant of the coefficient matrix of Equation (70) must be zero, that is,

$$\det \left[M_{ij}\left(\omega \right)\right]=0$$

(71)

this leads to the characteristic equation, which is a nonlinear algebraic equation of natural frequency $\omega$. By solving this characteristic equation numerically, we get the natural frequencies of the in-plane anti-phase resonant modes.

4.3. Results and Discussions

The geometric parameters of QTF are as follows. (1) The length, width, and thickness of the prong are $154\mathsf{\mu }\mathrm{m}$, $14\mathsf{\mu }\mathrm{m}$, and $7\mathsf{\mu }\mathrm{m}$, respectively. (2) The length, width, and thickness of the proof mass are $60\mathsf{\mu }\mathrm{m}$, $46\mathsf{\mu }\mathrm{m}$, and $7\mathsf{\mu }\mathrm{m}$, respectively. The gap between the two prongs is $8\mathsf{\mu }\mathrm{m}$. The Hamaker constant for quartz is $8.83\times {10}^{-20}\mathrm{J}$ [25], the radius of the tip is $8\mathrm{nm}$, the density of QTF is $\rho =2649\mathrm{kg}/\mathrm{m}$and the elastic stiffness coefficients are given by (10).

The natural frequency and frequency shift associated with various values of the spring constant are listed in

Table 1. Spring constants and natural frequencies of the first anti-phase mode and frequency shift.

$\mathit{k}=\frac{\mathbf{n}\mathbf{N}}{\mathbf{n}\mathbf{m}}$	Frequency (kHz)	Frequency Shift (Hz)	Frequency Shift Ratio (ppm)
1	451.219	21	46.5
2	451.24	42	93.1
3	451.262	64	141.8
4	451.284	86	190.6
5	451.307	109	241.5
6	451.33	132	292.5
7	451.353	155	343.4
8	451.377	179	396.6
9	451.402	204	451.9
10	451.426	228	505.1
11	451.452	254	562.6
12	451.477	279	618.0
13	451.504	306	677.7
14	451.53	332	735.3
15	451.558	360	797.2
16	451.586	388	859.2
17	451.614	416	921.1
18	451.643	445	985.3
19	451.673	475	1051.6
20	451.703	505	1118.0
21	451.734	536	1186.5
22	451.765	567	1255.1
23	451.797	599	1325.8
24	451.83	632	1398.8
25	451.864	666	1473.9
26	451.898	700	1549.0
27	451.934	736	1628.6
28	451.969	771	1705.9
29	452.006	808	1787.6
30	452.044	846	1871.5

through the numerical solution of (71), and are plotted in Figure 9. It is obvious based on Figure 9 that increasing the spring constant will increase the natural frequency, and hence there is a frequency shift, which is approximately a parabolic function of the spring constant.

Given an initial height $z$ denoted by $z_j$, the spring constant $k\left(z_j\right)$ and $F_0\left(z_j\right)$ can be evaluated from Equation (7). The solution of the static equilibrium state given by Equation (25) with the left end clamped is easily obtained as

$$w^{*}\left(z_j\right)=\frac{-F_0\left(z_j\right)}{\left(c_{22}{I}_2-\frac{k\left(z_j\right)L_b{}^3}{3}\right)}(\frac{x_2^3}{6}-\frac{x_2^2{L}_b}{2}).$$

(72)

During measurement, the vibration amplitude $\overline{w}$ of the prong is controlled as a small constant denoted by $A_{\mathrm{const}}$ (which is called amplitude modulation AFM, or AM-AFM) so as not to touch the sample, the distance $s\left(z_j\right)$ between the tip and sample is

$$s\left(z_j\right)=z_j-w^{\ast }\left(z_j\right)-A_{\mathrm{const}}$$

(73)

Then, by solving Equation (71) the natural frequency and frequency shift ($\Delta f$) can be obtained for this given initial height $z_j$. Varying the initial height means varying the distance $s\left(z_j\right)$ and the frequency shift can be plotted as function of $s\left(z_j\right)$, as shown in Figure 10 in the far-range regime $20\stackrel{\circ }{\mathrm{A}}\le s\le 60\stackrel{\circ }{\mathrm{A}}$. Including the shot-range repulsive regime $\left(3\stackrel{\circ }{\mathrm{A}}\le s\le 5\stackrel{\circ }{\mathrm{A}}\right)$, $\Delta f$ vs. $s$ is plotted in Figure 11.

As to the experimental work, due to the piezoelectric effect of the single-crystal quartz, QTF can sense the frequency shift. If the tip-sample distance $s\left(z\right)$ can be measured, the frequency shift $\Delta f$ as a function of $s\left(z_j\right)$ can be obtained experimentally as

$$\Delta f=g\left(s\right).$$

(74)

However, the relationship between the frequency shift and the force gradient is not available from QTF measurement. This hindered us finding the tip-sample force, but it can be obtained by solving the equations of motion established from proper physical modeling. From the equation $s=z-w^{\ast }\left(z\right)-\overline{w}$ we have a fixed value of $z$ and therefore the fixed value of static deflection $w^{\ast }\left(z\right)$,

$$\frac{\partial F}{\partial \overline{w}}=\frac{\partial F}{\partial s}\frac{\partial s}{\partial \overline{w}}=-\frac{\partial F}{\partial s}.$$

(75)

The solution of (71) gives

$$\Delta f=h\left(\partial F/\partial \overline{w}\right)=h\left(-\partial F/\partial s\right),\mathrm{or}\frac{\partial F}{\partial s}=-h^{-1}\left(\Delta f\right)$$

(76)

With the use of (74) and (75), Equation (76) becomes

$$\frac{\partial F}{\partial s}=-h^{-1}\left(g\left(s\right)\right)=p\left(s\right)$$

(77)

The integration of Equation (77) provides the tip-sample force as

$$F\left(s\right)\approx F\left(s\right)-F\left(\infty \right)={\displaystyle {\int }_{\infty }^{s}p\left(s\right)ds}$$

(78)

In short, the QTF experiment produces the frequency shift as function of distance $s$, such as in Equation (74). Physical modeling gives the frequency shift as function of the force gradient, like in Equation (76). Combining these two functions and performing integrations numerically, the tip-sample interaction force as a function of $s$ comes out. Since the interaction force is the negative gradient of the potential, therefore, the integration of Equation (78) again gives a negative potential for the tip-sample interaction.

4.4. Experimental Test for the Natural Frequency

To support the proposed physical model based on elasticity and dynamics, an experimental test for the natural frequency of the first flexural-bending and anti-phase mode was executed. If a QTF of tiny size as the one used in Section 4.3 may not be available at hand, a QTF of moderate size can be obtained from a wristwatch by etching out the epoxy packing. Since the geometrical dimensions influence the natural frequency significantly, the length, width, and gap between the two prongs are measured by microscope as shown in Figure 12, Figure 13 and Figure 14. The length and width of the proof mass are shown in Figure 15 and Figure 16. The thickness of the QTF is measured by using SEM as shown in Figure 17.

The results of analytical solution, finite element simulation using commercial code COMSOL, and the experimental measurements are listed in

Table 2. Natural frequency of the first in-plane anti-phase mode from experiment, FEM (finite element method), and theoretical analysis.

	Experiment	Analysis	FEM
Frequency (kHz)	32.768 (kHz)	33.367 (kHz)	31.466 (kHz)
Error (%)	-	1.83%	−3.97%

.

From Table 2, it is clear that the natural frequency predicted by theoretical analysis is closer to the experimental one than that obtained by finite element analysis. The reason for this is that the theoretical result is obtained by solving the equations of motion of the system exactly, so the eigenvalue (natural frequency) and eigenfunction (mode shape) should approach the true solutions. The only error comes from the small discrepancy between our equations of motion and the true equations of motion, while the method of finite element uses the admissible functions to replace the unknown eigenfunctions. The solution process is explained in Appendix C.

5. Conclusions

Although non-contact QTF can provide the frequency shift, due to the tip-sample interaction as a function of their distance, the high resolution images of the tested molecules, atoms, or living cells need precise physical modeling for providing the relationship between the frequency shift and the force gradient. The popular physical model contains a simple harmonic oscillator or two coupled oscillators. However, it encounters the problem of determining the spring constant correctly. This is still controversial among researchers. In the paper, we proposed a new physical model, which models the QTF as elastic bodies instead of point mass oscillator. Based on the elastic theory and analytical dynamics, we can completely solve the problem of finding the relationship between the frequency shift and force gradient.