Buckling Electrothermal NEMS Actuators: Analytic Design for Very Slender Beams Buckling Electrothermal NEMS Actuators: Analytic Design for Very Slender Beams

: Analytic approximations are presented for the response of buckling-mode electrothermal actuators with very slender beams with a width-to-length ratio of W / L ≤ 0.001 of the type found in nanoelectromechanical systems (NEMS). The results are found as closed-form solutions to the Euler beam bending theory rather than by an iterative numerical solution or a time-consuming ﬁnite element analysis. Expressions for transverse deﬂections and stiffness are presented for actuators with the common raised cosine and chevron pre-buckled shapes. The approximations are valid when the effects of bending dominate over those of axial compression. A few higher-order approximations are also presented for less slender beams with 0.001 ≤ W / L ≤ 0.01. Abstract: Analytic approximations are presented for the response of buckling-mode electrothermal actuators with very slender beams with a width-to-length ratio of 𝑊/𝐿 ≤ 0.001 of the type found in nanoelectromechanical systems (NEMS). The results are found as closed-form solutions to the Euler beam bending theory rather than by an iterative numerical solution or a time-consuming finite element analysis. Expressions for transverse deflections and stiffness are presented for actuators with the common raised cosine and chevron pre-buckled shapes. The approximations are valid when the effects of bending dominate over those of axial compression. A few higher-order approximations are also presented for less slender beams with 0.001 ≤ 𝑊/𝐿 ≤ 0.01 .


Introduction
Electrothermal actuators are widely used in microelectromechanical systems (MEMS) because a device with in-plane motion may be constructed simply by etching and undercutting a mechanical layer to form a suspended structure that is heated by passing a current between the anchors [1]. The most common arrangements are shape bimorph [2,3] and buckling [4][5][6] actuators. Here, we focus on the latter, which have the general arrangement as shown in Figure 1a. Here, an array of beams is supported on anchors at either end. The beams are pre-buckled to force motion in the direction shown, tied together with a crossbeam at their midpoint to ensure collective motion, and driven by constrained thermal expansion.

Introduction
Electrothermal actuators are widely used in microelectromechanical systems (MEMS) because a device with in-plane motion may be constructed simply by etching and undercutting a mechanical layer to form a suspended structure that is heated by passing a current between the anchors [1]. The most common arrangements are shape bimorph [2,3] and buckling [4][5][6] actuators. Here, we focus on the latter, which have the general arrangement as shown in Figure 1a. Here, an array of beams is supported on anchors at either end. The beams are pre-buckled to force motion in the direction shown, tied together with a crossbeam at their midpoint to ensure collective motion, and driven by constrained thermal expansion. If the beams are initially shaped into a chevron layout, the device is known as a V-beam actuator. The basis of actuation has been investigated using the Euler buckling theory and a finite element analysis [7,8], typically using the simplified single-beam model in Figure 1b. More recently, geometry effects [9], non-linearity [10], and dynamics [11,12] have been investigated, and there have been continuing attempts to improve the FEM analysis [13][14][15][16]. An alternative layout is the raised cosine pre-buckle in Figure 1c, which has been extensively analysed under the conditions leading to snap-through and bistability [17][18][19]. The thermal aspects of buckling actuators have also been modelled [20,21], and a new designthe Z-shaped actuator, which uses a stepped beam layout-has been developed [22][23][24]. More complex arrangements have been developed, including cascaded [7,[25][26][27][28][29] and outof-plane [30] actuators as well as actuators with a variable beam width [31,32], built-in strain gauges [33], and feedback control [34,35]. V-beam actuators have been used as the basis of linear [36][37][38] and rotary [36,39] stepping motors. Applications include electrical switches [40,41], optical alignment systems [42][43][44][45], tunable or scanning micro-optical devices [46][47][48], devices for cell manipulation [49][50][51], and movable neural microelectrode arrays [52,53]. The devices and applications are reviewed in [54]. Despite approximately twenty years of research, a few difficulties remain. Neither the analytic nor the finite element models of elastic deformation are very satisfactory. The former does not yield closed-form expressions and the latter requires excessive computation when the beams are very long and narrow and the number of elements is large. Apart from an analytic approximation for the deflection in [6], derived from purely geometric arguments, there are no closed-form design formulas so layouts are often chosen on an ad hoc basis. The aim of this paper is to remedy this deficiency by developing approximations to buckling theory that provide closed-form solutions. The solutions are valid for very slender beams with a width-to-length ratio in the approximate range of W/L ≤ 0.001 when bending dominates over axial compression; nanoelectromechanical systems (NEMS) actuators of this type capable of micron-scale deflections have been demonstrated by sidewall transfer lithography [55]. Raised cosine actuators and actuators with a chevron pre-buckle are analysed and discussed in Section 2. Additional results including alternative approximations, intrinsic stress, the relation between deflection and power, and transverse stiffness are considered in Section 3. Conclusions are drawn in Section 4.

Methods
We started by considering the raised cosine pre-buckle but followed the general approach used for V-beam actuators in [7]. We assumed a single beam with the shape shown in Figure 1c; the performance of the beam arrays could be estimated by scaling. The beam had a length L and central offset H, a rectangular cross-section with an in-plane width W and depth D, and was built-in at both ends. Its initial shape was: Forces F and moments M A modelled the effects of the anchors and the shape y(x) of the loaded beam could be found by solving the Euler buckling equation [56]: Here, M(x) = M A − Fy is the bending moment, I = W 3 D/12 is the second moment of the area of the beam, and E is the Young's modulus for the material used. An approximation implicit in the above was the replacement of the change in the beam curvature by the second derivative, which could clearly render the result inaccurate for a large pre-buckle or deflection. The analysis was best carried out by substituting y = y 0 + u, which yielded: Micro 2022, 2 56 Equation (3) may be solved in terms of a particular integral and a complementary function. Applying the boundary conditions u = du dx = 0 at x = 0 and x = L, the result was: Here, k = √ F/EI. This result demonstrated a special quality of the raised cosine pre-buckle: the deflected shape did not change. The midpoint deflection U = u(L/2) was: Equation (5) implies that U would be infinite when kL = 2π, the Euler buckling condition. Integrating the deflected beam shape allowed the change in length ∆L due to bending to be found as ∆L = δL − δL 0 , where we made the standard approximation δL = 1/2 L 0 dy dx 2 dx; δL 0 was the integral in the unloaded case (see, e.g., [9]). Integration gave: The force F was derived from the constrained thermal expansion. Assuming an average temperature rise ∆T avg , the beam must satisfy a compatibility condition [7]: Here, α is the linear thermal expansion coefficient of the beam material. The three terms in Equation (7) describe the changes in the length due to the axial compression, thermal expansion, and bending, respectively. This equation can also be written in the form: Equations (6) and (8) allowed ∆T avg to be found as a function of kL. As U is also a function of kL, it could be plotted as a function of ∆T avg . Alternatively, these coupled equations could be solved iteratively. For example, Figure 2 shows the variation of U with ∆T avg for a silicon beam with the following parameters: E = 170 × 10 9 Nm −2 , α = 2.6 × 10 −6 K −1 , L = 1 mm, H = 5 µm, D = 5 µm, and values of W ranging from 20 µm (MEMS domain) to 0.1 µm (NEMS). In each case, the variation was similar. When the temperature rise ∆T avg was negative (so the beam was artificially cooled and under tension), U was negative and tended to −H, straightening the beam. When ∆T avg was positive (so the beam was under compression), the deflection rose monotonically as the beam buckled. For values of W > 1 µm, the variations differed; however, for W ≤ 1 µm, the analytic solutions all tended to one another. Note, however, that for very large deflections, the approximation of the curvature changes by a second derivative in Equation (2) would no longer be valid.
Similar calculations were carried out using the commercial FEA solver, COMSOL ® Multiphysics 4.4. A 2D model of a beam with a raised cosine pre-buckle and the dimensions above was built in the x-y plane. The heat transfer module was used to apply a temperature change and the solid mechanics module to apply boundary constraints. Structural and thermal physics were coupled by thermal expansion in the multiphysics modelling interface. The plane stress approximation was used to eliminate the out-of-plane components of the stress tensor, and a stationary study was carried out using a geometrically non-linear analysis. The number of iterations was set to 50, with the solver automatically determining the damping factor for Newton's method. In addition to the previous parameters, a Poisson's ratio of ν = 0.28 was assumed. The discrete points in Figure 2 show the results for the four largest beam widths. The results were in good agreement with the analytic model for W = 20 µm and W = 10 µm but diverged for W = 5 µm and W = 1 µm; for smaller values of W, convergence was not reached after 50 iterations.
with ∆ for a silicon beam with the following parameters: = 170 × 10 Nm , = 2.6 × 10 K , = 1 mm, = 5 μm, = 5 μm, and values of ranging from 20 μm (MEMS domain) to 0.1 μm (NEMS). In each case, the variation was similar. When the temperature rise ∆ was negative (so the beam was artificially cooled and under tension), was negative and tended to − , straightening the beam. When ∆ was positive (so the beam was under compression), the deflection rose monotonically as the beam buckled. For values of 1 μm, the variations differed; however, for ≤ 1 μm, the analytic solutions all tended to one another. Note, however, that for very large deflections, the approximation of the curvature changes by a second derivative in Equation (2) would no longer be valid. A further difficulty was the increase in memory and run time as the W/L decreased. Table 1 shows the requirements for 37 temperature values using an Intel ® Core TM i7-2600 processor with a 3.4 GHz clock speed and 16 Gb RAM. The increase in the computation resource implied that the problem would rapidly become intractable, and these problems would worsen in 3D. Although complete, the analytic solution above was still unsatisfactory because it required a numerical evaluation. A closed-form solution for U in terms of ∆T avg could allow the performance to be estimated directly from the geometric parameters. We obtained such a solution by first comparing Equations (5) and (6). It is simple to show that: For slender beams with W/L 1, we neglected the second term on the right-hand side of Equation (8), effectively neglecting the effects of axial compression. In this case, ∆L/L = α∆T avg , and Equation (9) approximated to the quadratic: Equation (10) had two solutions. Retaining only the positive solution, the following expression for the deflection could be obtained: Equation (11) was the desired closed-form solution. It yielded a real result even when ∆T avg was negative until ∆T avg = −π 2 H 2 /4αL 2 . Figure 3 compares the exact (axial stress included) and approximate (axial stress ignored) variation of the deflection with the average temperature rise for the same parameters as Figure 2 but with only the smallest beam width (W = 0.1 µm, for which W/L = 10 −4 ). The agreement was excellent and the closed-form analytic formula tracked the exact solution closely, even when the deflection was large. Other curves are superimposed on this figure; these are discussed later.
Micro 2022, 2, FOR PEER REVIEW 5 closed-form analytic formula tracked the exact solution closely, even when the deflection was large. Other curves are superimposed on this figure; these are discussed later.
This result implied that the beam changed shape significantly as it deflected, in contrast to the raised cosine actuator. The maximum deflection was = ( − )| / or: Integrating the deflected beam shape again allowed the change in the length ∆ due to bending to be found. The result was [7]: Here, = tan . The compatibility equation was as before, and Equations (14) and (8) then allowed ∆ to be found as a function of . As Equation (13) was also a function of , could again be plotted in terms of ∆ . The results are superimposed on Figure 3 for the same parameters as before. The deflections were very similar to the variations obtained for the raised cosine actuator, suggesting that the deflection in the V-beam actuator was dominated by the excitation of the lowest-order buckling-mode. The results obtained using COMSOL (not shown) again agreed well but suffered from an even worse scaling of run times with / .
The use of the very slender beam approach-relating to ∆ / and then solving an approximate compatibility equation-could yield again an analytic approximation. This time, there were difficulties caused by the trigonometric functions in Equations (13) and (14). However, these could be circumvented using power series approximations. For example, substituting = /4 allowed Equation (13)   We repeated the analysis for the V-beam actuator shown in Figure 1b. The beam had a length L and tilt angle θ so the central offset was H = L 2 tan(θ). As y 0 is a linear function of x for 0 ≤ x ≤ L/2, Equation (2) could be written as d 2 y dx 2 + F EI y = M A EI in this range. The boundary conditions were y = 0 at x = 0 and dy dx = tan(θ) at x = 0 and x = L/2. The solution was [7]: This result implied that the beam changed shape significantly as it deflected, in contrast to the raised cosine actuator. The maximum deflection was U = (y − y 0 )| L/2 or: Integrating the deflected beam shape again allowed the change in the length ∆L due to bending to be found. The result was [7]: Here, G = tan kL 4 . The compatibility equation was as before, and Equations (14) and (8) then allowed ∆T avg to be found as a function of kL. As Equation (13) was also a function of kL, U could again be plotted in terms of ∆T avg . The results are superimposed on Figure 3 for the same parameters as before. The deflections were very similar to the variations obtained for the raised cosine actuator, suggesting that the deflection in the V-beam actuator was dominated by the excitation of the lowest-order buckling-mode. The results obtained using COMSOL (not shown) again agreed well but suffered from an even worse scaling of run times with W/L. The use of the very slender beam approach-relating U to ∆L/L and then solving an approximate compatibility equation-could yield again an analytic approximation. This time, there were difficulties caused by the trigonometric functions in Equations (13) and (14). However, these could be circumvented using power series approximations. For example, substituting λ = kL/4 allowed Equation (13) to be written as U L = tan(θ) f (λ), where: As tan(x) = x + x 3 3 + 2x 5 15 . . ., we could approximate f to order λ 4 as: In a similar way, Equation (14) could be written as ∆L L = tan 2 (θ)g(λ), where: To simplify the above, we first noted that: g(λ) could then be written as: With g in this form, we obtained a power series approximation to the order λ 4 as: These simple approximations clearly had limited validity. Figure 4 compares the variations of f , f 4 , g, and g 4 with λ; both approximations were accurate only up to λ ≈ 0.3 π, well below the buckling condition (when f and g both became infinite).
As tan( ) = + + …, we could approximate to order as: In a similar way, Equation (14) could be written as ∆ = tan ( ) ( ), where: To simplify the above, we first noted that:

tan( ) 1 − cos(4 )] + 1 − tan ( )]sin(4 ) = 4 tan( ).
( ) could then be written as: With in this form, we obtained a power series approximation to the order as: These simple approximations clearly had limited validity. Figure 4 compares the variations of , , , and with ; both approximations were accurate only up to ≈ 0.3 , well below the buckling condition (when and both became infinite). Despite this, they allowed useful closed-form solutions. For example, combining Equations (16) and (20), we obtained: Rearranging, we obtained the quadratic equation: Retaining only the positive root, we obtained:  Despite this, they allowed useful closed-form solutions. For example, combining Equations (16) and (20), we obtained: Micro 2022, 2

60
Rearranging, we obtained the quadratic equation: Retaining only the positive root, we obtained: Noting the original definitions of f and g, we then obtained the midpoint deflection U as: Neglecting the second term in the RHS of Equation (8) as before, we then obtained: This result can, of course, be written in the alternative form: Equation (26) was analogous to Equation (11) for the raised cosine pre-buckle. The predictions of this closed-form solution are superimposed on Figure 3 for the same parameters. Once again, there was a good agreement with the full model despite the apparent errors in f 4 and g 4 . We constructed higher-order expansions of both functions and found that the expressions above gave good results despite their simplicity.
The performances of the two designs were clearly similar; however, the V-beam actuator gave a slightly larger deflection for a given temperature rise. Equations (11) and (26) could both be written in the form U = A 1 + B∆T avg − 1 . Consequently, the initial sensitivity S = dU/d∆T avg 0 of the deflection to temperature was S = AB/2. For the two actuator types, the sensitivities S RC and S V were: Thus, the sensitivity was always proportional to L 2 and inversely proportional to H. For a common geometry, S V /S RC = π 2 /8, so the V-beam design was slightly more sensitive.

Results
In this section, we discuss several alternative results, power scaling laws and the related question of transverse stiffness. We started by considering the well-known formula for the deflection of a V-beam actuator derived in [6] by considering that the two beam sections remained straight and simply increased in length by an amount ∆L h due to thermal expansion. Using a first-order approximation we obtained: Here, L h = L/2. Rearranging and substituting H for L h sin(θ) and α∆T avg for ∆L h /L h we obtained: Equation (28) may also be written in the form U = A 1 + B∆T avg − 1 . Its predictions are also shown superimposed on Figure 3. It was clearly less accurate than Equation (26); however, the prediction was remarkable given that bending was entirely ignored.
We then considered alternative approximations for a raised cosine actuator that could be viable for less slender beams when the axial compression term in Equation (8) could not be neglected. Retaining this term, it was simple to show that Equation (10) becomes: Thus, a more exact relation between U and ∆T avg was a cubic equation. Clearly, Equation (30) had a numerical solution. However, in the spirit of the present MS, where the focus is on solutions that yielded an insight, we derived an analytic approximation. We assumed that U = U 0 + U 1 , where U 0 was the solution to Equation (10) and U 1 was an additional perturbation. Substituting and neglecting the higher-order terms, it was simple to obtain: Thus, the effect of axial compression was to reduce the deflection by an amount that depended on W and (for small deflections) was linearly proportional to U 0 . Consequently, the initial slope of the deflection characteristic must reduce for wider beams. To illustrate this, Figure 5 compares the exact solution, the first approximation (11), and the second approximation (31). The parameters were as before, but the beam width W was increased to 5 µm (so that W/L = 5 × 10 −3 ). This width was sufficiently large that the first approximation was no longer valid. However, the second approximation was a good match to the exact solution except when U 0 tended to −H. This behaviour was to be expected from Equation (31); it is unimportant for practical applications when positive temperature rises are the norm.
Thus, the effect of axial compression was to reduce the deflection by an amount that depended on and (for small deflections) was linearly proportional to . Consequently, the initial slope of the deflection characteristic must reduce for wider beams. To illustrate this, Figure 5 compares the exact solution, the first approximation (11), and the second approximation (31). The parameters were as before, but the beam width was increased to 5 μm (so that = 5 × 10 ⁄ ). This width was sufficiently large that the first approximation was no longer valid. However, the second approximation was a good match to the exact solution except when tended to − . This behaviour was to be expected from Equation (31); it is unimportant for practical applications when positive temperature rises are the norm. We noted that the addition of a tensile stress (which could arise during processing) inserted an axial strain term − / into Equation (8). The intrinsic stress could, therefore, be included by replacing ∆ with ∆ − / in the preceding expressions. The effect of stress was, therefore, to shift the deflection characteristics to the right. If measurements were made from the rest, the effect was to reduce the apparent deflection. Finally, we considered the well-known relation between temperature and power (see, e.g., [9,10]). If the beam was heated by a power and cooled only by thermal conduction to the anchors, the steady-state temperature rise ∆ ( ) was the solution to the heat conduction equation: We noted that the addition of a tensile stress σ (which could arise during processing) inserted an axial strain term −σ/E into Equation (8). The intrinsic stress could, therefore, be included by replacing ∆T avg with ∆T avg − σ/E in the preceding expressions. The effect of stress was, therefore, to shift the deflection characteristics to the right. If measurements were made from the rest, the effect was to reduce the apparent deflection.
Finally, we considered the well-known relation between temperature and power (see, e.g., [9,10]). If the beam was heated by a power P and cooled only by thermal conduction to the anchors, the steady-state temperature rise ∆T(x) was the solution to the heat conduction equation: Here, k th was the thermal conductivity of the beam material. The solution subject to the boundary conditions ∆T = 0 at x = 0 and x = L was: Here, ∆T max was the temperature at x = L/2. The average temperature could then be found by the integration of the temperature profile along the beam length as: The previous expressions for deflection could be converted into variations with power by making this substitution; for an N-beam array, P was simply divided by N. Similarly, for a device that was also cooled by gas conduction to the substrate, P was multiplied by an efficiency term, η.
We then considered the transverse stiffness of a buckling actuator and showed how the slender beam approximation could be used to evaluate this quantity for the raised cosine type. Following previous authors [7], we assumed that the stiffness was found from the transverse force F T needed to return the midpoint deflection to zero. Unfortunately, the application of F T altered the axial force to a new value, F , and altered the end moment to M A . At a certain point, the rising transverse force will lead to snap-through [9], so the analysis was carried out with caution. To find the deflection u (x), we solved: The variation of u was symmetric about the midpoint. For 0 ≤ x ≤ L/2, it could be found by standard methods as: Here: The midpoint deflection U = u (L/2) was then: Now, to return U to zero, we required a transverse force such that: As k was unknown, a further condition was required to determine F T . This could be obtained from the compatibility. As the temperature was unchanged by the transverse force, we had ∆L /L = ∆L/L, where ∆L /L = 1/L L/2 0 dy dx 2 dx. To proceed, we started with the deflected shape: Here, C 2 = H 4π 2 4π 2 −k 2 L 2 . Consequently: Using the above, the integral determining ∆L /L could be evaluated, numerically or analytically. However, the expressions obtained in the latter case were non-linear in k . In each case, the values of F T and F were, therefore, generally found by iteration to satisfy the conditions above for a given F. Figure 6 shows a typical result, where the initial shape y 0 (x), the actuated shape y(x), and the centrally loaded shape y (x) were plotted together for kL = π. The overall deflection increased as the actuator was driven, retaining a raised cosine shape; the midpoint deflection correctly returned to zero as the actuator was loaded.
A simple analytic solution for the transverse stiffness could be obtained by assuming that ′ was small and by making power series approximations to Equation (40). This process yielded: Squaring and integrating, we then obtained: At a similar level of approximation, we had: A substitution into the compatibility condition then yielded an approximation for ′ as:  Figure 6. Initial, deflected, and centrally loaded beam shapes y 0 , y, and y of a raised cosine actuator for kL/π = 1.
A simple analytic solution for the transverse stiffness could be obtained by assuming that k L was small and by making power series approximations to Equation (40). This process yielded: Squaring and integrating, we then obtained: Micro 2022, 2

64
At a similar level of approximation, we had: A substitution into the compatibility condition then yielded an approximation for k as: Thus, the action of forcing the actuator back to its starting point increased k significantly from the original value k due to electrothermal actuation alone. To find the transverse stiffness, we noted that F T and U could be approximated as: Consequently, the transverse stiffness k T = F T /U could be extracted as: The term k T0 = 192EI/L 3 was then the transverse stiffness of a straight, centrally loaded built-in beam. Shaping the beam and adding restraints at either end raised the stiffness by a factor of 69, as shown above. More generally, the transverse stiffness could always be written in the form k T = κk T0 , where κ was a coefficient that depended on the initial value of k. Figure 7 shows the numerically calculated variation of κ with kL. This fell monotonically from an initial value of 69, implying that the actuator became weaker and weaker as it deflected.

Conclusions
We presented analytic approximations to the results of the Euler theory that allowed the response of buckling-mode electrothermal actuators to be obtained in a closed-form, avoiding the need for an iterative solution of coupled analytic equations or a lengthy finite element analysis. The approximations were valid for actuators with very slender beams (such as NEMS designs) when the effects of bending dominated over those of axial compression. We compared the expressions with numerical solutions and showed them to be extremely accurate. We also presented higher-order corrections for less slender beams as found in MEMS actuators.
The analysis focused on two specific actuator layouts based on V-beam and raised cosine beam shapes. Each had the same set of material parameters ( , , and ) and geometric design parameters ( , , , and ). The careful choice of the latter allowed wellknown trade-offs between drive power, frequency response, and mechanical stiffness but the closed form Equations (11) and (26) and the results in Figure 3 showed that there was

Conclusions
We presented analytic approximations to the results of the Euler theory that allowed the response of buckling-mode electrothermal actuators to be obtained in a closed-form, avoiding the need for an iterative solution of coupled analytic equations or a lengthy finite element analysis. The approximations were valid for actuators with very slender beams (such as NEMS designs) when the effects of bending dominated over those of axial compression. We compared the expressions with numerical solutions and showed them to be extremely accurate. We also presented higher-order corrections for less slender beams as found in MEMS actuators.
The analysis focused on two specific actuator layouts based on V-beam and raised cosine beam shapes. Each had the same set of material parameters (E, α, and k th ) and geometric design parameters (W, L, D, and H). The careful choice of the latter allowed well-known trade-offs between drive power, frequency response, and mechanical stiffness but the closed form Equations (11) and (26) and the results in Figure 3 showed that there was little to choose between the two layouts for a fixed set of parameters. However, the approximation used (neglect of axial stress) may open the door to analytic design optimisation in more general cases; for example, when the beam width or other properties vary locally.