Nonlinear behaviour of self-excited microcantilevers in viscous fluids

Microcantilevers are increasingly being used to create sensitive sensors for rheology and mass sensing at the micro- and nano-scale. When operating in viscous liquids, the low quality factor of such resonant structures, translating to poor signal-to-noise ratio, is often manipulated by exploiting feedback strategies. However, the presence of feedback introduces poorly-understood dynamical behaviours that may severely degrade the sensor performance and reliability. In this paper, the dynamical behaviour of self-excited microcantilevers vibrating in viscous fluids is characterized experimentally and two complementary modelling approaches are proposed to explain and predict the behaviour of the closed-loop system. In particular, the delay introduced in the feedback loop is shown to cause surprising non-linear phenomena consisting of shifts and sudden-jumps of the oscillation frequency. The proposed dynamical models also suggest strategies for controlling such undesired phenomena.


Introduction
Microcantilever-based sensors have been recently proposed as a promising platform capable of achieving very high sensitivity in rheology and mass-sensing applications on extremely low volumes of fluids. Rheological properties and mass/ concentration of molecules of interest are often measured by shifts in resonance frequency of externally excited cantilevers [1]. In fact, attachment of target molecules to the cantilever or changes in fluid density or viscosity translate to an increase in the probe effective mass and damping that, in turn, induces changes in the cantilever resonance frequency. However, the main drawback of using an external excitation when operating in viscous fluids is the low quality factor (Q) of the resonator and, especially for standard acoustic excitation, the presence of unwanted modes of oscillations due to the fluid-probe interaction. Such unwanted modes make the identification of the resonance peak very challenging and require complicated tuning and setup procedures [2].
Such poor performance is due to the dramatic decrease of the microresonator quality factor when operating in liquid. Several strategies have been proposed to overcome such limitation via the exploitation of feedback loops aimed at making the frequency response more selective. Q-control is one of those strategies used to increase the selectivity of the probe frequency response (thus making the resonance peak more evident), but the circuit responsible for the enhancement of the Q-factor often requires ad-hoc tuning procedures [3].
Similarly, feedback strategies to effectively increase the Q-factor in atomic force microscopy imaging applications have also been proposed. In this case, highly resonant microcantilevers allow for a decrease of the interaction force with the sample [4][5][6]. Strategies as Q-control [7,8], parametric resonance [9,10] and self-excitation circuits [11,12] have allowed these setups to be successfully used for imaging and as rheological [13], mass [14] or chemical [8] sensors.
However, embedding the microcantilever in a feedback loop can induce undesired behaviours, and the nonlinearities that are either introduced in the loop or intrinsic of the canti lever mechanical response can generate poorly-understood phenomena and instabilities. Potential sources of nonlinearities include: nonlinear electronic components required for signal conditioning [11,12], nonlinear stress-strain relation of the cantilever material subject to large deformations [15,16] and nonlinear forces between the probe and the sample or surrounding viscous fluids [17]. As a result, the dynamics of these systems are often very complex and several cases of chaos [18] and bifurcations [19,20] on the response of the resonators have been reported.
Achieving fine control of the dynamics of the resonators can be extremely important to avoid undesired phenomena during device operation and to use the most advantageous features for a specific application. Therefore, it is crucial to develop accurate models to describe the resonator vibrating in closed-loops.
In this paper, self-sustained oscillations of microcantilevers operating in viscous fluids are induced using a feedback loop proposed in [11], in which a controllable phase shifter (PS) has been inserted in order to assess the influence of feedback delay on the cantilever response. This aspect has often been overlooked in the literature, but here it is shown to significantly affect the sensor mechanical response, sensitivity and reliability. The frequency of self-excited oscillation is studied as function of the delay introduced in the circuit by the PS and as a function of the rheological properties of the viscous fluid.
Smooth changes and sudden jumps on the frequency of oscillation are observed. Such nonlinear behaviour is generic and has been tested on three different cantilevers operating in two different viscous fluids (air and water). The observed jumps are expected to occur in any type of device based on this setup and, if used in a controlled way, may pave the way for developing more reliable and sensitive viscosity and mass sensors, or even provide new potentialities to different scanning probe techniques, such as acoustic force microscopy [21]. Two different models are proposed to explain the experimental observations and to predict the response of the self-excited microcantilever.
The paper is organized as follows: in section 2 the experimental setup is discussed. Special emphasis is given to the description and characterization of the PS, which is the critical component exploited here to highlight the nonlinear behaviours. Experimental results illustrating the nonlinear behaviour of the self-excited oscillation frequency are presented in section 3. In section 4 two mathematical models to describe the system are developed and discussed. Predictions provided by such models are compared to experimental results in section 5. Finally, some conclusions and outlooks for future opportunities are discussed in section 6.

Experimental setup
A schematic of the experimental setup used to characterise the behaviour of microcantilevers oscillating in different viscous media is presented in figure 1(A). The cantilever is acoustically excited using a dither piezo and the resulting deflection is measured using a four-quadrant detector connected to a R9 controller (RHK Technology).
The switch S in the diagram allows users to select between a classical externally excited configuration (amplitude modulation, AM) and a self-excitation setup (auto-tapping mode, AT). In the former modality the dither piezo is externally excited by a function generator implemented in the controller and the amplitude and phase of the deflection signal are read via a lock-in available in the controller itself. In AT mode, the measured deflection is fed to an electronic circuit (Elbatech srl) composed of an adjustable gain K, a saturator and a PS, whose output is connected to the dither piezo voltage input. The frequency of vibration of the cantilever is measured by feeding the deflection signal to a spectrum analyser embedded in the R9 controller (thermal spectrum function).
The experimental setup of figure 1(A) explicitly indicates the total delay, τ tot , of the signal along the feedback loop. As it will be discussed below, this total delay encloses three The cantilever is immersed in a viscous medium and excited by a dither piezo. In the amplitude mode (AM), the external excitation signal drives the dither piezo at different frequencies. In the autotapping mode (AT), the deflection signal passes through the feedback loop, where it is delayed by the phase shifter (PS), amplified by the gain and limited by the saturator, being finally fed back to the piezo as a voltage. (B) Detail of a single stage of the PS. Two stages, connected in series, were used, each one capable of shifting the signal by at most 180°. Values of C 1 and C 2 are 4.7 nF and 220 pF in each stage, respectively, and R 1 and R 2 are adjustable from 0-10.4 kΩ. The polarity of the voltage applied to the dither piezo can also be inverted.
individual contributions from elements of the loop and creates a natural phase shift between the excitation of the piezo and the motion of the tip of the cantilever.
In addition to the natural phase shift induced by the delay τ tot , a PS was inserted in the feedback loop to finely control the phase between dither piezo excitation and cantilever deflection. The complete PS consists of two stages connected in series, each one working as an all-pass filter able to shift the phase of the signal by at most 180°. Figure 1(B) shows the detail of the electrical scheme of a single stage. The two stages are individually operated by adjusting two potentiometers which control the value of two resistors, R 1 and R 2 , between 0-10.4 kΩ. The phase shift ϕ = −2atan (ωR i C i ) introduced by each stage of the PS depends on the frequency of oscillation of the loop, ω, and the values of the resistor and capacitor, R i , C i . In addition, the polarity of the voltage applied to the terminals of the dither piezo that excites the cantilever can also be inverted, allowing an extra phase shift of 180° on the signal. The two stages of the PS and the inversion of polarity in the piezo can be used to shift the signal along the feedback loop by a complete cycle (360°). This strategy will be used to assess the influence of feedback delay on the cantilever response.
A full characterization of each element in the experimental setup is described below.

Dynamics of acoustically excited microcantilevers immersed in viscous fluids in AM mode
Three different commercial cantilevers were used in this work, covering a wide range of resonance frequencies and stiffness. Their specifications and geometries are presented in table 1. In addition, two distinct viscous fluids were used to study the vibration of each cantilever: air and water.
Each cantilever was characterized at first by performing frequency sweeps in AM mode and measuring its frequency response in the viscous fluids. The recorded frequency response in AM mode, described by both amplitude and phase spectra, is discussed below, and will be used in section 4.1 to model the AT mode.
2.1.1. Amplitude spectra. One example of amplitude spectra of the cantilever CLFC-B immersed in air is presented in the top panel of figure 2. The inset shows the presence of two distinct resonance modes. The amplitude, A, of the measured amplitude spectrum can be fitted via a simple harmonic oscillator (SHO) model for each resonance peak [22]: (1) where ω is the excitation angular frequency and ω 01 , ω 02 , Q 1 and Q 2 are the natural angular frequencies and quality factors of the first and second modes, respectively. A 1 and A 2 are the amplitude of each individual resonance mode and depend on the applied external force and the effective mass of each mode. The total amplitude, A, is the sum of the amplitudes of each mode. For excitation frequencies ω far from ω 02 , the second term of equation (1) is negligible and the behaviour of the first mode is captured (and vice-versa).
The parameters ω 01 , ω 02 , Q 1 , Q 2 , A 1 and A 2 were used to fit the model of equation (1) to the experimental amplitude spectra of the three different cantilevers immersed in different viscous fluids. The results of the fits for each cantilever and medium are shown in table 2. In some experimental conditions the second mode is not observed (A 2 = 0) as it is either too damped or falls outside the frequency range of the lock-in used for this experiment. Moreover, the presence of multiple peaks in the cantilever frequency response, resulting from the coupling between the acoustic excitation of the piezo and spurious mechanical modes [23,24], may sometime make it difficult to isolate the real resonance peak. In these cases, equation (1) was fitted by starting from the estimates given by the thermal noise spectrum and then 'enveloping' all the forest of peaks measured in AM mode under one single SHO peak. The problem of the forest of peaks is thought to have its origin in resonances of the liquid cell or even of the dither piezo shaker. Strategies such as upgrading the setup with acoustic barriers [25] or optical isolation [26] to suppress the noise were proposed. In alternative, magnetic or photothermal excitation techniques [27,28] can also be used to avoid this problem.

Phase spectra.
The total delay block, τ tot , shown in figure 1(A), encloses three contributions of elements of the feedback loop. The first contribution is given by the dither piezo and cantilever, τ CT , and it is related to the time the acoustic waves take to propagate from the dither piezo to the cantilever through the material composing the cantilever holder. Such propagation translates to a delay between the excitation provided by the dither piezo and the deflection of the cantilever.
This acoustic delay can be estimated from the slope of the phase spectrum of the cantilever, away from the regions where the resonance jump occurs. In fact, far from the resonance peak, the phase shift due to the cantilever intrinsic resonance can be safely neglected and the phase shift (in radians) between dither piezo excitation and cantilever deflection simply reads where ω is the angular frequency of the dither piezo excitation. Note that in AM mode the loop is open, and that the measured delay does not contain information about the other components of the circuit. Figure 2(B) shows the AM mode phase spectrum for the cantilever CLFC-B in air, from where a delay τ CT = 6.6 µs is extracted. The expected jump of 180° around the resonance is also observed. Table 2 shows the delays associated with each individual cantilever, extracted using this method in the phase spectra measured in air. The values of the measured delays of the cantilevers are near one order of magnitude higher than those expected considering simply the speed of the acoustic waves in the silicon beam. This fact indicates that the delay introduced by the cantilever in the loop is mainly due to the connection between the dither piezo and the cantilever.

AT electronics
A dedicated electronic circuit has been developed by Elbatech srl based on the basic logic explained in [11] and augmented with the adjustable PS described in figure 1(B), which enables the phase of the AT feedback loop to be changed in a controlled way.
The PS is composed of two all-pass filters implemented using the operational amplifier circuit reported in figure 1(B). To characterize its frequency response, the output of the signal generator of the R9 controller was connected to the PS input, and the shifted output was connected to a lock-in input; this allows the user to measure amplitude and phase response of the circuit over a wide range of frequency. Figure 3 reports the measured phase shifts for different values of the resistors R 1 and R 2 . The experimental results presented in figure 3 (solid lines) are modelled by the transfer function where H i ( jω) = 1−jωRiCi 1+jωRiCi represents the transfer function of each individual stage of the circuit, and the parameter p is used to distinguish between non-inverted or inverted polarity on the terminals of the dither piezo (the convention p = 1 for Table 2. Parameters of the SHO model fitted to the amplitude spectra of the cantilevers vibrating in air and water in AM mode. The delay of each cantilever, τ CT , was experimentally estimated from the phase spectra of the cantilevers operating in air in AM mode, and used as the fitting parameter of the nyquist and simulink models to the experimental results.
AM mode amplitude SHO fit AM mode phase Nyquist model Simulink model  non-inverted polarity and p = −1 for inverted polarity will be followed). It was observed that the delay τ PS is required to accurately model the experimental results. The presence of such delay is due to the propagation of the signals through the electronic components and it represents the second contribution for the total delay, τ tot , shown in figure 1(A). The dotted lines of figure 3 show best fits of equation (3) to the experimental results with non-inverted polarity ( p = 1), resulting in a delay τ PS = 1.0 µs.
The saturation circuit is implemented using operational amplifiers, as described in [12,29]. Its response can be described by the function: where a represents the signal through the loop and σ is the saturation threshold, defined by the user. The solid wine line in figure 3 shows the measured phase delay introduced by the gain and saturator elements as function of the frequency of the input signal. Ideally such elements should not introduce any delay in the feedback loop, but experimental data suggest that a pure delay τ ET = 1.1 µs is required to correctly model the system. This delay is the third contribution for the total delay, τ tot , shown in figure 1(A).

Experimental results in AT mode
Experiments were performed by immersing the cantilevers in the two different viscous fluids and measuring the frequency of the self-excited cantilever oscillation. Such frequency was then studied as function of the delay introduced by the PS for different values of the potentiometers R 1 or R 2 . The frequency of oscillation was measured by feeding the deflection signal of the cantilever to a spectrum analyser and acquiring the FFT spectra (thermal mode available in the R9 controller).
When the values of R 2 were gradually increased (keeping a fixed polarization on the dither piezo and R 1 value) the oscillation frequency of the loop was observed to change progressively. A sudden jump on the oscillation frequency was then observed for very small variations of R 2 , and then the smooth changes on the frequency of oscillation resumed. The sudden jumps on the frequency of oscillation were observed either within the first resonance mode or between two adjacent resonance modes.
This behaviour is illustrated in figure 4, which shows a jump within the first mode of the ACST cantilever operating in air, when the values of R 2 are swept up, using non-inverted polarization on the piezo (p = 1) and R 1 = 0.01 kΩ. The frequency of oscillation decreases progressively with the increase of R 2 (from violet to green solid curves), before jumping to a much higher value (from green to black curves) with a very small change in R 2 value (from 1.08-1.10 kΩ). Finally, the frequency of oscillation keeps decreasing with the increase of R 2 (from black to red curves). It can also be observed that there is a component of motion of vibration at high frequencies that becomes more visible as the jump approaches. The inset of figure 4 shows the self-sustained oscillations of the loop as function of the values the potentiometer R 2 .
The measured frequencies of oscillation enclose the natural fundamental frequency of the cantilever ACST in air, as shown in table 2 (155.8 kHz). This general behaviour was observed for all the studied cantilevers and viscous fluids. The sudden jumps were observed to occur more often within the first resonance mode than between adjacent resonance modes, as the amplitude of the second mode is usually much smaller than the amplitude of first mode.
Such results clearly highlight the need of considering the phase shift when modelling any feedback loop used to excite the cantilever and the possibility to exploiting it to modify the cantilever mechanical response.

Mathematical modelling
In order to explain the nonlinear behaviour observed in the experiments described in section 3, here two distinct mathematical models are proposed. The first approach is based on a graphical method that exploits the Nyquist stability criterion-a classical result in feedback control theory-to determine the frequencies of oscillation of the circuit and to predict the sudden jumps. This strategy requires the use of the fitted parameters of the SHO model applied to the AM mode spectrum of each cantilever immersed in the viscous fluids (frequencies, quality factors and amplitudes shown in table 2). Although such model is capable of providing good predictions, it requires users to fit several parameters on experimental data and such data may not be always readily available.
The second approach approximates the cantilever behaviour with a single-degree-of-freedom equation of motion of a nonlinear resonator that can be easily simulated numerically. In this approach, the hydrodynamic force that acts on the vibrating cantilever is modelled using Sader's hydrodynamic function [30] and of the viscosity and density of each viscous medium are incorporated into the mass and damping parameters of the resonator. Therefore, it is possible to solve the equation of motion without any need of fitting experimental values.

Nyquist stability criterion
The schematic of the experimental setup shown in figure 1(A) is an example of Lure's system, where linear systems (cantilever dynamics, gain, delay and PS) are feedback connected via a nonlinear function (saturator). In such a system, the onset of stable self-sustained oscillations results from a competition between the feedback gain, which constantly amplifies the motion of the cantilever (inducing instability), and the presence of the nonlinear saturation, which limits the system trajectories and stabilizes the system dynamics on stable selfsustained oscillations. Given that the resonant response of the microcantilever effectively acts as a band-pass filter, the harmonic balance technique can be successfully exploited to better understand and predict the system behaviour (see [11] for an example of successful application of such technique for imaging applications). This technique assumes that every periodic output signal of a nonlinear block (subjected to a sinusoidal input signal) can be approximated by the first terms of associated Fourier series, i.e. the nonlinearity output is approximated by a sinusoidal wave having the same frequency as the input. This approximation is reasonable since the resonator possess intrinsic band-pass filter characteristics that attenuate the low frequencies and higher harmonics [31]. Within this framework, the saturator described in section 2.2 (equation (4)) can be replaced by an amplitude-dependent 'gain', called describing function, (5) Equation (5) shows that if the amplitude a of the input signal is smaller than the threshold σ defined by the user, the gain is unitary and the output signal is the same as the input signal. If the amplitude of the signal is higher than the threshold value, the output signal is decreased with respect to the input signal, which contributes to stabilizing the oscillation of the cantilever.
According to such approximation, a condition for existence of self-sustained oscillations on the feedback loop shown in figure 1(A) reads [31] where e −jωτtot is the transfer function of the total delay (τ tot = τ CT + τ PS + τ ET ), K represents the self-excitation loop gain, PS ( jω) is the transfer function of the PS (the same of equation (3) but omitting the delay τ PS , which is already contained in τ tot ), CT( jω) is the transfer function of the cantilever (equation (1)) and p = ±1, depending on the polarity applied on the terminals of the dither piezo. Condition (6) is equivalent to stating that the gain of the loop must be unitary and can be simplified to with G ( jω) = CT ( jω) e −jωτtot PS ( jω). Since the gain K and the describing function ψ(a) are real functions for each value of amplitude a and oscillation frequency, equation (7) can be re-written as [31]: which, by its turn, can be decomposed into two real equations [31,32]: The condition imposed by equation (9b) states that the total phase shift around the loop must be an integer multiple of 2π radians. The oscillation frequency is determined by this condition. Using the oscillation frequency, equation (9a) can then be solved to calculate the amplitude of vibration (this step is not needed for the results discussed in this paper). When the circuit is initially turned on, the frequency of the white noise that satisfies the condition of equation (9b) is initially amplified by the gain. Later, the gain and saturator will compete to stabilize an oscillation with an amplitude that also satisfies equation (9a). In the absence of the saturator, the trajectories of the resonator would be naturally limited by the force that the dither piezo can exert on the cantilever base or by the intrinsic mechanical nonlinearities of the cantilever [11,12].
Given the highly nonlinear nature of equation (9), they need to be solved numerically to predict the oscillation frequency and amplitude of vibration. In this case, equation (9b) is re-written by substituting the individual transfer functions of the cantilever, total delay and PS in G ( jω): Equation (10) is then solved using Matlab ® to find the frequency of self-sustained oscillation. The parameters appearing in this equation can be estimated as described in section 2 by the spectra obtained in AM mode. The delay associated with the piezo and cantilever, τ CT , is used as the fitting parameter of the model to the experimental data and the final results are shown in table 2. Nevertheless, some remarks on solving equation (10) graphically must still be done: equation (10) states that the frequency of self-oscillation will be such that the imaginary part of G ( jω) equals zero, i.e. the Nyquist diagram intersects the horizontal axis, but does not point out explicitly which crossing the frequency of oscillation corresponds to. However, the Nyquist stability criterion states that only the frequency ω auto corresponding to Im [G ( jω auto )] = 0 and Re [G ( jω auto )] > Re [G ( jω)] (ω > 0) corresponds to stable oscillations [31]. By comparing the predictions given by this model with the experimental data, it was also established that the solutions for non-inverted polarity on the piezo (p = 1) correspond to the positive outer crossing of equation (10), whereas solutions for inverted polarity (p = −1) correspond to the outer crossing with the negative horizontal axis. Figures 5 and 6 will be used to illustrate how this model is able to explain the observed phenomena. Figure 5 presents the Nyquist diagram of G ( jω) (equation (10)) for the case of the ACST cantilever vibrating in water (used parameters are shown in table 2). In this experiment, R 2 was swept up, as in the case discussed in figure 4, using non-inverted polarity (p = 1) and a fixed R 1 = 2.25 kΩ. Figure 5(A) shows that an increase in R 2 causes a clockwise rotation of the Nyquist diagram of G ( jω), as it increases the overall phase shift in G ( jω). These curves are parameterised by frequency and, in this specific case, only consider the first mode (A 2 = 0 in equation (10)). Figure 5(B) shows that the frequency associated with the positive outer crossing of the Nyquist diagram of G ( jω) with the zero of the Im axis changes as the value of the potentiometer R 2 increases. Furthermore, at a certain value of R 2 , the crossing moves from a low-frequency region of the G ( jω) curve (purple line) to a high-frequency region of the curve (green line). These two mechanisms explain the shifts and jumps on the oscillation frequency of the loop inside the first mode, as was experimentally observed and plotted in figure 4.
Finally, according to this model, there exists a specific value of R 2 that will cause both the low and high-frequency parts of the curve G ( jω) to cross simultaneously the zero of the Im axis. In this case, the oscillator will work in a bistable state [33]. This was also observed in figure 4, where, as the jump approached, the presence of the high-frequency motion started to be visible (green line). Figure 6 illustrates one of the less common cases in which the sudden jumps on the oscillation frequency of the loop occur between consecutive resonance modes, as it was observed for the cantilever CLFC-B operating in air, using non-inverted polarity (p = 1) and a fixed R 1 = 0.82 kΩ. In these cases, the transfer function G ( jω) has to contain information about the two resonance modes (A 2 ≠ 0 in equation (10)) and the resulting Nyquist diagram is slightly more complex: each resonance mode corresponds to a big circle in figure 6(A). Furthermore, the entire region of low-amplitude (the plateau between peaks shown in the inset of figure 2(A)) is barely visible in figure 6(A), since the values are very small and close to the origin. Figure 6(A) shows that, as before, an increase in the values of R 2 cause a clockwise rotation on the Nyquist diagram of G ( jω). Since these curves are parameterised by frequency, the frequency associated with the positive outer crossing of the Nyquist diagram with the zero of the Im axis changes as the values of the potentiometer R 2 increase, which explains shifts in the frequency of oscillation of the loop. In addition, as seen in figure 6(B), at a certain value of R 2 the outer crossing moves from the big circle on the curve of G ( jω), relative to the first mode (purple line), to the smaller circle, relative to the second mode (green line). This explains the observed jumps between consecutive resonance modes.
Finally, it is worth pointing out that for the jumps between different modes to occur, the second mode of the cantilever immersed in a specific viscous fluid must have a nonnegligible amplitude when compared with the amplitude of the first mode.

Dynamical numerical model: hydrodynamic force
The approach described in section 4.1 is able to correctly predict the cantilever behaviour but it requires several parameters to be estimated from experimental data. To avoid such fitting procedure, here the cantilever is approximated via a nonlinear single-degree-of-freedom resonator working in a feedback loop. The total force exerted by the fluid on the vibrating continuous beam depends on the viscosity and density of each viscous medium and is calculated using Sader's hydrodynamic function [30]. This total hydrodynamic force is thought as a combination of normal and tangential terms, correspondent to the pressure and viscous forces that act on every surface of the immersed beam. The first is an inertial term, described by the weight of the layer of fluid that the beam displaces as it moves. The second term is proportional to the velocity of the beam and accounts for the viscous drag force exerted by the fluid on the moving cantilever. These two components can therefore be described as an added mass, m A , and an added damping coefficient, c A [34,35]: where Γ = a 1 + a 2 δ(ω) W are expressions to approximate the hydrodynamic function of a rectangular cantilever [34], which depend on the constants, a 1 = 1.0553, a 2 = 3.7997 and b 1 = 3.8018 and b 2 = 2.7364 [34]. The parameter δ (ω) = 2η ρω is the thickness of the layer of fluid that surrounds the beam, which depends on the fluid viscosity, η, the fluid density, ρ, and on the angular frequency of oscillation, ω. L and W represent the length and width of the cantilever, respectively [35].
The single-degree-of-freedom equation of motion that describes the vibration of the cantilever immersed in a viscous fluid in the feedback loop is therefore given by [33] (m ct + m A )ÿ i (t) + (c i + c A )ẏ i (t) + k i y i (t) = F(t), (13) where the dot stands for time derivative and F(t) is the forcing term from the dither piezo. m CT = LWTρ c is the total mass of the beam, where T and ρ c are, respectively, the thickness of the beam and the density of the constituent material, while 2 m CT and y i are the intrinsic damping coefficient, stiffness and displacement associated to the i-th resonance mode, respectively. In the previous terms, f 0i is the natural frequency of each mode, which can be experimentally measured or estimated by analytical equations [36], while Q i is the quality factor of each mode. The total displacement of the resonator is the sum of the displacements of each mode, i.e. y(t) = y 1 (t) + y 2 (t), considering only the first two modes.
The force from the dither piezo that acts on the microcantilever in equation (13) is described by where τ tot is the self-excitation loop total delay, K is the feedback gain and y PS is the output of the PS. The function sat() describes the saturator detailed in section 2.2 by equation (4), with user-defined threshold value σ. A block diagram of the dynamical model is shown in figure 7. The red dashed arrows show the parameters that have to be used as inputs of the model. The value of the delay τ CT introduced by the piezo and cantilever is used as the fitting parameter of the model. Figure 7 shows that the mass of the cantilever, m CT , is first used to calculate the spring constant, k i , and intrinsic damping coefficient, c i , of each mode, using experimental or estimated values for the resonance frequencies, f 0i , and quality factors, Q i , respectively. Then the constants a 1 , a 2 , b 1 , b 2 , the geometry of the cantilever and the properties of the viscous fluid (η and ρ) are used to calculate the added mass and damping coefficient of the system, through equations (11) and (12).
The function of the saturator, described by equation (4) in section 2.2, is implemented with user-defined threshold value, σ. Similarly, the transfer function of the PS, described by equation (3) in section 2.2 (but omitting the delay τ PS , Figure 7. Block diagram of the dynamical model. All the parameters that are used as inputs are highlighted with a dashed red arrow. The hydrodynamic force that acts on the vibrating beam is calculated from the properties of the fluid and the cantilever geometry and each resonance mode is defined by its own spring constant. The delays associated with the PS and electronics were estimated from experimental data, and the delay introduced by the piezo and cantilever was used as the fitting parameter of the model. which is already contained in τ total ), is implemented with userdefined values of the potentiometers R 1 and R 2 , and polarity on the terminals of the piezo (p = 1 or p = −1). The values of C 1 and C 2 are 4.7 nF and 220 pF, respectively. The values of the delays associated to the PS, τ PS = 1.0 µs, and electronic components gain plus saturator, τ ET = 1.1 µs, were estimated as explained in section 2 and are fixed in the model, in the parameter τ total . Finally, the value of delay introduced by the piezo and cantilever, τ CT , is used as the fitting parameter of the model to the experimental results. The results are shown in table 2.
Equation (13) is then solved in time-domain using Simulink until the steady-state is reached, and the frequency of selfoscillation, f osc , is extracted.

Comparison between model predictions and experimental results
To highlight the crucial role played by the feedback phase shift, experiments were conducted by immersing the cantilevers in the two viscous fluids (air and water), and measuring the frequency of self-oscillation as a function of the values of R 1 , R 2 and the polarity of the dither piezo. The typical proto col was sweeping up the value of one of the potentiometers, while keeping the other parameters constant. The polarity of the piezo was then changed and the potentiometers sweep repeated. In this section, all the experimental results are compared with the results modelled by the two strategies discussed in section 4. The experimental and simulated frequencies of oscillation were plotted against the phase shift introduced by the PS: Equation (15) is a direct calculation of the phase of the PS transfer function given by equation (3), where ω is the measured or simulated angular frequency of oscillation in the loop, C 1 = 4.7 nF and C 2 = 220 pF are capacitances of the capacitors used in the architecture of each stage of the PS, R 1 and R 2 are the used values of the potentiometers in the experimental sweeps and in the models, τ PS is the delay associated with the PS (τ PS = 1.0 µs estimated in section 2.2) and the parameter P is used to describe the phase changes due to the inversion of polarity on the dither piezo (P = 1 if p = 1 or P = 0 if p = −1). Figure 8 shows the complete set of results for the cantilever CLFC-B immersed in air and water. The black circles represent the experimental measurements and the red triangles and blue diamonds represent the results obtained with the Nyquist approach (section 4.1) and the dynamical model (section 4.2), respectively. Figure 8(A) presents a case where the jump occurs between the first and second resonance modes, while figure 8(B) shows a jump occurring within the second resonance mode. The values of delay of the cantilever, τ CT , were used as the fitting parameter of the models to match the precise location of the jumps. Figure 9 presents additional results obtained with different cantilevers working in the viscous media. Figures 9(A) and (B) show jumps on the oscillation frequency within the first resonance mode.
As discussed, the condition for the existence of selfsustained oscillations is that the total phase shift around the feedback loop must be an integer multiple of 2π radians. The observed shift of the oscillation frequency corresponds to the cantilever adjusting its phase (and hence its oscillation frequency) so to compensate the delay imposed by the other components of the feedback loop. Therefore, it is expected that the observed nonlinear phenomena repeat themselves every 2π radians, being fairly independent of the delay introduced by the electronic components of the loop, τ ET , and by delay introduced by the piezo and cantilever, τ CT . Table 2 shows that the values used to fit the models in each viscous fluid are closely related to the values measured using the phase of AM mode spectra of each cantilever in air. In all conditions, the proposed models are able to correctly describe the nonlinear behaviour of the microcantilever within 15% error and can be used to predict the behaviour of the cantilevers in viscous fluids working in AT configuration.

Conclusions
In this work the non-linear behaviour of microcantilevers oscillating in viscous fluids and excited by a non-linear feedback loop is studied. The feedback loop comprises a linear gain, a non-linear saturator and a PS. The delay introduced by the PS in the loop can be controlled by two potentiometers and it is shown to play a crucial role in controlling the response of the self-excitation loop. Non-linear phenomena, such as sudden jumps in the oscillation frequency, are observed and modelled by two different strategies. The first strategy requires parameters obtained from the open-loop frequency response of the cantilevers vibrating in the viscous fluid and is based on feedback theory and the Nyquist stability criterion. The second strategy consists on numerically solving the equation of motion of the resonator in closed loop, exploiting the knowledge of the hydrodynamic function and without requiring prior experimental data fitting. Both models capture accurately the non-linear behaviour of the closed-loop and can be used to predict the response of this system when used as a mass or rheological sensor.
The experimental results and the models described in this paper highlight the importance of the delay and phase shifts that are intrinsically present whenever a feedback scheme is used to manipulate the probe dynamical response, but that are very often neglected in the literature. It is expected that similar behaviours will be present for other feedback approaches such as Q-control and parametric resonance.
Moreover, the understanding of the observed non-linear behaviour opens exciting opportunities for the development of a novel class of rheology, chemical and mass sensors using self-excited microcantilevers. In fact, the self-excitation mechanism analysed here allows for very high signal-to-noise ratio in frequency measurements and automatically tracks the oscillation frequency without requiring any external frequency sweeps and the associated forest of peaks. Furthermore, the location of the frequency jump can be easily controlled via the PS, thus opening the possibility of placing it at a particular fluid viscosity/added mass to realize an extremely sensitive threshold detector for such quantities. On the other hand, by placing such jump far away from the region of viscosity/mass of interests, the proposed microsensor output is smooth, as requested by most of the applications. wish to thank the anonymous reviewers for their constructive comments.