*4.2. Effect of Variable Over-Etch and Material Property Uncertainties*

In Figure 12, in addition to the effect of a variable Young's modulus according to the rationale described above, we took into account the scattering of the geometry due to the over-etch variation along the beam length. We assumed that the over-etch *o* along the beam was distributed as a Gaussian variable with a zero mean value and a standard deviation equal to 50 nm; for the *i*th domain along the beam, the corresponding moment of inertia was thus computed as *Ii* = *<sup>t</sup>*(*<sup>w</sup>* − <sup>2</sup> *oi*)3/12, where *oi* is the value of the over-etch sampled out of the aforementioned Gaussian distribution. *o* = 0 corresponds to the commonly-assumed over-etch reference value in the engineering practice: a constant value for a single MEMS structure or even for MEMS distributed at different positions on a wafer. The graphs in Figure 12 include, as a reference, the previous curves for FE and analytical results related to the effects due to the scattering of Young's modulus only: the scattering due to a variable moment of inertia was shown to be far larger, as evidenced by the lower slope of the dashed orange curve and by the corresponding standard deviation, which increased from 2.5–8.2 μm, while the mean value was about zero in all cases. Two additional curves, in blue lines, were also added to the graph: they represent the scattering of the offset, due to the uncertainty on the material only, when the two springs in Figure 1 differed in width exactly by a value of 150 nm, which is representative of the 95% confidence interval from the target during the manufacturing process [12,43]. It can be seen that most of the offsets fell within the range dictated by the two blue curves.

**Figure 12.** Effect on the offset of the variable (over-etch affected) moment of inertia *I* and Young's modulus along the beam length (200 <sup>×</sup> <sup>2</sup> <sup>μ</sup>m<sup>2</sup> case, uniform strain BCs). Blue lines refer to fixed values for the over-etch equal to ±150 nm.

### **5. Conclusions**

We proposed a stochastic method to evaluate the offset from the rest position for a statically-indeterminate MEMS typology, featuring a proof mass supported by two opposite springs. When the polysilicon production process was pushed to the current limit, the uncertainties about the stiffnesses of critical MEMS parts, such as nominally-identical supporting beams, generated an offset from the designed position when a residual stress distribution was also present inside the structure. Both material heterogeneity and fluctuations in critical geometry dimensions, i.e. the beam width, could be responsible for such an offset. For the first cause, we addressed the definition of the elastic properties accounting for the polycrystalline morphology by an artificially-reconstructed Voronoi tessellation in terms of statistically-representative elements, whose dimensions were equal to a characteristic dimension of the spring, namely their width, and were comparable with the silicon grain. For the second cause, we considered the effect of the over-etch on the moment of inertia of the MEMS supporting beams.

The classical homogenization procedure, under the hypothesis of statistical independence of the variables, provided a stochastic insight into the relevant mechanical properties, in particular the elastic moduli and Poisson's ratio, through numerical cumulated distribution functions that were subsequently approximated by an analytical, lognormal distribution. The stiffness of each beam supporting the moving mass was then estimated via the analytical lognormal distribution according to a simple reasoning based on the principle of virtual work applied to an exemplary beam.

By considering a statistically-indeterminate beam-moving mass-beam chain, we assumed as known the internal force due to residual stresses: it generated an offset from the rest position even in the presence of a stiffness varying only according to the polycrystalline material. A variable over-etch along the supporting beams, however, was responsible for a higher offset than the one induced by uncertainties on the material properties only. An important conclusion was that the evaluation of the offset depended on the correct characterization of the variance of Young's modulus and the moment of inertia, not on their mean value.

The approach shown here for a very simple class of structures could be extended to other statically-indeterminate MEMS configurations, provided that they are modeled as beam systems, in case the heterogeneity of the polysilicon is a matter of concern.

**Author Contributions:** The authors contributed equally to this work.

**Funding:** We gratefully acknowledge financial support through the project MaRe (Material Reliability) by STMicroelectronics.

**Acknowledgments:** The authors are also indebted to Marco Geninazzi for the results related to the stochastic homogenization procedure.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Abbreviations**

The following abbreviations are used in this manuscript:


### **Appendix A**

To quantify the stiffness of the supporting springs in the reference configuration of Figure 9, we relied upon the principle of virtual work for slender beams with a Bernoulli–Euler kinematics. When a force *F* acts perpendicularly to the longitudinal axis of the beam, it causes an internal bending moment M. At equilibrium, the external work *We* equals the internal work *Wi*, which read:

$$\mathcal{W}\_{\mathfrak{C}} = \mathcal{F}u \tag{A1}$$

$$\mathcal{W}\_{\mathbb{I}} = \int\_{0}^{L} \mathcal{M}(\mathbb{J}) \, \chi(\mathbb{J}) \, \mathrm{d}\mathbb{J} \tag{A2}$$

where *u* is the virtual displacement of force *F* in the same direction (hence, perpendicular to the beam axis), *χ* is the virtual curvature of the beam axis, and *ζ* is the axis aligned with the longitudinal direction of the spring.

For a statically-indeterminate structure, if both the elastic modulus *E* and the moment of inertia *I* are constant all along the beam axis, the solution to *We* = *Wi* can be found in any textbook of structural mechanics. In the present case, where both *E* and *I* can vary stochastically as functions of *ζ*, M(*ζ*) and *χ*(*ζ*) need to be determined. The variation of the bending moment along the spring axis reads M = *M* − *Fζ*, where *M* is the value of the top end of the beam, as induced by the local constraint that prevents any in-plane rotation. Hence, in compliance with the mentioned constraint, we have:

$$0 = \int\_0^L \frac{M - F\_\sharp^\gamma}{E(\zeta)|I(\zeta)|} \, d\zeta. \tag{A3}$$

The integral was solved numerically, by subdividing the beam into *N* disjoint domains of the same longitudinal length Δ*ζ* = *L*/*N*; within each of them, we assumed that the flexural rigidity was locally constant, i.e. stochastic fluctuations were locally neglected. Accordingly, we arrive at:

$$\sum\_{i=1}^{N} \frac{1}{E\_i \ I\_i} \int\_{\tilde{\zeta}\_i}^{\tilde{\zeta}\_i + \Delta \tilde{\zeta}} \left(M - F\zeta\right) \,\mathrm{d}\tilde{\zeta} = 0\tag{A4}$$

where *Ei* and *Ii* are the Young's modulus and the moment of inertia for the *i*th domain, respectively. By solving for *M*, we obtain:

$$\mathcal{M} = \left[ \frac{\sum\_{i=1}^{N} \frac{\Delta\_{2}^{\tau} + 2\bar{I}\_{i}}{2 \cdot E\_{i} \bar{I}\_{i}}}{\sum\_{i=1}^{N} \frac{1}{E\_{i} \bar{I}\_{i}}} \right] F. \tag{A5}$$

In compliance with the adopted subdivision of the beam, we also have *ζ<sup>i</sup>* = (*i* − 1)Δ*ζ*, and Equation (A5) becomes:

$$M = \frac{FL}{2}\gamma\tag{A6}$$

where:

$$\gamma = \frac{\sum\_{i=1}^{N} \frac{2i - 1}{N \left(E\_i \overline{I\_i}\right)}}{\sum\_{i=1}^{N} \frac{1}{E\_i \overline{I\_i}}} \tag{A7}$$

is a factor that allows for the non-uniform value of the flexural rigidity of the beam and modifies the value *M* = *F L*/2 relevant to the uniform beam case.

Now that the internal moments are fully given as a function of the external force *F*, we can solve *We* = *Wi* for the displacement *u* to obtain:

$$\mu = \int\_0^L \frac{F\left(\frac{1}{2}\gamma - \zeta\right)^2}{E(\zeta)I(\zeta)} \,\mathrm{d}\zeta\,\tag{A8}$$

being *χ* = *<sup>M</sup> <sup>E</sup>*(*ζ*) *<sup>I</sup>*(*ζ*). With the help of the same subdivision of the beam into *N* domains adopted to compute *M*, we finally obtain Equation (7), here rewritten for convenience:

$$k = \frac{F}{\mu} = \frac{1}{L^3} \left( \sum\_{i=1}^{N} \frac{\psi\_i}{E\_i \ I\_i} \right)^{-1} \tag{A9}$$

where:

$$
\psi\_i = \left(\frac{i-1}{N} - \frac{\gamma}{2}\right)^2 \frac{1}{N} + \left(\frac{i-1}{N} - \frac{\gamma}{2}\right) \frac{1}{N^2} + \frac{1}{3N^3}.\tag{A10}
$$

### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
