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.


Introduction
Coulomb-actuated microbeams are beams of microscopic dimensions deformed by the application of electrostatic forces. In addition to membranes, they are among the most important components used for capacitive actuation and sensing in microelectromechanical systems (MEMS) [1][2][3] . Recent advances in MEMS modeling and design enabled pioneering devices for medical applications (Bio-MEMS) 4 , communication systems (radiofrequency MEMS, RFMEMS) 5 , environmental sensing 2 , and consumer products 6 . Prominent examples are microresonators 7-9 , micropumps 10,11 , accelerometers 12,13 , gyroscopes 14,15 , capacitive micromachined ultrasonic transducers (CMUTs) 16,17 , microphones 18,19 , and microloudspeakers 20,21 . Global trends such as the 5G Internet of Things (5G-IoT) 22 , augmented reality 23 , and Green ICT (information and communications technology) 24,25 drive a system complexity that can only be handled by introducing lumped-parameter modeling at the MEMS level; an ab initio approach seems impractical. This requires a deep understanding of the electromechanics of the MEMS components involved.
Coulomb-actuated microbeams are commonly modeled as Euler-Bernoulli beams 7,26 in regard to lumped-parameter modeling. The models of Nayfeh, Younis, and Rahman [27][28][29][30][31][32][33] represent the current state-of-the-art in this field. These authors proposed a modal decomposition technique to analyze the static and dynamic behavior. Unfortunately, Younis et al. encountered convergence issues (for example, reported in Fig. 3 in ref. 29 ). At least three symmetrical modes are needed in their approach to reach decent accuracy. More modes do not necessarily mean more accuracy. Furthermore, all of this also comes at high computational complexity compared to numerical solvers 29 .
Classical mechanics suggests to the contrary that there hardly is a need to include higher spatial modes into lumped-parameter models for the type of beams and load situations relevant to MEMS components. The deflection profile of an Euler-Bernoulli beam, clamped at both ends and subject to a symmetric unidirectional load, is almost independent of the spatial distribution of load magnitude 34 . Within very small error margins, it resembles the shape of the lowest eigenfunction (zero mode) of the simple Euler-Bernoulli eigenvalue equation 33,35,36 d 4 dξ 4 ψ n ξ ð Þ ¼ λ n ψ n ξ ð Þ ð1Þ where ψ n (ξ) is the nth eigenfunction, λ n is the nth eigenvalue, and ξ is the coordinate along the beam. This fact is illustrated in Fig. 1a, comparing the Euler-Bernoulli zero-mode with the bending profiles of the two limiting load cases: a constant distributed load and fully concentrated load (point load). The explanation for this phenomenon is simple as well as far reaching. The bending profile w(ξ) resulting from a load q(ξ) is the solution to the equation satisfying appropriate boundary conditions 33,34 . Upon decomposition of the load function with respect to a set of normalized Euler-Bernoulli eigenfunctions 33,36 , we find for the bending profile: In case the load does not change its sign, we always haveq 0 ≠ 0. The key observation is now that for an Euler-Bernoulli beam, we have largely separated eigenvalues: λ 1 λ 0 % 29 and λ n grow with n 4 (see Supplementary Section 1, Table S1, and Eq. (S3)). Therefore, the Euler-Bernoulli beam acts as a kind of "low-pass filter", allowing only the first spatial mode, the zero mode, to by far dominate the picture. Note that the eigenvalues are directly related to the elastic energy required to excite the respective eigenmode. The eigenvalue separation therefore also implies a dominant impact of the zero  Does this also hold true for Coulomb-actuated microbeams? We certainly have a symmetric load, pulling in only one direction.
If the zero mode in fact dominates the bending of Coulomb-actuated microbeams to very much the same degree as is the case in pure elastomechanics, this implies major simplifications with respect to efficient modeling of the actuator physics. This would mean, for instance, that only one spatial degree of freedom would be enough to simulate the harmonic distortion of a loudspeaker based on such actuators 21 .
The notion that a single mode of oscillation suffices to describe the physics of a Coulomb-actuated microbeam is not self-evident. This stands in sharp contrast with the findings of Younis et al. 29 . Moreover, the singular nature of the Coulomb force itself nurtures substantial doubt that the lesson taken from elastomechanics outlined above is of any use in the MEMS world. After all, the Coulomb force gives rise to the coexistence of stable and unstable equilibrium states, separated by bifurcation points. Does the saddle-node bifurcation, called pull-in 37,38 and marking the condition for a potential electromechanical collapse, resemble a kind of elastomechanical buckling? This would destroy all prospects to arrive at a powerful zero-mode-based lumped-parameter model that is amenable to direct physical interpretation. Recently, Spitz et al. observed experimentally that the dynamics of a fairly complex cantilever design can in fact be well described by a lumped-parameter model with only a single degree of freedom 39 . This motivated us to revisit the physics of Coulomb-actuated prismatic microbeams analytically, numerically, and experimentally.
The paper is organized as follows. We devise a finite element method (FEM) strategy that can handle the entire bifurcation diagram, including the pull-in instability and all unstable equilibrium states. Next, our FEM approach is validated by applying it to the reference case presented by Gilbert et al. 40 . We then turn to the main question of how much the Euler-Bernoulli zero-mode actually determines the Coulombactuated bending profile compared to the contributions of the three next higher modes. The difference between a pure mechanical load and Coulomb actuation, as well as the impact of the pull-in instability, is investigated by varying the beam thickness. In fact, we find that the physics of a Coulomb-actuated, clamped-clamped prismatic microbeam is by far dominated by a single mode, the Euler-Bernoulli zero-mode. Finally, we compare our numerical findings with those of a microbeam experiment.

Results
The Coulomb-actuated prismatic Euler-Bernoulli beam A prismatic Coulomb-actuated microbeam (see Fig. 1b) can be modeled according to the following equation 32 : where ξ, τ, c, and w(ξ) are dimensionless quantities, denoting the coordinate along the beam, the time, the damping, and the lateral displacement, respectively. In addition, α 1 and α 2 are geometry-dependent parameters, N is a dimensionless external tensile axial force, and v(τ) is the driving voltage. The dimensionless quantities are defined as follows: where x are the coordinates along the beam and l is the length of the beam in absolute terms; u being the displacement in the direction of the fixed electrode and g the electrode gap at zero deflection in absolute terms; t being the beam thickness in absolute terms; and where ∈ is the dielectric constant of the gap medium and E is Young's modulus of the beam 32 . For a clamped-clamped beam, we have the boundary conditions In the following, we consider the static deflection of a beam without the presence of an external axial stress. All time derivatives are set to zero, and the equation reduces to the following nonlinear nonhomogeneous equation of fourth order, containing γ, which is a nonlocal functional of the deflection w(ξ) For small deflections, the effect of the stress stiffening can be neglected, and this equation further simplifies to Note that upon reducing the drive voltage v to zero, a Coulomb-actuated beam, originally in a stable state, will return to zero deflection. This means that the effect of the Coulomb force in Eq. (14) approaches that of a constant load of magnitude α 2 v 2 . For further use, we note here that the same situation occurs with thicker beams beyond the scope of Euler-Bernoulli theory.
The load at the contact singularity A little bit less obvious than the constant load situation is the load situation reached when starting from an unstable equilibrium and reducing the drive voltage v to zero while keeping forces in equilibrium. The critical bending profile w c (ξ) of the equilibrium state eventually reaches (the upper bifurcation point) and touches the counter electrode at the beam center (contact singularity 29 ): This causes a strictly positive difference in the shear forces ΔS 0 at the clamping points In fact, the contact singularity is reached when the shear force develops a discontinuity at the beam center. This is shown in Supplementary Section 3. To evaluate the implied load situation, we observe and find upon using the boundary conditions Eq. (11) and Eq. (16) On the other hand, Eq. (12) together with Eq. (15) implies In other words, the nonnegative integrand in Eq. (19) has a nonzero Lebesgue measure and is according to Eq. (20) supported in only a single point ξ = 0. This is sufficient to conclude that it exists only in the sense of tempered distributions and that it is proportional to the Dirac delta function 36,41 . Therefore, the mechanical load at the upper bifurcation point is a fully concentrated load, The mechanical load experienced by a Coulombactuated Euler-Bernoulli beam therefore ranges between a constant load and a concentrated load at low voltages. However, the question remains whether the saddle-node bifurcation at the pull-in voltage changes the picture away from the template of the pure elastomechanic scenario.

Numerical analysis of the beam used by Gilbert et al.
To validate our numerical methods, we first reproduced the results of Gilbert et al. 40 Fig. 2a-c as cyan crosses. This includes two regions: (i) actuation on the stable branch between 0 and 0.3 μm and (ii) maximum displacement, which is constant over the drive voltage above a certain threshold. The latter occurs in the Gilbert case, when the movement of the beam center is stopped at u y = 0.6 μm by a spacer.
The comparison of the results obtained from applying the collocation solver and our two-dimensional (2D) FEM approach to Gilbert's geometry, without and with stress stiffening, is displayed in Fig. 2a, b. For the collocation method, the iteration scheme from Eq. (S14) is used; the initial values need to be set carefully to achieve convergence. For the 2D FEM model, the nonlinear geometric effects are activated in ANSYS accordingly. Both methods are capable of solving for unstable states beyond pull-in. It is remarkable that the collocation method can deliver convergent solutions for unstable equilibrium states under all these circumstances. However, the arclength solver achieves a solution significantly closer to the contact singularity than the collocation solver. Our numerical results are perfectly in line with Gilbert's findings on the stable branch.
Gilbert's method generally does not yield points on the unstable branch. The only exception is the leftmost equilibrium state on the upper branch, where a spacer, used in Gilbert's case, stops the beam. This point lies exactly on the unstable branch of the bifurcation diagram obtained by our FEM approach. Upon lowering the voltage below the threshold corresponding to this contact point, the flexible electrode is released from sticking (pullout). In principle, spacers can be used to explore more points on the unstable branch.
To analyze the stable contact branch, we extended our 2D FEM model to a three-dimensional (3D) model with nonlinear contact between the microbeam and a spacer, as used by Gilbert. Again, we apply the arc-length solver. The resulting deflection curve is shown in Fig. 2c. Gilbert's reference curve is now exactly reproduced in full detail. In addition, the complete unstable branch is obtained. This allows us to determine the pull-in and pullout points much more precisely compared to tracking stable states up to the voltage where a Newton-Raphson solver ceases to converge.

Modal decomposition as a function of thickness
The relative contributions of Euler-Bernoulli eigenmodes to the deflection profile were computed by means of the FEM and the arc-length method. This calculation was based on Eq. (S10). The impact of the pull-in was studied by a broad variation of the beam thickness of Gilbert's reference geometry 40 .
The voltage-deflection curves normalized by the pull-in voltage v PI and the electrode gap g are shown in Fig. 3a for different beam thicknesses t in the range between 0.12 and 2 μm. Apparently, the shape of the normalized bifurcation diagram varies significantly with the beam thickness t. Thinner beams are shown as orange to red curves. They are more affected by the stress-stiffening effect than thicker beams. Stress stiffening introduces an additional nonlinear restoring force that shifts the pull-in point to higher deflections and changes the shape of the curve. Note that thinner microbeams exhibit an almost linear region between v v PI ¼ 0:2 and v v PI ¼ 0:8. Such behavior is of interest in many MEMS applications, e.g., in microscale loudspeakers 21,42 .
Only even modes contribute to this symmetric load situation. The relative weights b 0 , b 2 , b 4 , and b 6 of the first four even modes according to Eq. (S10) are shown in Fig. 3b-e as a function of the normalized driving voltage. The Euler-Bernoulli zero-mode has by far the highest relative weight, amounting to at least b 0 ≈ 0.9974. This "worst" case emerges for the thinnest beam (t = 0.12 μm) close to the pull-in. The maximum absolute values for the second, fourth, and sixth modes are 0.068, 0.051, and 0.0078, respectively. All these maximums occur for thin beams, while the absolute values significantly drop when t increases.
Comparison of the FEM results with the analytical formula at the contact singularity Does the FEM solver provide the correct shape of the bending profile after tracking the bifurcation diagram dimensionless quantity α 1 (see Eqs. (6) and (9)). The thinnest FEM-simulated beam has a thickness t = 120 nm and a length l = 80 μm and behaves closely to a string when stress-stiffening effects are considered. This is reflected in Fig. 3f-i (red line), where the modal contribution of the upper bifurcation point (unstable branch) deviates from linear Euler-Bernoulli theory starting from α 1 = 10 2 . In the case of a lower bifurcation point on the stable branch ( Fig. 3f-i, blue line), the modal contribution deviates below α 1 = 10 0 , where a transition to Timoshenko theory occurs. The location of the break-off towards Timoshenko theory is not fixed on the α 1 -axis and is in fact determined by the parameter θ in Supplementary Section 3 in Eq. (S25). Although the Euler-Bernoulli assumptions are violated for low and high α 1 values, we observe in Fig. 3f-i that for both limiting cases in the huge range 10 −1 < α 1 < 10 3 (four orders of magnitude), the zero-mode contribution is far above 99%.

Experimental results
A scanning electron microscope image of the manufactured clamped-clamped microbeam is shown in Fig. 4a, and the optical experimental method for measuring the bending profile is illustrated in Fig. 4b-d. The experimental raw data were fitted to a linear superposition of the first four even Euler-Bernoulli eigenmodes. A subset of the raw data used and the respective fitted bending profiles are shown in Fig. 4e, while a smaller subset of microscopic images is illustrated in Fig. 4f. The coefficients were used to plot the deflection curve (beam center) in Fig. 4g and the modal contributions as functions of the normalized voltage. The fit statistics were used to compute the respective error bars.
Since single-crystal silicon used for manufacturing has unneglectable anisotropy, the 2D FEM model from the previous section was extended by introducing the anisotropic elasticity of crystalline silicon, as provided in the literature (see Eq. (S30) in Supplementary Section 5) 43 . To accurately match the experiment, the manufactured geometry was updated: t = 2.47 μm and g = 10.15 μm. A small compressive stress of 1.4 MPa was included in the 2D FEM model, which agrees with the experimental observations of Younis et al. 29 . Furthermore, the anchor geometry was included in the model to better approach the real boundary conditions of the microbeam. The modal contribution from experimental data and from FEM simulations (magenta solid and dashed lines) is displayed in Fig. 5. All experimental data are not shown here to avoid unnecessary complexity of the graphs.
The experimental bending profile in Fig. 4g (beam center) and the experimental modal contributions b 0 , b 2 , and b 4 in Fig. 5a-d are well matched by the FEM simulation. Furthermore, the experiment and 2D FEM simulation are in good agreement regarding a value for the pull-in voltage of approximately v PI = 115 V, the pull-in deflection around u y (x = 0) = 6.5 μm, and the global shape of the curve (see Fig. 4g). The slight deviations in the deflection could be the result of cumulative effects of local geometric imperfections and 3D-distributed compliance of the anchors, which are attached to the handle on the bottom surface only. Note that the error bars for the deflection curve in Fig. 4g and for b 0 in Fig. 5a are extremely small. The data for b 6 show substantially more scatter (see Fig. 5d), in line with significantly larger error bars. This indicates that tracking more than four modes is beyond the capability of this experiment. According to our FEM analysis, the nature of the anchor design causes small distortions of the Euler-Bernoulli boundary conditions of Eq. (11). Adding geometrical details of the anchor design results in, e.g., the pronounced minimum observed in the b 2 graph (see Fig. 6b). Figure 6 shows a cumulative plot of the considered modal contributions illustrating the tendency of decreasing contribution with increasing n. The experiment confirms the dominance of the Euler-Bernoulli zero-mode, despite minor deviations from a pure Euler-Bernoulli situation at the boundary conditions of a real MEMS microbeam.

Discussion
In Section "Numerical analysis of the beam used by Gilbert et al.", we demonstrated numerical results for a Coulomb-actuated microbeam with FEM combined with the arc-length method and with the collocation method. Comparison of our analysis with the results of Gilbert et al. 40 shows perfect agreement, while the FEM in combination with the arc-length solver delivers all bifurcation points and all branches in the 2D and 3D setups. To the best of our knowledge, none of the FEM studies on electrostatically actuated MEMS published so far considers all physically relevant bifurcation points 40,[44][45][46][47][48][49][50][51] . The challenge is that at the saddle-node bifurcation (called pull-in instability), the effective stiffness of the system is equal to zero 37,38 , and the widely used Newton-Raphson methods tend to fail accordingly 52,53 . From a computational point of view, this behavior is similar to pure structural instabilities, which can be handled by arclength solvers [54][55][56] . We open up this approach here to the evaluation of unstable states of Coulomb-actuated microbeams and other comparable systems. In fact, this method is also applicable to more complex beam designs, which exhibit the coexistence of a multiplicity of stable equilibrium states 42 .
In Section "Modal decomposition as a function of thickness", we proceeded to apply our FEM-based modal analysis to prismatic beams within a wide thickness range from 0.12 to 2.0 μm. Owing to the accessibility of all bifurcation points, it was possible to numerically observe the modal contribution at the stable and unstable     branches, including the saddle-node bifurcation (pull-in) and the state close to contact singularity (v → 0 and u y (0) → 1). In Section "Comparison of the FEM results with the analytical formula at the contact singularity", the existence of analytical solutions for the contact singularity and for the constant distributed load case is highlighted. These analytical solutions confirm the FEM results, and the domination of the zero mode is claimed to hold true for the entire applicable voltage range with a relative zeromode weight factor of more than 99% (zero-mode hypothesis). The experimental results reported in Section "Experimental results" substantially differ from the measurements of Coulomb-actuated MEMS presented in the literature 7,8,45,51,57 . The proposed and easily reproducible experimental technique enables accurate static analysis of the beam-bending profile realizable with a common microscope and a voltage supply system. This accuracy allows for the analysis of up to four Euler-Bernoulli modes, which shows excellent agreement with the FEM results on the stable branch. Anisotropic elasticity was considered for comparison of the results; however, in the case of the pure Euler-Bernoulli setting, the elastic properties degenerate to a single elasticity parameter, which is considered in the quantity α 2 and does not change the outcomes. Furthermore, the experimental results confirm that the zero-mode hypothesis holds even in the case of a mild but noticeable violation of the Euler-Bernoulli model assumptions. Figure 7 illustrates the boundaries of the realm of Euler-Bernoulli theory. Reducing the beam thickness eventually affects the transition from a Coulomb-actuated Euler-Bernoulli beam to a Coulomb-actuated string. The Coulomb-actuated string can be analyzed along the lines presented above. This analysis shows that the contact singularity of the Coulomb-actuated string again is the case of a fully concentrated load. The respective bending profile is triangular, as is well known from mechanics 58 . A low electric drive voltage at a low deflection amplitude causes a constant load as the limiting case. The respective bending profile is a parabola 58 . The general normalized bending profiles lie between these two limiting profiles, with the stable profiles being close to the parabola. This allows us to compute the upper and lower bounds of b 0 to be 0.9792 < b 0 < 0.9958 for all values of α 1 . The pink area in Fig. 7 lies within the bounds and illustrates the range covered by b 0 when v > 0. However, the FEM results in Fig. 3b show that b 0 falls below its value at the contact singularity only above α 1 = 6 × 10 1 at the verge of applicability of the Euler-Bernoulli theory. The value of the lower bound reflects the fact that the spectral distance of the eigenmodes of a string is much smaller than that of an Euler-Bernoulli beam growing with n 2 rather than with n 4 . Therefore, our argument, based on Eq. (5), presented in the introduction does not readily apply to a string. However, the pure string-like case is less of a challenge for modeling since the general bending profiles of the Coulomb-actuated string, including stress stiffening, can be computed analytically. However,  it may be of interest in the future to analyze the transitional region more closely.
In the case of thick beams, i.e., in the realm of Timoshenko theory, the smallest value of b 0 (most unfavorable situation) is obtained in the limit of low drive voltages (see Fig. 3b). These cases can be solved analytically and yield polynomials of up to third degree, depending on the analyzed situation. We obtain the b 0 boundaries of 0.9792 < b 0 < 0.9999. Noticeable deviations from the Euler-Bernoulli scenario are observed for θ values above θ = 5 × 10 −3 . Thick Timoshenko-like beams are less relevant to the design of electrostatic actuators since they require excessive drive voltages to work properly. In summary, the Euler-Bernoulli zero-mode approximation seems useful, beyond the strict limits of the Euler-Bernoulli theory.
The purpose of this paper is to demonstrate that the bending profile of a Coulomb-activated prismatic microbeam clamped at both ends is almost identical to the shape of the first Euler-Bernoulli eigenmode as defined in Eq. (1). To the best of our knowledge, this aspect is not considered in the state-of-the-art modeling of microbeams, as proposed by Nayfeh, Younis, and Rahman [27][28][29][30][31][32][33] . The issues encountered by Younis et al. 29 originate in their treatment of the Coulomb term in Eq. (6). By either multiplying the entire equation with the Coulomb denominator or alternatively by approximating the Coulomb term by a polynomial in w(ξ) before employing a Galerkin-type procedure, Younis et al. 29 change the singular nature of the differential equation. A nonperturbative method is required here to resolve the issue. In the sequel to this paper, we plan to propose such an approach based on techniques originating from analytical probability theory.

Conclusion
Stable and unstable states of prismatic Coulombactuated Euler-Bernoulli beams can be successfully simulated by combining the FEM with arc-length solvers. The resulting model predictions can be experimentally verified by combining direct optical observations with modal analysis. Both approaches confirm that the shape of the bending profile of a Coulomb-activated prismatic microbeam, clamped at both ends, is almost identical to the shape of the Euler-Bernoulli zero-mode. This is true for the entire applicable voltage range with a relative zeromode weight of more than 99%. This observation paves the way for lumped-parameter models of high accuracy that only use a single spatial degree of freedom, amenable to direct physical interpretation.
In a forthcoming paper 59 , we will present a suitable method to carefully treat the Coulomb singularity in Eq. (6), thereby resolving the apparent discrepancy between our findings and the literature mentioned above.
The fruits of all this work are an easy-to-handle and surprisingly accurate zero-mode approximation that has numerous applications in the modeling of the physics of electrostatic MEMS actuators.

FEM analysis using an arc-length solver
The deflection solution of the microbeam can be determined either by numerically solving Eq. (12) using, e.g., the collocation method (see Supplementary Section 2), or by utilizing a more general method called FEM. The latter is a common and very versatile approach to model micromechanical systems and avoids the assumptions made to achieve Eq. (12). Commercial FEM codes routinely handle certain nonlinear phenomena such as large deflections and stress stiffening. Typically, these codes efficiently solve the resulting large system of nonlinear equations employing the Newton-Raphson algorithm or one of its modifications. At the pull-in instability, the tangent matrix approaches an infinite slope, and the Newton-Raphson algorithms tend to fail. Therefore, the use of static FEM is usually confined to the analysis of stable equilibrium states, leaving unstable states inaccessible. Even the accurate determination of the pull-in voltage is cumbersome. In practical applications, often the last convergent voltage is accepted as a surrogate for the actual pull-in voltage.
Purely structural (elastomechanic) instabilities, such as buckling or snap-through, share the phenomenon of the tangent matrix approaching an infinite slope. For elastomechanics, Riks 54 developed in late 1970s a modification of the Newton method allowing for a continuation beyond these instabilities. Crisfield 55 improved Riks' method to allow for easy FEM implementation. Today, this method is referred to as the arc-length method. It is a standard method included in commercial codes such as ANSYS and ABAQUS 56 . In this paper, we use the arclength method to solve a coupled-field problem including electrostatics and elastomechanics. This allows us to compute the entire bifurcation diagram of a Coulombactuated prismatic Euler-Bernoulli beam, including all unstable states.

The experimental setup
To experimentally determine the beam deflection, a Coulomb-actuated MEMS microbeam was manufactured on a bonded silicon on insulator (BSOI) wafer using deep reactive ion etching (Bosch process); see Fig. 4a, b. The 200 mm BSOI wafer consisted of a 75 μm device wafer p-doped with boron at a concentration of 10 18 cm −3 , a 1 μm buried oxide layer (SiO 2 ), and a 650 μm handle wafer. The prismatic beam design had a beam length of l = 1000 μm, a beam width of w = 75 μm, and a beam thickness of t = 3 μm. The beam is shown in pink in Fig. 4b and was designed for in-plane movement upon activation by a planar counter electrode (violet) placed at a distance of g = 10 μm. The beam was clamped in two prismatic anchors with square cross-sections (50 μm × 67.7 μm on the left, 237 μm × 228 μm on the right) and isolated from the rest of the chip by a trench of width 3 μm. To allow for beam movement, the buried oxide layer was removed underneath the beam by means of HF (hydrogen fluoride) vapor etching, and the beam was operated under ambient atmospheric conditions. The actual geometry after fabrication was determined by microscopy, revealing an actual average thickness of t = 2.47 μm and an electrostatic gap of g = 10.15 μm. To quantify the oxide underetching underneath the anchors of estimated 13 μm, beams were detached from the MEMS chip and inspected from the bottom. A resulting scallop pattern height of approximately 200 nm was determined using scanning electron microscopy. Since the scallop grooves are oriented along the beam axis and their size is small compared to the microbeam dimensions, it was concluded that their impact essentially is a contribution to the average beam thickness t, included in the value given above.
The bending profile of the microbeam was recorded using a Leica DM8000 M microscope with an objective 20×/0,4 N PLAN L. All pictures were captured using Bandicam 4.5.6 1647 software in MP4 format at a rate of 30 frames per second. Since the objective could not map the entire beam into one picture frame, only the deflection of one half of the microbeam was recorded (see visible area in Fig. 4b). The location of the beam center was marked on the fixed electrode during manufacturing. A driving voltage sweep from 0 to 120 V was supplied by a Keysight B2902A source unit, set to a slew rate of 1 V s −1 , small enough to avoid dynamic effects. The driving voltage was applied to the MEMS microbeam contact pads by two Imina MiBot microprobe stations equipped with ST-20-0.5 wolfram needles.
The deflection data were extracted directly from the recorded images, and a typical image is shown in Fig. 4c. A total of 3685 frames were recorded, out of which 34 were used. Per frame, 866 cross-sections were evaluated. At every position x along the beam, the normalized color channel data for the respective section across the beam were used to determine the beam deflection at this position. The red line in Fig. 4c marks a typical cross-section used. The corresponding color channel data are shown in Fig. 4d. To extract the geometric positions, the heuristic function L(y), L y ð Þ ¼ G y À y G ð ÞþH y H À y ð Þ; was fitted to the color channel data. The coordinate y is the coordinate along the beam cross-section. G(y) is the and H(y) is the Heaviside step function, The fit parameters were σ, the position y G of the middle of the microbeam along the respective cross-section and the position y H of the relevant edge of the fixed electrode. The respective fit is also shown in Fig. 4d. This finally leads to the raw data for the experimental deflection profile