Coulomb-actuated microbeams revisited: experimental and numerical modal decomposition of the saddle-node bifurcation

Electrostatic micromechanical actuators have numerous applications in science and technology. In many applications, they are operated in a narrow frequency range close to resonance and at a drive voltage of low variation. Recently, new applications, such as microelectromechanical systems (MEMS) microspeakers (µSpeakers), have emerged that require operation over a wide frequency and dynamic range. Simulating the dynamic performance under such circumstances is still highly cumbersome. State-of-the-art finite element analysis struggles with pull-in instability and does not deliver the necessary information about unstable equilibrium states accordingly. Convincing lumped-parameter models amenable to direct physical interpretation are missing. This inhibits the indispensable in-depth analysis of the dynamic stability of such systems. In this paper, we take a major step towards mending the situation. By combining the finite element method (FEM) with an arc-length solver, we obtain the full bifurcation diagram for electrostatic actuators based on prismatic Euler-Bernoulli beams. A subsequent modal analysis then shows that within very narrow error margins, it is exclusively the lowest Euler-Bernoulli eigenmode that dominates the beam physics over the entire relevant drive voltage range. An experiment directly recording the deflection profile of a MEMS microbeam is performed and confirms the numerical findings with astonishing precision. This enables modeling the system using a single spatial degree of freedom.


SUPPL. 1. SYMMETRIC EULER-BERNOULLI EIGENMODES
For a prismatic beam, clamped at both ends, the Euler-Bernoulli eigenmodes are defined by the eigenvalue equation Eq. (1) subject to the boundary conditions Eq. (11). The eigensystem for even indices, λ 2n and ψ 2n (ξ), are given by , (S1) where the β 2n are the zeros of and can be approximated using β 2n = 2n + 3 2 π+ (S4) Table S1 illustrates numerically determined β 2n from Eq. (S3) in comparison to the approximation using only the first term in Eq. (S4). It is readily verified that the function in Eq. (S1) satisfies Eq. (1) and the boundary conditions in Eq. (11). Note that the eigenmodes as defined above are ortho-normal, Moreover the even eigenmodes form a complete ortho-nomal base of the Hilbert space of symmetric square integrable functions over the interval − 1 2 , + 1 2 . This means that we can expand any symmetric deflection profile w(ξ) in terms of these eigenmodes where we have due to ortho-normalitŷ Note that from these definitions we also get Parseval's equation To asses the relative contribution of the individual eigenmodes to the deflection profile Finally, to be able to apply the above formula for the evaluation of the experimental data, we use Parseval's equation (S10) gorithm. We decided to use a fourth order collocation algorithm provided by the SciPy library 1,2 , which requires a system of first order differential equations as input. To this end we cast Eq. (12) into a set of four ordinary equations, containing however the integral γ (S11) In order to deal with γ, we introduce a new function The function Γ(ξ) can be directly obtained during the solution procedure by adding the differential equation to the system in Eq. (S11). We then define a simple iteration scheme starting with γ 0 = 0 and solving the new system according to where i is the iteration index and γ i is updated using The iteration is terminated upon meeting the criterion where γ is the targeted accuracy. Depending on parameter settings, a damping constant c ID < 1, adjusting the step size, may be required to reach convergence

SUPPL. 3. EULER-BERNOULLI BEAM SUBJECT TO A CONCENTRATED LOAD
In this section we compile the solution to the concentrated load equation (Eq. (21)) including stress-stiffening γ according to Eq. (13) and subject to the boundary conditions Eq. (11). We start by noticing that outside the beam center we need to solve Eq. (20) It is worth emphasising that Eq. (20) holds, irrespective of our claims about the shear force.
These can therefore be verified by inspecting the results below. The respective bending profile, compatible with the boundary conditions in Eq. (11) and satisfying Eq. (15) is where the the parameters c a , c b and c c are given by (S19) This may be readily verified by direct computation upon inserting Eq. (S18) into Eq. (20) and into Eq. (15).
The bending profile w c (ξ), its first and second derivative are continuous at the beam center ξ = 0. However the third derivative w This establishes that w c (ξ) as defined in Eq. (S18) is the bending profile resulting from a concentrated load at the beam center, i.e. it is the solution to Eq. (21). In the limit of a vanishing stress stiffening we obtain To establish the relation between γ and α 1 , required to analytically deal with stressstiffening at the upper bifurcation point, we use Eq. (13). Performing the necessary integration and solving for α 1 yields For practical applications it is helpful noticing Finally we compute the contributions of the Euler-Bernoulli eigenmodes to w c (ξ) according to Eq. (S10), which admittedly is a bit tedious

SUPPL. 4. TIMOSHENKO BEAM SUBJECT TO A CONSTANT LOAD
At very low drive voltages the Coulomb force generates a constant beam load of magnitude α 1 v 2 . In this section we will compile the relevant formulae. Actually we will do this for the more general case of a Timoshenko beam. This helps us exploring the limits of Euler-Bernoulli theory when varying the beam thickness. In contrast to the Euler-Bernoulli assumptions, Timoshenko beam theory allows for a rotation of the normal to the mid-surface of the beam. The rotation angle is denoted here by φ(ξ). In the formulae below, w(ξ) is as usual, the displacement in dimensionless form. The Timoshenko beam equations for a prismatic beam of rectangular cross section, subject to a constant unit load then read 3 where α 1 is taken from Eq. (9) and the value for ν zx is given in Eq.