Najib Kacem
email: najib.kacem@cea.fr

Sébastien Baguet
email: sebastien.baguet@insa-lyon.fr

Sébastien Hentz
email: sebastien.hentz@cea.fr

Régis Dufour
email: regis.dufour@insa-lyon.fr

Nonlinear phenomena in nanomechanical resonators: mechanical behaviors and physical limitations

Keywords:

In order to overcome the loss of performances issue when scaling resonant sensors down to NEMS, it proves extremely useful to study the behavior of resonators up to large displacements and hence high nonlinearities. A comprehensive nonlinear multiphysics model based on the Euler-Bernoulli equation which includes both mechanical and electrostatic nonlinearities in the case of a capacitive doubly clamped beam is presented. This purely analytical model captures all the nonlinear phenomena present in NEMS resonators electrostatically actuated including bistability, multistability which can lead to several physical limitations such as noise mixing, frequency stability deterioration as well as dynamic pull-in. Moreover, close-form expressions of the critical amplitudes and pull-in domain initiation amplitude are provided which can potentially serve for NEMS designers as quick design rules.

Resumé

Afin de compenser la détérioration des performances des capteurs résonants lors de la réduction de leurs tailles au niveau NEMS, il s'avère extrêmement utile d'étudier le comportement des résonateurs en grands déplacements et par conséquent sous fortes nonlinéarités. Le modèle de poutre présenté est non-linéaire et multiphysique. Ce capteur capacitif est basé sur l'équation d'Euler-Bernoulli et contient des non-linéarités mécaniques et électrostatiques. Ce modèle, purement analytique, capture tous les phénomènes non-linéaires principaux présents dans des nanorésonateurs capacitifs comme la bistabilité et la multistabilité qui peuvent mener à plusieurs limitations physiques telles que le "noise-mixing", la détérioration de la stabilité en fréquence aussi bien que le "pull-in". De plus, des expressions analytiques des amplitudes critiques et de l'amplitude d'initiation du domaine "pull-in" sont obtenues, ce qui peut potentiellement servir aux concepteurs comme guides de pré-dimensionnement analytique de NEMS.

Mot s clés: résonateur NEMS \ Dynamique non-linéaire \ Amplitude critique \ Pull-in \ Comportement mixte

Introduction

Micro and nanoelectromechanical (MEMS/NEMS) devices have been the subject of extensive research for a number of years and have generated much excitement as their use in commercial applications has increased. Indeed, MEMS technology has opened up a wide variety of potential applications not only in the inertial measurement sector, but also spanning areas such as communications (filters, relays, oscillators, LC passives, optical switches), biomedicine (point-of-care medical instrumentation, microarrays for DNA detection and high throughput screening of drug targets, immunoassays, in-vitro characterization of molecular interactions), computer peripherals (memory, new I/O interfaces, read-write heads for magnetic disks) and other miscellaneous areas such as in projection displays, gas detection and mass flow detection.

NEMS are the natural successor to MEMS as the size of the devices is scaled down to the submicron domain. This transition is well adapted with the resonant sensing technique for a large panel of applications. One reason of scaling down resonant sensors to the NEMS size is the ability to detect very small physical quantities by increasing their sensitivity [START_REF] Ekinci | Nanoelectromechanical systems[END_REF]. In particular, NEMS have been proposed for use in ultrasensitive mass detection [START_REF] Huang Ekinci | Ultrasensitive nanoelectromechanical mass detection[END_REF][START_REF] Jensen | An atomic-resolution nanomechanical mass sensor[END_REF], radio frequency (RF) signal processing [START_REF] Nguyen | Micromechanical components for miniaturized low-power communications[END_REF][START_REF] Nguyen | Tunable, switchable, high-q vhf microelectromechanical bandpass filters[END_REF], and as a model system for exploring quantum phenomena in macroscopic systems [START_REF] Cho | Physics -researchers race to put the quantum into mechanics[END_REF][START_REF] Lahaye | Approaching the quantum limit of a nanomechanical resonator[END_REF].

Unfortunately, the nonlinear regime for nanomechanical resonators is easily reachable [START_REF] Kacem | Nonlinear dynamics of nanomechanical beam resonators: improving the performance of nems-based sensors[END_REF][START_REF] Kacem | Dynamic range enhancement of nonlinear nanomechanical resonant cantilevers for high sensitive NEMS gas/mass sensors applications[END_REF][START_REF] Sumpter | The onset of instability in nanostructures: The role of nonlinear resonance[END_REF], so that the useful linear dynamic range of the smallest NEMS devices is severely limited. In fact, many applications we are hoping for in the near future will involve operation in the nonlinear regime, where the response to the stimulus is suppressed and frequency is pulled away from the original resonant frequency.

Actually, it is a challenge to achieve large-amplitude motion of NEMS resonators without deteriorating their frequency stability [START_REF] Feng | Phase noise and frequency stability of very-high frequency silicon nanowire nanomechanical resonators[END_REF]. The relative frequency noise spectral density [START_REF] Robins | Phase Noise in Signal Sources[END_REF] of a NEMS resonator is given by:

S f = 1 2Q 2 S x P 0 ( 1
)
where Q is the quality factor, S x is the displacement spectral density and P 0 is the displacement carrier power, ie the RMS drive amplitude of the resonator 1 2 A 2 . Remarkably, driving the resonator at large oscillation amplitude leads to better signal to noise ratio (SNR) and, thus, simplifies the design of the electronic feedback loop. However, doing so in the nonlinear regime reduces the sensor performances since the frequency instability of a nonlinear resonator is proportional to its oscillation amplitude. Moreover, even when NEMS resonators are used as oscillators in closed-loop, a large part of noise mixing [START_REF] Howe Roessig | Nonlinear mixing in surface-micromachined tuning fork oscillators[END_REF][START_REF] Kaajakari | Phase noise in capacitively coupled micromechanical oscillators[END_REF] due to nonlinearities drastically reduces their dynamic range and alters their detection limit. Indeed, the combination of the mechanical and electrostatic nonlinearities causes a direct voltage-to-frequency effect at the electrodes, leading to a mixing of the low-frequency 1 f close to the carrier power signal [START_REF] Howe Roessig | Nonlinear mixing in surface-micromachined tuning fork oscillators[END_REF]. In this paper, the main ideas of [START_REF] Kacem | Nonlinear dynamics of nanomechanical beam resonators: improving the performance of nems-based sensors[END_REF] are highlighted in order to investigate in details some nonlinear phenomena which limit the performances of resonant sensors. Based on a compact and analytical model including all main sources of nonlinearities of a two-port resonator, the bifurcation points are used to provide close-form solutions of the critical amplitude and dynamic pull in voltage.

Model of a nonlinear two-port resonator

A clamped-clamped beam is considered (Figure 1) subject to a viscous damping with coefficient c per unit length which models all dissipation sources including ohmic losses [START_REF] Seoánez | Dissipation in graphene and nanotube resonators[END_REF][START_REF] Sazonova | A tunable carbon nanotube resonator[END_REF]. The device is actuated by an electric load v(t) = -V dc + V ac cos( Ωt ), where V dc is the DC polarization voltage, V ac the amplitude of the applied AC voltage, and Ω the excitation frequency.

2

Produit par HAL -19 Sep 2014

Equation of motion

Let E be the Young's modulus, I, the geometrical moment of inertia of the cross section, ρ, the mass density, l, the length, b, the width, h, the thickness of the resonator. Its out-of-plane deflection w(x, t) which depends on c, the viscous damping, Ñ the axial load, and on the forces provided by the two electrodes is governed by the following nonlinear Euler-Bernoulli equation which is the commonly used approximate equation of motion for a thin beam [START_REF] Landau | Theory of Elasticity[END_REF].

ρbh ∂ 2 w( x, t) ∂ t2 + c ∂ w( x, t) ∂ t + EI ∂ 4 w( x, t) ∂ x4 -        Ñ + Ebh 2l l 0 ∂ w( x, t) ∂ x 2 d x       ∂ 2 w( x, t) ∂ x2 = 1 2 ε 0 bC n1 V ac cos( Ωt ) -V dc 2 (g a -w( x, t)) 2 H 1 ( x) - 1 2 ε 0 bC n2 [V s -V dc ] 2 (g d + w( x, t)) 2 H 2 ( x) (2)
H 1 ( x) = H x - l + l a 2 -H x - l -l a 2 (3)
H 2 ( x) = H x - l + l d 2 -H x - l -l d 2 (4)
where ε 0 is the dielectric constant of the gap medium. The two terms of the right-hand member in Equation (2) represent an approximation of the electric forces assuming a resonator design with two stationary electrodes : electrode 1 for the actuation and electrode 2 for the sensing including the fringing field effect [START_REF] Nishiyama | Capacitance of a strip capacitor[END_REF] using the coefficients C ni . The boundary conditions are:

w(0, t) = w(l, t) = ∂ w ∂ x (0, t) = ∂ w ∂ x (l, t) = 0 (5)

Nondimensionalisation

For convenience, let the following nondimensional variables be introduced:

w = w g d , x = x l , t = t τ (6)
where τ = 2l 2 h 3ρ E . Substituting Equation (6) into Equations ( 2), (3), ( 4) and ( 5), gives:

∂ 2 w ∂t 2 + c ∂w ∂t + ∂ 4 w ∂x 4 -        N + α 1 1 0 ∂w ∂x 2 dx        ∂ 2 w ∂x 2 + α 2 C n2 [V s -V dc ] 2 (1 + w) 2 H 2 (x) = α 2 C n1 [V ac cos(Ωt) -V dc ] 2 (R g -w) 2 H 1 (x) (7) w(0, t) = w(1, t) = ∂w ∂x (0, t) = ∂w ∂x (1, t) = 0 (8)
The parameters appearing in Equation ( 7) are:

H 2 (x) = H(x - l + l a 2l ) -H(x - l -l a 2l ) (9)
H 2 (x) = H(x - l + l d 2l ) -H(x - l -l d 2l ) ( 10
) c = cl 4 EIτ , N = Ñl 2 EI , α 1 = 6 g a h 2 (11) R g = g a g d , α 2 = 6 ε 0 l 4 Eh 3 g 3 a , Ω = Ωτ (12)

Reduced order model

A reduced-order model is generated by modal decomposition transforming Equation ( 7) into a discrete system consisting of ordinary differential equations in time. For this end, the Galerkin procedure is used and the deflection is expanded on the undamped linear mode shapes of the straight beam:

w(x, t) = n k=1 a k (t)φ k (x) ( 13
)
where n is the number of retained modes, a k (t) is the k th generalized coordinate and φ k (x) is the k th linear undamped mode shape of the straight non-prestressed beam, normalized such that

1 0 φ k φ j = δ k j where δ k j = 0 if k j and δ k j = 1 if k = j.
The linear undamped mode shapes φ k (x) are governed by:

d 4 φ k (x) dx 4 = λ 2 k φ k (x) (14)
φ k (0) = φ k (1) = φ ′ k (0) = φ ′ k (1) (15)
Let Equation ( 7) be multiplied by

φ k (x) (1 + w)(R g -w)
2 in order to include the complete contribution of the nonlinear electrostatic forces in the resonator dynamics without approximation. This method, which has been used by Nayfeh et al. [START_REF] Nayfeh | Dynamic pull-in phenomenon in mems resonators[END_REF] for a one-port resonator, has some disadvantages like the non orthogonality of the operator w 4 ∂ 4 w ∂x 4 with respect to the undamped linear mode shapes of the resonator, the increase of the nonlinearity level in the normalized equation of motion [START_REF] Lahaye | Approaching the quantum limit of a nanomechanical resonator[END_REF] as well as the incorporation of new nonlinear terms such as the Van der Pol damping. Nevertheless, the resulting equation contains less parametric terms than if the nonlinear electrostatic forces were expanded in Taylor series and the solution of nonlinear problem is valid for large displacements of the beam up to the sensing gap. 4

Produit par HAL -19 Sep 2014

Next, we substitute Equation ( 13) into the resulting equation, use Equation ( 14) to eliminate

d 4 φ k (x)
dx 4 , integrate the outcome from x = 0 to 1, and obtain a system of coupled ordinary differential equations in time.

For several nonlinear mechanical systems such as shells [START_REF] Thouzé | Non-linear behaviour of free-edge shallow spherical shells : Effect of the geometry[END_REF], many modes are needed in order to build the frequency response without losing any physical information concerning the nonlinear coupling between the modes.

Nevertheless, in our case (simple beam resonators), the first mode can be considered as the dominant mode of the system and the higher modes could be neglected which has been demonstrated numerically using time integration [START_REF] Kacem | Nonlinear dynamics of nanomechanical beam resonators: improving the performance of nems-based sensors[END_REF] as well as shooting and harmonic balance method coupled with an asymptotic numerical continuation technique [START_REF] Kacem | Nonlinear dynamics of M&NEMS resonant sensors: design strategies for performance enhancement[END_REF]. Thus, only one mode is considered (n = 1) and the following equation is obtained:

ä1 + cȧ 1 + ω 1 2 a 1 + µ 1 a 1 ä1 + µ 2 a 1 2 ä1 + µ 3 a 1 3 ä1 + µ 4 a 1 4 ä1 + cµ 1 a 1 ȧ1 + cµ 2 a 1 2 ȧ1 +cµ 3 a 1 3 ȧ1 + cµ 4 a 1 4 ȧ1 + χ 2 a 1 2 + χ 3 a 1 3 + χ 4 a 1 4 + χ 5 a 1 5 + χ 6 a 1 6 + χ 7 a 1 7 + ν +ζ 0 cos(Ωt) + ζ 1 a 1 cos(Ωt) + ζ 2 a 1 2 cos(Ωt) + ζ 3 cos(2Ωt) +ζ 4 a 1 cos(2Ωt) + ζ 5 a 1 2 cos(2Ωt) = 0 ( 16
)
where ω 1 is the natural frequency of the first mode and ν is the static force. χ i∈{2,3,4,5,6,7} are the nonlinear stiffnesses and ζ i∈{0,3} are the amplitudes of the time varying linear forces. µ i∈{1,2,3,4} and cµ i∈{1,2,3,4} are the nonlinear acceleration and damping terms which come from the multiplication of Equation ( 7) by the common denominator of the electrostatic forces and ζ i∈{1,2,4,5} are the amplitudes of the parametric forces coming from the coupling between the actuation and sensing electrodes. Equation ( 16) contains canonical nonlinear terms such as nonlinear stiffness (χ 3 a 1 3 ), nonlinear damping (cµ 2 a 1 2 a ′ 1 ) as well as parametric excitation (ζ 1 a 1 cos(Ωt)). Therefore, Equation ( 16) governs the behavior of a Duffing-Van der Pol-Mathieu oscillator. Moreover, the presence of other high-level nonlinearities in Equation ( 16) makes the described system as a forced nonlinear resonator under multifrequency parametric excitation. For a two-port resonator, this kind of equation is similar to the one used by Nayfeh et al. [START_REF] Nayfeh | Dynamic pull-in phenomenon in mems resonators[END_REF] for a one-port resonator and includes terms coming from the coupling between the mechanical and the electrostatic nonlinearities as well as the nonlinear coupling between both electrostatic forces. The expressions of all the integration parameters presented in Equation ( 16) are listed in [START_REF] Kacem | Nonlinear dynamics of nanomechanical beam resonators: improving the performance of nems-based sensors[END_REF]. They can be easily computed with any computational software.

To analyze this equation of motion, we use perturbation techniques which are well adapted to "small" excitation and damping (Q > 10), typically valid in NEMS resonators [START_REF] Husain | Nanowirebased very-high-frequency electromechanical resonator[END_REF]. Since the device operates near resonance, a detuning parameter is introduced as Ω = ω n + εσ and the averaging method [START_REF] Nayfeh | Nonlinear Oscillations[END_REF] is used in order to transform Equation ( 16) into two first order non-linear ordinary-differential equations that modulate the amplitude A and the phase β of the transformation

a 1 = A(t) cos Ωt + β(t) : Ȧ = ε ζ 0 sin β 2ω n + ε A 2 ζ 2 sin β 8ω n -ε Ac 2 -ε A 3 cµ 2 8 -ε A 5 cµ 4 16 + O(ε 2 ) (17) β = εσ -ε 3A 2 χ 3 8ω n -ε 5A 4 χ 5 16ω n -ε 35A 6 χ 7 128ω n -ε ζ 0 cos β 2Aω n -ε 3Aζ 2 cos β 8ω n + ε 3 8 A 2 ω n µ 2 + ε 5 16 A 4 ω n µ 4 + O(ε 2 ) ( 18
)
The steady-state motions occur when Ȧ = β = 0, which corresponds to the singular points of Equations ( 17) and [START_REF] Nishiyama | Capacitance of a strip capacitor[END_REF].Thus, the frequency-response equation can be written in its parametric form {A = K 1 (β), Ω = K 2 (β)} 5

Produit par HAL -19 Sep 2014

as a function of the phase β as a set of 2 equations easy to introduce in Matlab or Mathematica. This ability makes the model suitable for NEMS designers as a quick tool of resonant sensor performance optimization.

For the sake of clarity, the frequency-response equation can be written in its implicit form as:

A 2 {Λ 1 (A) -8ω n (16Ω + Λ 2 (A))} 2 256 4ζ 0 + 3A 2 ζ 2 2 + 8cA + 2A 3 cµ 2 + A 5 cµ 4 2 ω 2 n 4 4ζ 0 + A 2 ζ 2 2 = 1 (19) Λ 1 (A) = 48A 2 χ 3 + 40A 4 χ 5 + 35A 6 χ 7 ( 20
)
Λ 2 (A) = 6A 2 µ 2 + 5A 4 µ 4 -16 ω n ( 21
)

Typical behaviors and physical limitations due to nonlinear phenomena

All the numerical simulations were carried out for a device with the following parameters: l = 20 µm, b = 4 µm, h = 2 µm, l a = 20 µm, l d = 112 µm and an axial residual stress on the beam material σ r = 15 MPa. V s is set to zero volt which has no impact on the generality of the analytical investigations since the electrostatic force is proportional to the square of the electric potential difference. Both gaps can be used to tune the electrostatic nonlinear effect. Nevertheless, for convenience, the gap between the beam and the actuation electrode g a is set to 1 µm. V ac , V dc and g d were used for parametric analysis.

These dimensions are extracted from the design of a resonant accelerometer which has been fabricated in LETI clean rooms in order to investigate the nonlinear dynamics of resonant sensors within the framework of the transition from MEMS to NEMS [START_REF] Kacem | Nonlinear dynamics of nanomechanical beam resonators: improving the performance of nems-based sensors[END_REF].

Hardening behavior

At the micro and nanoscale, the spring hardening is the most classical effect observed in clampedclamped resonators electrostatically actuated [START_REF] Gui | Nonlinearity and hysteresis of resonant strain gauges[END_REF][START_REF] Shao | The nonlinearity cancellation phenomenon in micromechanical resonators[END_REF]. The mechanical nonlinearities due to mid-plane stretching dominate the resonator dynamics and the frequency response peak is hysteretic and shifted to the high frequencies (Figure 2). Furthermore, when the mechanical nonlinearities are preponderant, the dynamics of one and two ports resonators are equivalent.

Softening behavior

In order to increase the softening electrostatic nonlinearities, the resonator designs have to involve very narrow gaps with respect to the beam width which is of great difficulty with Top-Down technology for NEMS. Therefore the softening behavior is difficult to obtain in clamped-clamped NEMS resonators electrostatically actuated. For the considered resonator, the softening behavior is possible for a sensing gap g d = 400 nm as shown in Figure 3. The frequency response curve is hysteretic and shifted to the low frequencies.

Mixed behavior

Both mechanical and electrostatic nonlinearities are always operating into the system. However, in some configurations one kind of nonlinearity is negligible with respect to the second one. Practically, when h/g d << 1, then the dynamics is dominated by the hardening nonlinearities and in the opposite case (h/g d >> 1), the frequency response is nonlinearly softening. Between these two configurations and for the typical fabricated resonator described in Figure 1, the hardening behavior could occur for low DC voltages. Then, when V dc is increased up to 8V, the nonlinear softening stiffness coming from the electrostatic becomes sufficiently large to create a mixed hardening-softening behavior described in Figure 4 except for an optimal gap g d = 500 nm for which both nonlinearities could be perfectly equilibrated [START_REF] Kacem | Nonlinear dynamics of nanomechanical beam resonators: improving the performance of nems-based sensors[END_REF]. The mixed 6

Produit par HAL -19 Sep 2014 behavior is characterized by a four-bifurcation points frequency response and up to 5 amplitude for a given frequency [START_REF] Kacem | Bifurcation topology tuning of a mixed behavior in nonlinear micromechanical resonators[END_REF]. In the present case of Figure 4, the bifurcation point 3 can not be reached since the P point is at a higher frequency than the bifurcation point 1 and thus, at this point, the jump up in sweep down frequency brings the solution at the point d 3 of the lower branch where the basin of attractions are quite large. Hence, the upper branch between d 1 and 3 is unstable. Indeed, the mixed behavior is highly unstable, dangerous and thus undesirable for MEMS and NEMS designers.

Critical amplitude

The critical amplitude is the oscillation amplitude A c above which bistability occurs. Thus, A c is the transition amplitude from the linear to the nonlinear behavior (see Figure 5). At the critical drive, the resonance curve exhibits a point of infinite slope, called the critical point. Moreover, at the same point, the phase curve also exhibits an infinite slope at the same detuning as the resonance curve itself. Nayfeh studied the stability of an excited Duffing oscillator [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. Kaajakari [START_REF] Kaajakari | Nonlinear mechanical effects in silicon longitudinal mode beam resonators[END_REF] provided close form expression for the critical amplitude using a reduced order model including the crystalline direction of beam resonators. However, the complete contribution of the electrostatic nonlinearities was not incorporated. In other words, he simply provided the mechanical critical amplitude.

Mechanical critical amplitude

For the primary resonance of a clamped-clamped beam, A c is defined as the oscillation amplitude for which the equation dσ dβ = 0 has a unique solution. In order to explain how to deduce the mechanical 9 19) written in its parametric form, let the nonlinear electrostatic effects be neglected (h/g d << 1). The parametric form of the frequency response can be written as:

σ = 1 8        3ζ 2 0 χ 3 c 2 ω 3 1 sin 2 β -4c cot β        (22) A = ζ 0 cω n sin β (23)
Using Equation ( 22) permits the differentiation of the detuning parameter σ with respect to the phase β:

dσ dβ = 6χ 3 ζ 2 0 sin 2β -3χ 3 ζ 2 0 sin 4β + 16c 3 ω 3 1 32c 2 ω 3 1 sin 2 β (24)
One can search the condition for which the equation 6χ 3 ζ 2 0 sin 2β -3χ 3 ζ 2 0 sin 4β + 16c 3 ω 3 1 = 0 has solutions. This equation can be written as ∆(2 sin 2βsin 4β) + 1 = 0 where ∆ = 3ζ 2 0 χ 3 16c 3 ω 3 1 and then, transformed to a fourth order polynomial equation by using the change of variable X = sin 2β.

1 + 4∆X + 4∆ 2 X 4 = 0 (25)
Equation ( 25) can be solved using the Ferrari's method and the following cases are obtained:

             I f ∆ < 2 3 √ 3 =⇒ No solutions Linear behavior I f ∆ = 2 3 √ 3 =⇒ A unique solution β c = 2π 3 Critical behavior I f ∆ > 2 3 √ 3 =⇒ 2 distinct solutions β 1 et β 2 Nonlinear behavior
The condition ∆ = 2 3 √ 3 implies the following relation

ζ 0c = 4 √ 2c 3/2 ω 3/2 1 3 4 √ 3 √ χ 3 (26)
ζ 0c corresponds to the critical drive for which the resonance curve exhibits an infinite slop at the phase β c = 2π 3 and the two bifurcation points B 1 and B 2 are superposed. When this condition is satisfied, the mechanical critical amplitude is computed at the resonance peak. It is the value of Equation ( 23) at the phase β = π 2 (for a critical behavior) multiplied by 1.588 which correspond to the amplification coefficient of the first mode at the middle of the clamped-clamped beam x = 1 2 . Substituting Equation [START_REF] Kacem | Bifurcation topology tuning of a mixed behavior in nonlinear micromechanical resonators[END_REF] in the resulting equation and mutiplying by the gap g d to obtain the dimensional value, yields:

A c m = 2.275g d √ c √ ω 1 √ χ 3 (27)
Substituting ω 1 , c and χ 3 expressions in Equation ( 27), yields:

A c m = 1.685 h √ Q ( 28
)
where Q = ρbh ω 1 τc is the quality factor of a resonant clamped-clamped beam under its first bending vibration mode. By using the reduced order model and the first order averaging method, the mechanical critical amplitude was easily deduced and the obtained close-form solution is comparable to those defined in [START_REF] Kaajakari | Nonlinear mechanical effects in silicon longitudinal mode beam resonators[END_REF][START_REF] Shao | The nonlinearity cancellation phenomenon in micromechanical resonators[END_REF]. 11

Produit par HAL -19 Sep 2014

Global critical amplitude

The close-form expression of the mechanical critical amplitude A c m has been provided using the nonlinear model of a clamped-clamped resonator while neglecting the effects of the electrostatic nonlinearities. Indeed, Equation [START_REF] Osterberg | M-test: A test chip for mems material property measurement using electrostatically actuated test structures[END_REF] represents the classical form of the critical amplitude for a Duffing resonator (only mechanical nonlinearity) [START_REF] Osterberg | M-test: A test chip for mems material property measurement using electrostatically actuated test structures[END_REF]. It shows that the critical amplitude is only determined by the beam thickness in the direction of vibration h and the quality factor Q and does not depend on the beam length l. This information has been observed experimentally by Shao et al. [START_REF] Shao | The nonlinearity cancellation phenomenon in micromechanical resonators[END_REF] for micromechanical clamped-clamped beam resonators using stroboscopic SEM. Our model allows the computation of the global critical amplitude when all sources of nonlinearities are included: it can be deduced using the same method as explained in the simplified case for the mechanical critical amplitude.

A c = θ 1 h 2 -θ 2 + θ 3 -θ 4 h 2 + θ 5 h 4 + θ 6 h 2 Q ( 29
)
θ 1 = 3.6327 * 10 -3 (30)
θ 6 = 2.5683g d 2 (34)
Remarkably the mechanical critical amplitude (Equation ( 28)) depends only on the thickness of the beam h and the quality factor Q. However, in addition to these two parameters, the global critical amplitude (Equations ( 29)-( 35)) depends also on the detection gap g d .

It can easily checked that lim

g d →∞ A c = A c m .
For example, the critical amplitude of a resonator having a quality factor of 10 4 designed with 100 nm of thickness in the direction of vibration and a sensing gap of 200 nm, is about 1.68 nm.

In order to stay linear, the restrictive amplitude is 0.84% of the sensing gap, which leads to a very weak signal to noise ratio (SNR) and thus a very poor resolution.

Pull-in

As it can be observed in Equation ( 2), the electrostatic force is inversely proportional to the square of the gap between the beam and the electrode. As the gap decreases, the generated attractive force increases quadratically. The only opposing force to the electrostatic loading is the mechanical restoring force. If the voltage increases, the gap decreases generating an incremented force. At some point the mechanical forces defined by the spring cannot balance this force anymore. Once this state is reached, the beam snaps against the electrode, and in most cases, the system is permanently disabled. Consequently, the electrostatic loading has an upper limit beyond which the mechanical force can no longer resist the opposing electrostatic force, thereby leading to the collapse of the structure. This actuation instability phenomenon is known as pull-in, and the associated critical voltage is called the Pull-in Voltage.

Static pull-in

To derive the expression for pull-in, let us consider a one-port resonator [START_REF] Kacem | From mems to nems: Modelling and characterization of the non linear dynamics of resonators, a way to enhance the dynamic range[END_REF]. Neglecting the mechanical nonlinearities, the total potential energy in the system can be written as follows:

E p = - 1 2 ǫ 0 bl g -u V 2 dc + 1 2 ku 2 ( 36
)
where u is the average displacement of the resonator and g is the gap between the beam and the electrode. The first term in Equation ( 36) is the electrostatic potential of the deformable capacitor (the resonator) and the second term is due to the mechanical energy stored in the spring (k is the effective spring constant of the resonator). The force acting on the movable beam is obtained by deriving Equation (36):

F = - ∂E p ∂u = 1 2 ǫ 0 bl (g -u) 2 V 2 dc -ku (37)
At equilibrium, the electrostatic force and spring force cancels (F = 0) and Equation (37) gives:

ku = 1 2 ǫ 0 bl (g -u) 2 V 2 dc (38)
Equation ( 38) can be solved for the equilibrium beam position w as a function of applied voltage V dc . Above the pull-in voltage V P , Equation (38) has no solutions. A simple expression for the pull-in point is deduced by deriving Equation (37) to obtain the stiffness of the system:

∂F ∂u = ǫ 0 bl (g -u) 3 V 2 dc -k (39)
Substituting Equation (38) gives the stiffness around the equilibrium point:

∂F ∂u = 2ku g -u -k (40)
With no applied voltage Equation ( 40) is simply ∂F ∂u = -k; a virtual positive deflection δu induces a negative restoring force ∂F ∂u δu = -kδu. Increasing the bias voltage V dc makes the stiffness less negative. The unstable point is given by ∂F ∂u = 0 giving u = 1 3 g (41)

Beyond this point the stiffness becomes positive and the system is unstable: a virtual positive deflection δu induces a positive force that increases u. Substituting Equation (41) to Equation (38) gives the pull-in voltage at which the system becomes unstable

V P = 8kg 3 27ǫ 0 bl (42)

Dynamic pull-in

The pull-in amplitude is the oscillation amplitude above which the resonator position becomes unstable and collapses. The dynamic pull-in is the collapse of the beam subjected to a time-varying electrostatic force. It is completely different from the static pull-in [START_REF] Osterberg | M-test: A test chip for mems material property measurement using electrostatically actuated test structures[END_REF] where the electrostatic force depends only on the gap. 13

Produit par HAL -19 Sep 2014

In the general case, pull-in can occur for multistable hardening and softening behaviors. Nevertheless, for practical cases of nanoresonators which are designed with gaps and width in the direction of vibration of the same order of magnitude, the mixed behavior occurs before pull-in is reached.

Consequently, and like other studies [START_REF] Ouakad | Nonlinear dynamics of electrically-actuated carbon nanotube resonator[END_REF], we define point P as the initiation of an unstability domain and thus an upper bound of possible drive, beyond which dynamic pull-in is likely to occur. Thus, the pull-in domain initiation amplitude is located at the bifurcation point P where the effect of the nonlinear electrostatic stiffness becomes significantly important, and where the domain of attraction of stable branches are small, making the jump of the system to these stable branches quite hard physically [START_REF] Nayfeh | Dynamic pull-in phenomenon in mems resonators[END_REF].

Following [START_REF] Kacem | Nonlinear dynamics of nanomechanical beam resonators: improving the performance of nems-based sensors[END_REF], the initiation of the dynamic pull-in domain can be located via the transition from two to three bifurcation points. The third bifurcation point is situated at the phase β = π 2 which corresponds to the initiation of the mixed behavior. The latter is characterized by a slope approaching infinity at point P as shown in Figure 4 and it has been demonstrated analytically and experimentally that the amplitude of this bifurcation point is set by the geometry of the resonator [START_REF] Kacem | Bifurcation topology tuning of a mixed behavior in nonlinear micromechanical resonators[END_REF]. Hence, the following expression of the pull-in voltage V p :

V p = 3lπcω n ∆ -χ 3 + µ 2 ω 2

Conclusion

This work has detailed the development of an analytical model to qualitatively and quantitatively assess the nonlinear phenomena of NEMS resonators. This model includes the main sources of nonlinearities, in particular the electrostatic ones, without approximation and is based on the modal decomposition using the Galerkin procedure combined with the averaging perturbation technique method.

It has been shown that the model is able to capture all the nonlinear phenomena in the resonator dynamics including bistability and multistability. Beyond the occurrence of such phenomena, the sensor frequency stability is extremely altered. Combined with the noise mixing issue [START_REF] Howe Roessig | Nonlinear mixing in surface-micromachined tuning fork oscillators[END_REF][START_REF] Kaajakari | Phase noise in capacitively coupled micromechanical oscillators[END_REF], this leads to drastic losses in the resonator performances. The developed analytical model proves to be a practical tool for NEMS designers since it provides close-form expressions of the critical amplitudes as well as the pull-in domain initiation amplitude.

Figure 1 :

 1 Figure 1: Two-ports resonator layout.

Figure 2 :

 2 Figure 2: Hardening behavior of the typical resonator described in Figure 1. σ is the detuning parameter, W max is the displacement of the beam normalized by the gap g d at its middle point l 2 and {1, 2} are the bifurcation points.

Figure 3 :

 3 Figure 3: Softening behavior of the typical resonator described in Figure 1. σ is the detuning parameter, W max is the displacement of the beam normalized by the gap g d at its middle point l 2 and {1, 2} are the bifurcation points.

Figure 4 :

 4 Figure 4: Mixed behavior of the typical resonator described in Figure 1 for g d = 750 nm, V dc = 8 V and V ac = 10 -2 V dc . σ is the detuning parameter,W max is the displacement of the beam normalized by the gap g d at its middle point l 2 and {1, 2, 3, P} are the different bifurcation points.

Figure 5 :

 5 Figure 5: Forced frequency responses of the resonator. f a is the dimensionless frequency and W max is the displacement of the beam normalized by the gap g d at its middle point l 2 . A c is the mechanical critical amplitude and {B 1 , B 2 } are the two bifurcation points of a typical hardening behavior.

θ 5 =

 5 1.3197 * 10 -5

n 1 2 χ 5 -µ 4 ω 2 2 C∆ = 200χ 2 5 + µ 4 -21µ 2 2 + 200µ 4 ω 4 n -30µ 2 χ 3 χ 5 +2 6µ 2 µ 4 χ 3 + 5 3µ 2 2 -40µ 4 χ 5 ω 2 n + 9µ 4 χ 2 3

 152252453222 n1 l 2 cos π l sin πl d l 2 -πl d (43)

Produit par HAL -19 Sep 2014

Acknowledgments

The authors are indebted of the Carnot Institute CEA-LETI and I@L for their supports.


N Kacem

S Baguet

S Hentz

R Dufour

Computational and quasi-analytical models for nonlinear vibrations of resonant MEMS and NEMS sensors

Keywords: MEMS, resonator, microbeam, nonlinear dynamics, reduced-order model, harmonic balance method, asymptotic numerical method

Large-amplitude nonlinear vibrations of micro and nano electromechanical resonant sensors around their primary resonance are investigated. A comprehensive multiphysics model based on the Galerkin decomposition method coupled with the averaging method is developed in the case of electrostatically actuated clamped-clamped resonators. The model is purely analytical and includes the main sources of nonlinearities as well as fringing field effects. The influence of the higher modes and the validation of the model are demonstrated with respect to the shooting method as well as the harmonic balance coupled with the asymptotic numerical method. This model allows designers to investigate the sensitivity variation of resonant sensors in the nonlinear regime with respect to the electrostatic forcing.

Introduction

Today, micro-electro-mechanical systems (MEMS) are in the heart of new technologies which are widely used in aerospace, automotive, biotechnology, instrumentation, robotics or manufacturing. Overcoming current technological challenges requires more and more reduced sizes. The next step in miniaturization is the so-called nano-electro-mechanical systems (NEMS). Operating at the nanoscale, these systems come with many interesting physical attributes which makes them the best candidates for ultrasensitive mass and force detection [START_REF] Yang | Zeptogram-scale nanomechanical mass sensing[END_REF][START_REF] Li | Ultra-sensitive nems-based cantilevers for sensing, scanned probe and very highfrequency applications[END_REF][START_REF] Jensen | An atomic-resolution nanomechanical mass sensor[END_REF]. The transition from MEMS to NEMS has implied many reconsiderations concerning the sensing techniques. Besides, in order to compensate the loss of performances when sensors are scaled down to the NEMS level, the resonant sensing has been widely implemented in nanosensors. This type of sensing basically consists in detecting a frequency shift due to an applied axial force to be measured while electrostatically exciting the resonator at its fundamental frequency.

The resolution of such a sensor is given by the minimum detectable frequency shift which is inversely proportional to the resonator drive amplitude. It is limited by the onset of nonlinearities which occur proportionally sooner with respect to the device size [START_REF] Kacem | Nonlinear dynamics of nanomechanical beam resonators: improving the performance of nems-based sensors[END_REF][START_REF] Kacem | Dynamic range enhancement of nonlinear nanomechanical resonant cantilevers for high sensitive nems gas/mass sensors applications[END_REF]. Moreover, below the bistability limit, it is extremely difficult to detect the oscillations of such small sensors. Hence, it is important to model the nonlinear dynamics of NEMS-based resonant sensors at large amplitudes.

The dynamic behavior of MEMS has been investigated by several authors so far. A survey of the electrostatic force modeling, a comprehensive review on MEMS modeling and an overview of the nonlinear dynamics of MEMS can be found in [START_REF] Zhang | Stability, nonlinearity and reliability of electrostatically actuated mems devices[END_REF], [START_REF] Batra | Review of modeling electrostatically actuated microelectromechanical systems[END_REF] and [START_REF] Lifshitz | Reviews of Nonlinear Dynamics and Complexity[END_REF]. Existing models are of variable complexity, depending on assumptions concerning the mechanical and electrostatic nonlinearities. For example a simple lumped spring-mass model was considered in [START_REF] Pamidighantam | Pull-in voltage analysis of electrostatically actuated beam structures with fixed-fixed and fixed-free end conditions[END_REF][START_REF] Shao | The nonlinearity cancellation phenomenon in micromechanical resonators[END_REF] and the continuous linear Euler-Bernoulli beam theory was used in [START_REF] Tilmans | Electrostatically driven vacuum-encapsulated polysilicon resonators: Part ii. theory and performance[END_REF]. Nevertheless, several studies pointed out the importance of deformation dependent mid-plane stretching for axially constrained microbeams undergoing large-amplitude vibration [START_REF] Abdel-Rahman | Characterization of the mechanical behavior of an electrically actuated microbeam[END_REF] and the nonlinear Euler-Bernoulli beam theory is therefore used in most recent studies [START_REF] Younis | A reduced-order model for electrically actuated microbeam-based mems[END_REF][START_REF] Younis | A study of the nonlinear response of a resonant microbeam to an electric actuation[END_REF][START_REF] Xie | Nonlinear dynamic analysis of mems switches by nonlinear modal analysis[END_REF][START_REF] Kuang | Dynamic characteristics of shaped micro-actuators solved using the differential quadrature method[END_REF][START_REF] Najar | Modeling and design of variable-geometry electrostatic microactuators[END_REF][START_REF] Nayfeh | Dynamic pull-in phenomenon in mems resonators[END_REF][START_REF] Krylov | Lyapunov exponents as a criterion for the dynamic pull-in instability of electrostatically actuated microstructures[END_REF][START_REF] Batra | Vibrations of narrow microbeams predeformed by an electric field[END_REF][START_REF] Moghimi Zand | Semi-analytic solutions to nonlinear vibrations of microbeams under suddenly applied voltages[END_REF]. The electrostatic force is usually modeled using the parallel-plate approximation together with correction terms accounting for the fringing field effects due to beam's finite width [START_REF] Nishiyama | Capacitance of a strip capacitor[END_REF]. In order to facilitate the use of numerical or analytical methods, the nonlinearity of the electrostatic forcing term was further simplified by means of a high-order Taylor expansion in [START_REF] Shao | The nonlinearity cancellation phenomenon in micromechanical resonators[END_REF][START_REF] Younis | A reduced-order model for electrically actuated microbeam-based mems[END_REF][START_REF] Moghimi Zand | Semi-analytic solutions to nonlinear vibrations of microbeams under suddenly applied voltages[END_REF][START_REF] Chao | Dc dynamic pull-in predictions for a generalized clamped-clamped microbeam based on a continuous model and bifurcation analysis[END_REF] or a least squares approximation in [START_REF] Xie | Nonlinear dynamic analysis of mems switches by nonlinear modal analysis[END_REF].

Except for simplified models [START_REF] Pamidighantam | Pull-in voltage analysis of electrostatically actuated beam structures with fixed-fixed and fixed-free end conditions[END_REF][START_REF] Shao | The nonlinearity cancellation phenomenon in micromechanical resonators[END_REF][START_REF] Batra | Vibrations of narrow microbeams predeformed by an electric field[END_REF], analytical solutions are usually not obtainable for the continuous model. Some authors used finite-element simulations to study the vibrational behavior of MEMS either by means of a specific fully coupled finite element formulation [START_REF] Rochus | Electrostatic coupling of mems structures: transient simulations and dynamic pull-in[END_REF] or a commercial software [START_REF] Batra | Vibrations of narrow microbeams predeformed by an electric field[END_REF]. Reduced-order models are generally preferred due to their limited computational cost. The Galerkin method is used in [START_REF] Younis | A reduced-order model for electrically actuated microbeam-based mems[END_REF][START_REF] Xie | Nonlinear dynamic analysis of mems switches by nonlinear modal analysis[END_REF][START_REF] Nayfeh | Dynamic pull-in phenomenon in mems resonators[END_REF][START_REF] Krylov | Lyapunov exponents as a criterion for the dynamic pull-in instability of electrostatically actuated microstructures[END_REF][START_REF] Chao | Dc dynamic pull-in predictions for a generalized clamped-clamped microbeam based on a continuous model and bifurcation analysis[END_REF], the simplified 1-mode Galerkin method in [START_REF] Batra | Vibrations of narrow microbeams predeformed by an electric field[END_REF][START_REF] Moghimi Zand | Semi-analytic solutions to nonlinear vibrations of microbeams under suddenly applied voltages[END_REF] and the differential quadrature method (DQM) in [START_REF] Kuang | Dynamic characteristics of shaped micro-actuators solved using the differential quadrature method[END_REF][START_REF] Najar | Modeling and design of variable-geometry electrostatic microactuators[END_REF]. Several techniques have been proposed for solving the resulting nonlinear equations of motion, among which time integration [START_REF] Younis | A reduced-order model for electrically actuated microbeam-based mems[END_REF][START_REF] Krylov | Lyapunov exponents as a criterion for the dynamic pull-in instability of electrostatically actuated microstructures[END_REF][START_REF] Chao | Dc dynamic pull-in predictions for a generalized clamped-clamped microbeam based on a continuous model and bifurcation analysis[END_REF][START_REF] Rochus | Electrostatic coupling of mems structures: transient simulations and dynamic pull-in[END_REF], perturbation methods [START_REF] Younis | A study of the nonlinear response of a resonant microbeam to an electric actuation[END_REF], the shooting method [START_REF] Nayfeh | Dynamic pull-in phenomenon in mems resonators[END_REF], the invariant manifold method (or method of nonlinear normal forms) [START_REF] Xie | Nonlinear dynamic analysis of mems switches by nonlinear modal analysis[END_REF] or the homotopy analysis method which provides semi-analytic solutions [START_REF] Moghimi Zand | Semi-analytic solutions to nonlinear vibrations of microbeams under suddenly applied voltages[END_REF].

In this paper, the nonlinear dynamics of micro and nano electromechanical resonant sensors are investigated at large amplitudes around their primary resonance. A comprehensive multiphysics model accounting for main sources of nonlinearities as well as fringing field effects and based on the Galerkin discretization is developed in the case of electrostatically actuated clamped-clamped resonators. An efficient numerical procedure for the continuation of periodic solutions based on the harmonic balance method (HBM) coupled with the asymptotic numerical method (ANM) is used to obtain the response curves of the resonator. A reduced-order purely analytical solution is then derived and validated with respect to the numerical HBM+ANM solution for several designs and working parameters. Finally, the interest of such an analytical solution is demonstrated though a study of the sensitivity of the resonator to axial force for acceleration sensing purpose.

Model

A clamped-clamped micro-beam as sketched in Fig. 1 is considered in the present study.

It is actuated by a voltage v( t) = V dc + V ac cos( Ωt ), where V dc is the DC polarization voltage, V ac is the amplitude of the applied AC voltage, and Ω is the forcing frequency. A distributed electrostatic force is created between the beam and the rigid electrode, thus leading to the microbeam deformation and vibration.

Equation of motion

The equation of motion that governs the transverse deflection w( x, t) along Z axis is written as

ρbh ∂ 2 w ∂ t2 + c∂ w ∂ t + EI ∂ 4 w ∂ x4 -Ñ(t) ∂ 2 w ∂ x2 = Fe (1)
where x is the coordinate along the length l of the micro-beam which is assumed sufficiently narrow, so that the cylindrical deformations can be neglected, t is time, E and ρ are the effective

Ñ(t) = Ñ0 + Ebh 2l l 0 ∂ w( x, t) ∂ x 2 d x (2)
which consists of a constant applied tensile axial force Ñ0 and a nonlinear part due to stretching effects [START_REF] Michon | Parametric instability of an axially moving belt subjected to multi-frequency excitations: experiments and analytical validation[END_REF]. Ñ0 results from the residual stress on the silicon or the effect of the applied axial load which controls the natural frequency of the micro-beam. The right-hand side of Eq. ( 1) corresponds to the parallel-plate approximation of the electric force including the fringing field effects with respect to the beam width by means of the coefficient C n = 1 + 1.9861 g b 0.8258 [START_REF] Nishiyama | Capacitance of a strip capacitor[END_REF] which has been validated using three-dimensional finite element simulations [START_REF] Kacem | Nonlinear dynamics of nanomechanical beam resonators: improving the performance of nems-based sensors[END_REF], where the assumption of neglected effect of the resonator thickness with respect to its transverse displacement has been verified.

Fe ( x, t) = 1 2 ε 0 bC n V dc + V ac cos( Ωt ) 2 (g -w) 2 (3)
The beam is subject to the boundary conditions

w(0, t) = w(l, t) = ∂ w ∂ x (0, t) = ∂ w ∂ x (l, t) = 0 ( 4 )

Normalization

It is convenient to rescale equations by introducing the non-dimensional variables

w = w g , x = x l , t = t τ (5)
where τ = 2l 2 h 3ρ E . Substituting Eq. ( 5) into Eqs. ( 1) and (4), we obtain

∂ 2 w ∂t 2 + c ∂w ∂t + ∂ 4 w ∂x 4 - ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ N 0 + α 1 1 0 ∂w ∂x 2 dx ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ∂ 2 w ∂x 2 = α 2 [V dc + V ac cos(Ωt)] 2 (1 -w) 2 (6) w(0, t) = w(1, t) = ∂w ∂x (0, t) = ∂w ∂x (1, t) = 0 ( 7 )
where

c = cl 4 EIτ , N 0 = Ñ0 l 2 EI , α 1 = 6 g h 2 , α 2 = 6ε 0 C n l 4 Eh 3 g 3 , Ω = Ωτ (8)

Modal decomposition

The beam deflection w(x, t) can be written as a sum of a static DC displacement w s (x) and a time-varying AC displacement w d (x, t). However, for the considered devices, it is easy to check that the static deflection is negligible. Typically, the measured quality factors Q are in the range [10 4 -5.10 4 ] and V dc ≤ 200V ac . Thus, the static/dynamic deflections ratio is

w s (x) w d (x, t) ≈ V dc 2Q.V ac ≤ 1% (9)
The Galerkin decomposition method based on linear undamped mode shapes is used to eliminate the spatial dependence and transform Eq. ( 6) into a multiple degree-of-freedom system (MDOF) consisting of ordinary differential equations in time. The deflection of the beam is approximated by the series

w(x, t) = N m k=1 a k (t)φ k (x) (10)
where a k (t) is the k th time varying generalized coordinate and φ k (x) is the k th linear undamped eigenmode solution of

d 4 φ k (x) dx 4 = λ 4 k φ k (x) (11)
Thus, the analytical form of the eigenmodes is given by

φ k (x) = A k cos λ k x -cosh λ k x + cosh λ k -cos λ k sin λ k -sinh λ k [sin λ k x -sinh λ k x] ( 12
)
with λ k solutions of the transcendental equation

1 -cos λ k cosh λ k = 0 ( 13
)
These functions are a modal basis for the scalar product

u, v = 1 0 u(x)v(x)dx ( 14
)
and the coefficients A k are chosen to normalize the eigenmodes such that φ i , φ j = δ i j .

Galerkin procedure

In many models in the literature, the electrostatic nonlinear forcing term is approximated by means of a Taylor development in order to simplify the Galerkin procedure [START_REF] Younis | A reduced-order model for electrically actuated microbeam-based mems[END_REF][START_REF] Xie | Nonlinear dynamic analysis of mems switches by nonlinear modal analysis[END_REF]. However, when vibration amplitudes become large, such a development is no longer valid. Here, the complete contribution of the nonlinear electrostatic forces is included in the resonator dynamics without approximation since Eq. ( 6) is multiplied by

φ i (x)(1 -w) 2 .
This method has some disadvantages like the non orthogonality of the operator w 4 ∂ 4 w ∂x 4 with respect to the undamped linear mode shapes of the resonator, the increase of the nonlinearity level in the normalized equation of motion [START_REF] Zhang | Stability, nonlinearity and reliability of electrostatically actuated mems devices[END_REF] as well as the incorporation of new nonlinear terms such as the Van der Pol damping. Nevertheless, the resulting equation contains less parametric terms than if the nonlinear electrostatic forces were expanded in Taylor series and the solution of nonlinear problem is valid for large displacements of the beam up to the gap. The modal projection consists in substituting Eq. [START_REF] Shao | The nonlinearity cancellation phenomenon in micromechanical resonators[END_REF] in Eq. ( 6), using Eq. ( 11) to eliminate d 4 φ k (x)/dx 4 and integrating the outcome from x = 0 to 1. Doing so, Eq. ( 6) becomes äi -2

N m j=1 N m k=1 1 0 φ k φ j φ i dx a k ä j + N m j=1 N m k=1 N m l=1 1 0 φ l φ k φ j φ i dx a l a k ä j +c i ȧi -2 N m j=1 N m k=1 1 0 φ k φ j φ i dx c j a k ȧ j + N m j=1 N m k=1 N m l=1 1 0 φ l φ k φ j φ i dx c j a l a k ȧ j +λ 4 i a i -2 N m j=1 N m k=1 1 0 φ k φ j φ i dx λ 4 j a k a j + N m j=1 N m k=1 N m l=1 1 0 φ l φ k φ j φ i dx λ 4 j a l a k a j (15) - ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ N 0 + α 1 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ N m m=1 N m n=1 1 0 φ n φ m dx a n a m ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ N m j=1 1 0 φ j φ i dx a j -2 N m j=1 N m k=1 1 0 φ k φ j φ i dx a k a j + N m j=1 N m k=1 N m l=1 1 0 φ l φ k φ j φ i dx a l a k a j =α 2 (V dc + V ac cosΩt) 2 1 0 φ i dx for i = 1, 2, . . . , N m
Eq. ( 15) can be written in matrix-vector form as

[M 0 + M 1 (a) + M 2 (a)] ä + [C 0 + C 1 (a) + C 2 (a)] ȧ + [K 0 + K 1 (a) + K 2 (a)] a (16)
-

[N 0 + α 1 T 2 (a)] [K T + K T 1 (a) + K T 2 (a)] a = α 2 (V dc + V ac cosΩt) 2 F with a(t) = [a 1 (t), a 2 (t), . . . , a N m (t)] T . The components of matrices M 0 , M 1 ,M 2 , C 0 , C 1 ,C 2 , K 0 , K 1 ,K 2 , K T , K T 1 and K T 2 are respectively M i j , M 1i j , M 2i j , C i j , C 1i j , C 2i j , K i j , K 1i j , K 2i j , K T i j , K T 1i j , K T 2i j with M 0i j = δ i j M 1i j = -2 N m k=1 1 0 φ k φ j φ i dx a k M 2i j = N m k=1 N m l=1 1 0 φ l φ k φ j φ i dx a l a k C 0i j = c i δ i j C 1i j = c j M 1i j C 2i j = c j M 2i j K 0i j = λ 4 i δ i j K 1i j = λ 4 j M 1i j K 2i j = λ 4 j M 2i j K T i j = 1 0 φ j φ i dx K T 1i j = -2 N m k=1 1 0 φ k φ j φ i dx a k K T 2i j = N m k=1 N m l=1 1 0 φ l φ k φ j φ i dx a l a k ( 17
)
The scalar T 2 (a) and the entries of vector F are respectively

T 2 (a) = N m m=1 N m n=1 1 0 φ n φ m dx a n a m F i = 1 0 φ i dx (18)
Indices 1 and 2 for matrices M, C, K, K T and for T 2 denote nonlinearities of order 1 and 2 with respect to a. Once all these coefficients have been computed, system ( 16) can be solved with a suitable numerical procedure.

High order Harmonic Balance Method combined with the Asymptotic Numerical Method (HBM+ANM)

The Harmonic Balance Method (HBM) is commonly used for computing periodic solutions.

It consists in assuming a time solution in the form of a Fourier series and comparing/balancing the coefficients of the same harmonic components. In this way, non linear differential equations in the space variables and time are transformed into a nonlinear algebraic system in the space variables and frequency. However, when nonlinearities are complicated, the derivation of the algebraic system becomes very cumbersome. Alternative methods have been proposed to overcome these shortcomings, such as the incremental harmonic balance method (IHBM) [START_REF] Lau | Amplitude incremental variational principle for nonlinear vibration of elastic systems[END_REF] or the alternating frequency/time domain harmonic balance method (AFT) [START_REF] Nacivet | A dynamic lagrangian frequency-time method for the vibration of dry-friction-damped systems[END_REF] but they are very demanding from a computational point of view.

Recently, Cochelin et al. [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF] have proposed another strategy for applying the classical HBM with a large number of harmonics. The basic idea consists in recasting the original system [START_REF] Kuang | Dynamic characteristics of shaped micro-actuators solved using the differential quadrature method[END_REF] into a new system where nonlinearities are at most quadratic polynomials by introducing as many new variables as needed. This leads to an augmented, but quadratic only, nonlinear system for which the application of the HBM is quite straightforward. Furthermore, this quadratic framework makes it possible to use the so-called Asymptotic Numerical Method (ANM) for the continuation of solutions. The ANM consists in computing power series expansions of solution branches and presents several advantages: it provides continuous solutions, the continuation is very robust, and the control of the step length is automatic and always optimal [START_REF] Azrar | An asymptotic-numerical method to compute the postbuckling behaviour of elastic plates and shells[END_REF][START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF][START_REF] Baguet | On the behaviour of the anm in the presence of bifurcations[END_REF][START_REF] Cochelin | Méthode asymptotique Numérique, Collection méthodes numériques[END_REF].

This method is detailled in [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF]. Its application to our nonlinear differential system ( 16) is detailled hereafter. The association of the HBM and the ANM will be subsequently referred to as HBM+ANM.

Quadratic recast

The key point of this method lies in the quadratic recast of Eq. ( 16) by introducing the following set of auxiliary variables,

y = ȧ (size N m ) z = ä = ẏ (size N m ) M tot = M 1 (a) + M 2 (a) ( s i z e N 2 m ) (19)
K T tot = K T 1 (a) + K T 2 (a) ( s i z e N 2 m ) S = K T tot a (size N m ) T = T 2 (a) ( s i z e 1 )
system ( 16) can be rewritten as

ȧ = y ẏ = z 0 = M tot -M 1 (a) -M 2 (a) 0 = K T tot -K T 1 (a) -K T 2 (a) 0 = S -K T tot a 0 = T -T 2 (a) 0 = -α 2 (V dc + +M 0 z + C 0 y +M tot (z + C 0 y + K 0 a) V ac cosΩt) 2 F +K 0 a -N 0 K T a -N K T tot a -α 1 T a -N 0 S -α 1 T S m( Ẋ) = l 0 (Ω) + l(X) + q(X, X) (20)
where X = (a, y, z, M tot , K T tot , S, T ) T is the unknown vector of size N eq = 2N 2 m + 4N m + 1, in which matrices M tot and K T tot are reshaped as vectors. l 0 is a constant vector with respect to X, l(•) and m(•) are linear vector valued operators with respect to X, and q(•, •) is a quadratic vector valued operator.

Application of the HBM

The HBM consists in decomposing X(t) into a truncated Fourier series

X(t) = X 0 + H k=1 X c k cos (kΩt) + H k=1 X s k sin (kΩt) (21)
and inserting this expansion in the nonlinear differential system [START_REF] Batra | Vibrations of narrow microbeams predeformed by an electric field[END_REF]. By balancing the first 2H + 1 harmonic terms and neglecting the higher order harmonics, the following nonlinear algebraic system of size (2H + 1) × N eq is obtained

Ω M(U) = L 0 + L(U) + Q(U, U) (22)
where

U = [X 0 , X c 1 T , X s 1 T , . . . , X c H T , X s H T ]
T contains the components of the Fourier series [START_REF] Moghimi Zand | Semi-analytic solutions to nonlinear vibrations of microbeams under suddenly applied voltages[END_REF]. The new operators L 0 , M(•), L(•) and Q(•, •) depend only on the operators l 0 , m(•), l(•) and q(•, •). Their explicit expressions can be found in Cochelin et al. [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF].

Continuation by the ANM

In order to apply the Asymptotic Numerical Method, Eq. ( 22) is reformulated as

R(U, Ω) = L 0 + L(U) + Q(U, U) -Ω M(U) = 0 ( 23
)
Assuming that a regular solution (U 0 , Ω 0 ) of Eq. ( 22) is known, the branch of solution starting at (U 0 , Ω 0 ) is represented as a power series expansion (truncated at order N) with respect to the path-parameter η [32]

U(η) = U 0 + ηU 1 + η 2 U 2 + . . . + η N U N Ω(η) = Ω 0 + ηΩ 1 + η 2 Ω 2 + . . . + η N Ω N ( 24
)
If the pseudo-arclength definition is chosen, the path-parameter a can be expressed as [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF]

η = U 1 .(U -U 0 ) T + Ω 1 (Ω -Ω 0 ) (25)
where (U 1 , Ω 1 ) is the tangent vector at (U 0 , Ω 0 ). Introducing the series expansions (24) into Eqs.

(23) and ( 25) and keeping the power terms up to order N leads to

R(η) = R 0 + ηR 1 + η 2 R 2 + . . . + η N R N = 0 ( 26
)
Equating the terms R i (1≤i≤N) to zero permits transforming the nonlinear system [START_REF] Chao | Dc dynamic pull-in predictions for a generalized clamped-clamped microbeam based on a continuous model and bifurcation analysis[END_REF] into N successive linear systems of N eq equations, that are then solved recursively in order to obtain U i and Ω i

(1 ≤ i ≤ N). This is very efficient from a computational point of view since all the linear systems share the same matrix which corresponds to the Jacobian of R evaluated at (U 0 , Ω 0 ). Only the right-hand side vector changes with the order. Moreover, the range of validity a max of the series can be approximated a priori by

η max = ε R N+1 1 N+1 (27)
where ε is a user-defined tolerance parameter. Thus a part of the solution branch is obtained by following Eq. ( 24) until η = η max . This end point is then used as a new starting point (U 0 , Ω 0 ) and the next part of the solution branch is obtained by restarting the continuation process. As a consequence, the step length is naturally adaptive and optimal [START_REF] Baguet | On the behaviour of the anm in the presence of bifurcations[END_REF] and the continuation is very robust.

Stability analysis

A stability analysis can be conducted during the continuation by studying the complex eigenvalues (Floquet multipliers) of the Monodromy matrix M associated with the periodic solution of system [START_REF] Batra | Vibrations of narrow microbeams predeformed by an electric field[END_REF]. This Monodromy matrix is obtained by time integration over one period of the differential system

Ṁ = J M ( 28
)
where J is the Jacobian matrix of [START_REF] Batra | Vibrations of narrow microbeams predeformed by an electric field[END_REF] and with the identity matrix as initial condition. The solution is stable if all the multipliers lie inside the complex unit circle. It is unstable otherwise [START_REF] Nayfeh | Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods[END_REF][START_REF] Dufour | Parametric instability of a beam due to axial excitations and to boundary conditions[END_REF].

Simplified analytical model

Mode properties and approximate integrals

Using the change of variable z = x -1 2 for all the eigenmodes, it is easy to check that Δ : z = 0 represents an axis of symmetry for all the odd modes and an axis of antisymmetry for all the even modes. Thus, the odd modes are even functions and the even modes are odd functions ∀ z ∈ - 1 2 , 1 2 . Consequently:

∀ k ∈ N * , ∀ i ∈ N * and ∀ z ∈ - 1 2 , 1 2
1 0 Φ 2i dx = 1 0 Φ 2i Φ k 2 dx = 1 0 Φ 2i Φ 3 2k-1 dx = 0 1 0 Φ 2i-1 Φ 3 2k dx = 1 0 Φ 2i Φ k Φ k dx = 0 1 0 Φ 2i Φ 2k-1 dx = 1 0 Φ 2i-1 Φ 2k dx = 0 1 0 Φ 2i Φ 2 2k-1 Φ 2k-1 dx = 1 0 Φ 2i-1 Φ 2 2k Φ 2k dx = 0 ( 29
)
As a first approximation, the eigenvalues λ k which are solutions of the transcendental Eq. ( 13), can be written as:

∀k ≥ 1 λ k [2k + 1]π 2 (30)
Consequently:

∀k ≥ 1 cos[λ k ] sin [2λ k ] 0 sin [λ k ] [-1] k cos [2λ k ] -1 (31)
Moreover, cosh [λ k ] ∼ sinh [λ k ] in [λ 1 , +∞].
One can check that the validity of this equivalence starts from the first mode. In fact, for k = 1, λ k 4.73 and cosh [λ 1 ] sinh [λ 1 ] 56.65.

In order to simplify the different integrals of Eqs. ( 17)-( 18), Eqs. ( 30) and ( 31) are used combined with the equivalence cosh [λ k ] ∼ sinh [λ k ]. Then, these integrals are approximated as equivalent to their limit when the unbounded function cosh [λ k ] → +∞. Thus, approximate analytical closedform expressions of the different integrals are obtained with respect to the mode number (see Appendix).

Reduced order model

Assuming that the first mode should be the dominant mode of the system, only one mode is retained (N m = 1) and the other modes are neglected. Equation ( 16) becomes :

ä1 + (500.564 + 12.3N 0 )a 1 + (927

+ 28N 0 + 151α 1 ) a 1 3 +347α 1 a 1 5 + (1330.9 + 38.3N 0 )a 1 2 + 471α 1 a 1 4 +2.66c 1 a 1 ȧ1 + 1.85c 1 a 1 2 ȧ1 + c 1 ȧ1 + 2.66a 1 ä1 +1.85a 1 2 ä1 = - 8 3π α 2 [V dc + V ac cos(Ωt)] 2 (32)
To analyse the equation of motion [START_REF] Cochelin | Méthode asymptotique Numérique, Collection méthodes numériques[END_REF], it proves convenient to invoke perturbation techniques which work well with the assumptions of "small" excitation and damping, typically valid in MEMS resonators. To facilitate the perturbation approach, in this case the method of averaging [START_REF] Nayfeh | Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods[END_REF], a standard constrained coordinate transformation is introduced by assuming a slowly time-varying amplitude:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ a 1 = A(t) cos Ωt + β(t) ȧ1 = -A(t)Ω sin Ωt + β(t) ä1 = -A(t)Ω 2 cos Ωt + β(t) (33)
In addition, since near-resonant behavior is the principal operating regime of the proposed system, a detuning parameter σ is introduced, as given by:

Ω = ω n + εσ (34)
Separating the resulting equations and averaging them over the period 2π Ω in the t-domain results in the system's averaged equations, in terms of amplitude and phase, which are given by:

Ȧ = - 1 2 εξ 0 A - 1 8 εξ 2 A 3 - 1 2 ε κ ω n sin β + O(ε 2 ) ( 35
)
A β = Aσε - 3 8 ε χ 3 ω n A 3 - 5 16 ε χ 5 ω n A 5 + 7 10 εω n A 3 + 1 2 ε κ ω n cos β + O(ε 2 ) ( 36
)
where ω n = 500.564 + 12.3N 0 [START_REF] Kacem | Bifurcation topology tuning of a mixed behavior in nonlinear micromechanical resonators[END_REF] and the other parameters are ξ 0 = c 1 , ξ 2 = 1.85c 1 , χ 3 = 927 + 28N + 151α 1 , χ 5 = 347α 1 and κ = 16 3π α 2 V ac V dc . The steady-state motions occur when Ȧ = β = 0, which corresponds to the singular points of Eqs. ( 35) and [START_REF] Sazonova | A tunable carbon nanotube resonator[END_REF]. Thus, the frequency-response equation can be written in its implicit form as:

3χ 3 4ω n A 2 + 5χ 5 8ω n A 4 - 7ω n 5 A 2 -2σ 2 + ξ 0 + ξ 2 4 A 2 2 = κ Aω n 2 (38)
The normalized displacement W max with respect to the gap at the middle of the beam and the forcing frequency Ω can be expressed in function of the phase β. Thus, the frequency response curve can be plotted parametrically.

This analytical model is a lightweight, easy-to-use and very efficient tool for MEMS designers.

Although it is very simple, it is able to capture all the nonlinear phenomena in the resonator dynamics (hardening and softening behaviors) and describe the competition between them.

Results and discussion

For all the numerical simulations with the ANM+HBM, only symmetric modes were considered due to the symmetry of the problem. It was verified also numerically that antisymmetric modes do not participate in the global dynamics of the resonator. In other words, a computation with 3 modes uses the odd modes 1, 3 and 5. Two resonator designs were used, which have different size and quality factor Q as listed in Table 1. The ratio between the AC and DC voltages was set at V ac = 0.1V dc for all simulations.

Resonator L(μm) b(μm) h(μm) g(μm) Q

Design 1 400 10 10 2 10000

Design 2 50 1 1 0.4 1000

Table 1: Design parameters of investigated resonators.

Validation of the numerical solution

The dynamic response of the resonator obtained by HBM+ANM continuation is plotted in Fig. 2 for design 1 and an actuation voltage V dc =9V, thus V ac =0.9V. W max is the displacement of the beam normalized by the gap g at its middle point x = l/2. For this set of parameters, the resonator exhibits a strongly hardening behavior. The solution was computed with N m =1 mode for the Galerkin procedure, H=7 harmonics for the HBM and the series truncated at order N=30 for the ANM. It can be observed that the high order polynomial expansions of the ANM continuation permit the description of the complete resonance curve with only 15 continuation steps and going through the sharp turning point with no difficulty. Unlike the classical continuation methods based on the Newton-Raphson prediction-correction procedure, the solution curve can be accurately evaluated and plotted between two consecutive ends of step. Moreover, the stability analysis can also be conducted for any value of the continuation parameter η during one step and therefore the transition between stable and unstable branches can be precisely observed even if it occurs between two consecutive ends of step.

In order to validate the HBM+ANM, the shooting method was also implemented to be used as a reference solution. For the sake of conciseness, the implementation of the well-known shooting method is not detailed herein but can be found in [START_REF] Nayfeh | Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods[END_REF] for instance. The confrontation is shown in Figs. 3 and4 on design 1 described in Table 1 for electrostatic forcing leading to slightly (V dc =5V, V ac =0.5V) and strongly (V dc =9V, V ac =0.9V) nonlinear behaviors respectively. For each case, the number of modes N m retained in the simulation and the number of harmonics H used for the HBM are varied in order to study both the influence of N m and the convergence of the HBM+ANM according to H. In the slightly nonlinear case of Fig. 3, the convergence of the HBM+ANM is achieved with only 3 harmonics and results are in very good agreement with the shooting method.

The slight difference between the curves obtained by shooting and HBM+ANM is due to the numerical precision used for the time-integration scheme of the shooting procedure. Concerning the number of modes, the results are identical with 2, 3 or 4 modes and the difference between 1 and 2 modes is negligible when compared to the peak frequency. In the strongly nonlinear case of Fig. 4, a small difference between the curves for 1 mode and 2, 3 or 4 modes is noticeable. Nevertheless, it is less than 0.1% with respect to the peak frequency, which is negligible compared to the frequency shifts induced by the fabrication tolerances. Again, the results are in good agreement between HBM+ANM and shooting methods, and 3 harmonics are sufficient to achieve a very good precision. It is important to underline that despite the higher number of equations to be solved (255 equations for the HBM+ANM instead of 2 equations for the shooting in the case N m =2 H=7), the HBM+ANM is much faster than shooting and its continuation also proves to be more robust as it does not require the development of specific adaptive step length strategies.

Validation of the analytical model

As previously demonstrated with the HBM+ANM, the errors made by using only the first mode in the Galerkin method are negligible with respect to errors due to others constraints such as fabrication, temperature, pressure. . . Thus one can expect that the reduced order analytical model should give satisfying results. For this purpose, the validation is ensured through comparisons with the HBM+ANM results.

Figure 5 shows a first confrontation between both methods on design 1 of Table 1 at a slight nonlinear regime corresponding to a DC voltage V dc =5V. Analytical and numerical results are in good agreement. Particularly, at this nonlinear level, there is no notable difference between one, two or three modes for the HBM+ANM model. Indeed, the slight variations reported here between the resonance curves are below 15 parts per million (ppm) in resonance frequency and less than 1% in peak amplitude, which are very low compared to 0.5% in frequency due to residual stress [START_REF] Brueckner | Strainand pressure-dependent RF response of microelectromechanical resonators for sensing applications[END_REF] or micromachining imperfections as well as 5% in peak amplitude due to the variation of the quality factor by ohmic losses [START_REF] Sazonova | A tunable carbon nanotube resonator[END_REF] for a reasonable DC voltage range. Obviously, the coupling between the different modes is extremely low in this case and the influence of higher modes is negligible.

Then, the polarization voltage is increased up to 9V. Consequently, the dynamic behavior of the resonator becomes strongly nonlinear as shown in Fig. 6 for peak amplitudes between 20% and 35% of the gap. In these configurations, the coupling between the modes is strongly amplified.

Nevertheless, the error between the analytical model and the HBM+ANM model is still negligible, even with respect to the computational solution with 3 modes (frequency shift< 0.1%).

Finally, for confirmation, the same investigations have been made on a smaller resonator with a different geometry described in Table 1 (Design 2) . Figure 7 displays the confrontation of both models at a high nonlinear regime for a polarization voltage V dc =5V (V ac =0.5V). The same conclusions as before can be drawn concerning accuracy, which completes the analytical model validation.

The analytical model has the advantage to be fast and accurate. Furthermore, with this model it is possible to perform analytical parametric investigations with respect to the phase of the resonator oscillation and thus to derive analytical expressions that can be used as design rules by MEMS and NEMS designers in order to enhance the performances of resonant sensors.

Sensitivity to axial force -Acceleration sensing

Resonators can be used for acceleration sensing. Basically, as shown in Fig. 8, a resonant accelerometer consists in a proof mass M, suspended on the substrate by anchors (made of narrow beams), and a resonator placed close to the intersection of the anchors. The considered resonator is a clamped-clamped micro/nanobeam electrostatically actuated. An external acceleration γ that is applied to the proof mass along the sensitive X-axis of the device, results in a force communicated axially onto the resonator. The applied axial force induces a shift in the resonator frequency due to a change in the nominal stored potential energy of the system. Thus, by evaluating the frequency shift, the acceleration applied to the device can be estimated.

The sensitivity or scale factor of the device is the constant of proportionality relating the input (acceleration) and output (frequency shift). Since the relation between the acceleration γ and the communicated force to the resonator is linear, the scale factor in term of axial force can be used for parametric investigations with respect to the drive voltage when the resonator is driven in the nonlinear regime. By using the analytical model in its parametric form with respect to the phase β, the relative sensitivity of the resonator to an axial force can be written as:

S F = Δ f / f Δ Ñ0 = 1 f (β = π 2 , Ñ0 = 0) d f (β = π 2 , Ñ0 ) d Ñ0 (39)
with the forcing frequency f = Ω/2π, Ω computed from Eqs. ( 34)- [START_REF] Sazonova | A tunable carbon nanotube resonator[END_REF]. Figure 9 shows the variation of the resonator sensitivity (Design1) with respect to the drive DC voltage for two axial forces Ñ0 = 10 -6 N and Ñ0 = 10 -3 N. The difference between both curves is negligible which ensures a large dynamic range inside which the scale factor is linear and is similar to the first order Taylor series expansion of Eq. ( 37). However, remarkably, the scale factor increases with respect to the drive voltage from 4.82N -1 for V dc = 1V up to 5.1N -1 for V dc = 15V (W Max ≈ 0.95) which is due to the variation of the nonlinear stiffness with respect to applied axial force Ñ0 at large amplitudes. Indeed, the resonator sensitivity varies with respect to the level of nonlinearities as opposed to the linear case for which the scale factor of the resonant sensor is constant.

As expected and due to the scaling effects, the sensitivity of the second resonator is much higher than the first one which is shown in Fig. 10. Interestingly, the same investigations as before demonstrate the decrease of the scale factor with respect to the DC voltage from 750N -1

for V dc = 1V down to 630N -1 for V dc = 8V. This can be explained by the onset of a typical mixed hardening-softening behavior [START_REF] Kacem | Nonlinear dynamics of nanomechanical beam resonators: improving the performance of nems-based sensors[END_REF][START_REF] Kacem | Dynamic range enhancement of nonlinear nanomechanical resonant cantilevers for high sensitive nems gas/mass sensors applications[END_REF][START_REF] Kacem | Bifurcation topology tuning of a mixed behavior in nonlinear micromechanical resonators[END_REF] when the resonator is actuated using a DC voltage beyond 5V. 21 As a conclusion, for smaller ratio g/h the sensitivity of the resonator to axial force is improved but can be drastically reduced at large displacements i.e. high actuation voltages.

Conclusions

In this paper, the nonlinear dynamics of MEMS and NEMS clamped-clamped beam resonators electrostatically actuated was modeled including the main sources of nonlinearities as well as the fringing field effects. The modal Galerkin decomposition method was used in order to transform the nonlinear Euler-Bernoulli PDE that governs the resonator motion into a system of a coupled nonlinear ODE.

First, a multimodal numerical approach was developed by using the harmonic balance method (HBM) coupled with the so-called asymptotic numerical method (ANM), and validated with respect to a reference solution obtained by shooting.

Then, it was shown that the error made using only one mode is negligible compared to fabrication errors. A reduced order model neglecting the higher modes was therefore developed. To this end, the mode properties were investigated in order to simplify the matrices of the considered nonlinear system and to provide simplified analytical expressions for the Galerkin coefficients. Then, a perturbation method, namely the averaging method, was used to obtain closed-form expressions for the resonator dynamics which proved to be in very good agreement with the numerical results obtained with the HBM+ANM.

Moreover, parametric investigations can be easily performed. Analysis concerning the variation of the resonator sensitivity to axial force (in the case of a resonant accelerometer) demonstrates that when working in the nonlinear regime the scale factor is improved when the resonator is scaled down to the NEMS level, but it can be significantly reduced at high actuation voltages due to the soon onset of the mixed behavior. Remarkably, the variation of the scale factor with respect to the nonlinearity level was demonstrated. Consequently and in addition to the noise mixing issue [START_REF] Kaajakari | Phase noise in capacitively coupled micromechanical oscillators, Ultrasonics, Ferroelectrics and Frequency Control[END_REF],

this is another motivation for the nonlinearity cancellation strategy [START_REF] Kacem | Nonlinear dynamics of nanomechanical beam resonators: improving the performance of nems-based sensors[END_REF][START_REF] Kacem | Bifurcation topology tuning of a mixed behavior in nonlinear micromechanical resonators[END_REF] in resonant MEMS and NEMS sensors.

More generally, this reduced-order analytical model is a fast but accurate tool for the prediction of the nonlinear behavior of MEMS and NEMS resonators. It can be used to derive simple design rules in order to enhance the performances of resonant sensors, which is of high interest for MEMS designers. On the other hand, although it is computationally more demanding, the numerical model is a powerful tool which can be used even when the dynamic response of the resonator is more complex.

Appendix

Φ i Φ 3 k dx 1 π × ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 3(1+[-1] i+k ) 1+i+k + 1+[-1] i+3k 2+i+3k + 6(1+[-1] i+k )(i-k) 1+2i+2i 2 +2k+2k 2 - 3(1+[-1] i+k )(1+2k) 2+2i+i 2 +6k+2ik+5k 2 + 4-2i+10k+3(1-i+3k)[-1] i+k +(1+i+k)[-1] i+3k 2(-1+i-3k)(1+i+k) - 6(1+[-1] i+3k )(1+i+k) 5+6i+2i 2 +14k+8ik+10k 2 - 6(1+[-1] i+3k )(1+2k)(-1+4i+2i 2 -8k+4ik-10k 2 ) (1-2i+2i 2 +6k-8ik+10k 2 )(5+6i+2i 2 +14k+8ik+10k 2 ) + 2(1+[-1] i+3k )(-1+i-3k) 5+2i+2i 2 +18k+18k 2 - 2(1+[-1] i+k )(-1+5i+2i 3 -11k+20ik-6i 2 k-32k 2 +26ik 2 -30k 3 ) (1+2i+2i 2 +2k+2k
2 )(5+2i+2i 2 +18k+18k 2 ) -3(1+[-1] i+3k )(-1-i+i 3 -5k-4ik-3i 2 k-8k 2 -ik 2 -5k 3 ) (1+i 2 +4k-2ik+5k 2 )(2+2i+i 2 +6k+2ik+5k 2 )

⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ ∀(i, k) ∈ N * × N * and ∀i k 1 0 Φ i Φ 2 k Φ k dx π(1 + 2k) 2 4 × ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1+[-1] i+3k 2+i+3k - 2(1+[-1] i+k )(i-k) 1+2i+2i 2 +2k+2k 2 + 2(1+[-1] i+3k )(-1+i-3k) 5+2i+2i 2 +18k+18k 2 -(1+i+k)[-1] i+3k -3(-1+i-3k)[-1] i+k -2(-2+i-5k) 2(-1+i-3k)(1+i+k) -1+[-1] i+k 1+i+k + (1+[-1] i+k )(1+2k) 2+2i+i 2 +6k+2ik+5k 2 - 2(1+[-1] i+3k )(1+i+k) 5+6i+2i 2 +14k+8ik+10k 2 + 2(1+[-1] i+3k )(1+2k)(-1+4i+2i 2 -8k+4ik-10k 2 ) (1-2i+2i 2 +6k-8ik+10k 2 )(5+6i+2i 2 +14k+8ik+10k 2 ) + 2(1+[-1] i+k )(-1+5i+2i 3 -11k+20ik-6i 2 k-32k 2 +26ik 2 -30k 3 ) (1+2i+2i 2 +2k+2k 2 )(5+2i+2i 2 +18k+18k 2 ) - (1+[-1] i+3k )(-1-i+i 3 -5k-4ik-3i 2 k-8k 2 -ik 2 -5k 3 ) (1+i 2 +4k-2ik+5k 2 )(2+2i+i 2 +6k+2ik+5k 2 ) ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭

Figure 1 :

 1 Figure 1: Schema of an electrically actuated micro-beam

Figure 2 :

 2 Figure 2: HBM+ANM continuation on a strongly nonlinear behavior (Design 1, V dc =9V, V ac =0.9V), N m =1, H=7, N=30.

Figure 4 :

 4 Figure 4: Second confrontation Shooting/HBM+ANM on a strongly nonlinear behavior (Design 1, V dc =9V, V ac =0.9V).

Figure 5 :

 5 Figure 5: Confrontation HBM+ANM/Analytical model on a slightly nonlinear behavior (Design 1, V dc =5V, V ac =0.5V).

Figure 6 :

 6 Figure 6: Confrontation HBM+ANM/Analytical model on a strongly nonlinear behavior (Design 1, V dc =9V, V ac =0.9V).

Figure 7 :

 7 Figure 7: Confrontation HBM+ANM/Analytical model on a strongly nonlinear behavior (Design 2, V dc =5V, V ac =0.5V).

Figure 8 :

 8 Figure 8: Resonant accelerometer structure

Figure 9 :Figure 10 :

 910 Figure 9: Variation of the resonator scale factor with respect to the drive DC voltage for design 1 (V ac = 0.1V dc )

4 ∀i ∈ N * 1 0Φ

 41 + 1)[(2k + 1)π -4] 2i-1 dx 8 π [1 -4i]

1 0

 1

Preprint submitted to Non-Linear MechanicsDecember[START_REF] Baguet | On the behaviour of the anm in the presence of bifurcations[END_REF] 2010

∀(i, k) ∈ N * × N * and ∀i k