*3.1. Apparent Elastic Properties of Polysilicon Films*

Let us start by allowing for the uncertainties linked to the microstructure of the polysilicon film, namely the morphology of the polycrystal. As already stated, we focus here on the in-plane motion of the structure; therefore, a two-dimensional geometry was considered. To assess the impact of the microstructure on the dispersion of the results, an MC procedure was adopted; within such a stochastic procedure, each realization or sample of the polysilicon film had the shape of a square of side-length -. This length-scale was assumed to be linked to the characteristic size of the already mentioned critical details of the geometry of the movable structure. In the simple geometry of Figure 1, as the proof mass was assumed rigid, had to be driven by the springs' features: out of their length and in-plane width, the latter obviously drives scattering at the material level and is therefore considered in what follows. As far as the length was instead concerned, the way to handle it will be detailed in Section 3.2.

The microstructure inside each square sample was obtained via a (regularized) artificial Voronoi tessellation, by imposing an average grain size *g* = 1 μm [17] measured as the average distance between the centroids of the grains in the tessellation. Due to the columnar morphology of the polysilicon film [26], with a texture almost aligned with the out-of-plane axis for all the grains, the randomness of crystal lattice was assumed to hold in-plane only [19]. Because of the small -/*g* ratio, the length-scale separation assumption of homogenization did not hold true, and relevant asymptotic bounds on the elastic properties of the RVEs did not actually cover the full range of values of microstructure-affected apparent elastic properties: outliers therefore happened to fall out of the bilateral limits, as observed in [14]. The sampling of values out of the range provided by asymptotic, analytical bounds was assumed to undermine the assumption of representativeness, and the results thus had to be interpreted in a statistical context. To understand the above reasoning, in Figure 2 a few examples of the microstructural samples are shown, reporting also the orientation of the in-plane axes of elasticity for each grain.

In the numerical, Finite Element (FE) procedure feeding the MC analysis, the whole sample was meshed with quadratic triangular elements featuring a characteristic size of 125 nm, to assure accuracy in the results. This element size was already proven in [17,27] to lead to mesh-independent results, in terms of the homogenized elastic property of the polycrystalline samples. To decrease the computational burden of the simulations, especially in view of the number of samples to be handled in each MC analysis, some authors proposed the use of (grain) boundary element formulations; see e.g. [28–31], wherein only the grain boundary network has to be discretized to reduce the dimensionality of the problem.

**Figure 2.** The 2 <sup>×</sup> <sup>2</sup> <sup>μ</sup>m2 SVE examples of the in-plane polycrystalline morphology adopted in the MC analysis.

The reasons for relying on a numerical strategy with numerical simulations to feed an MC analysis were already discussed in [27]. If one is interested in the reference or mean values of the mechanical properties of the polysilicon film, asymptotic approaches will suffice; if the scattering induced by micromechanical features is instead the main focus of the analysis, a proper statistical description of all the stochastic variables and a proper sampling out of the relevant statistical distributions are indeed necessary. The aforementioned need is strongly enhanced in the considered case of -/*g* values on the order of 2–3. To this end an SVE, as handled in the analysis, was therefore defined as follows [32,33]: it is a domain whose dimension is smaller than a conventional RVE, but larger than the characteristic length scale (i.e. the grain size), so that the elastic properties obtained from the homogenization procedure are not to be considered as the effective ones to be adopted for a homogeneous deterministic device-level analysis, but as representative of the actual microstructure randomness. These elastic properties have to be considered as the "apparent" ones, according to the definition provided in [32]: i.e. the apparent Young's modulus obtained from each SVE is simply a realization of a random variable, which accounts for the relevant microstructure, and it should not be used in a deterministic analysis at the device-level scale.

For the MC analysis, two different sets of SVEs were handled with - = 2, 3 μm. These two values were selected in accordance with our previous studies [12–14]: to emphasize the microstructural effects, -/*g* has to be minimized in compliance with the microfabrication constraints. This rationale justified the value - = 2 μm, whereas - = 3 μm was selected on the basis of the dimension of beams in the structure described in [34].

The apparent elastic moduli of each sample were obtained according to the procedure discussed in [17]. A bilateral bounding scheme was adopted by imposing either uniform stress or uniform strain Boundary Conditions (BCs) all over the boundary of the SVE. Periodic BCs were not adopted, since a periodic microstructure can never be observed in real polysilicon films. For each type of BCs, three analyses were run, inducing in turn only one non-zero component of the macroscopic in-plane stress or strain; the results of these three analyses then allowed estimating the microstructure-affected values of *E*, *G*, and *ν*. The latter values turned out to depend on the film morphology due to the local intensification of the stress and strain fields caused by the misorientation of adjacent grains. Additional details, which are beyond the scope of this work, can be found in [17]; see also [27].

The outcomes of this statistical analysis are reported with continuous lines in Figures 3–5, in terms of the Cumulative Distribution Functions (CDFs) of *E*, *G*, and *ν*, as obtained with the two types of BCs and for the two SVE sizes. As is known, the uniform strain BCs lead to an upper bound on the actual value of *E* and *G* and a lower bound on the actual value of *ν*; vice versa, the uniform stress BCs provide a lower bound for *E* and *G* and an upper bound for the actual *ν*. The plots allow also to compare easily the 2 × <sup>2</sup> <sup>μ</sup>m<sup>2</sup> and the 3 × <sup>3</sup> <sup>μ</sup>m<sup>2</sup> cases and to quantify the effect of a smaller device-level size with respect to the average grain size: it marginally affects the mean values of the elastic moduli, while the scattering around them is strongly modified. The same results are also collectively reported in Figure 6 as box-whisker plots, to show in a compact way the mean values and the dispersion around them. No limits were assumed to exclude possible outliers in the data, so the full range of values of the results was included in the plots.

The mean value of the apparent Young's modulus was slightly dependent on the BCs (as expected, see [32] and see Table 1): for the 2 × <sup>2</sup> <sup>μ</sup>m2 case, it amounted to *<sup>E</sup>* = 149.98 GPa for uniform strain BCs, while it read *E* = 148.11 GPa for uniform stress BCs. It was almost negligibly affected (on the order of 1%) by the dimension of the SVE (- = 2 or 3 μm) for the assumed grain size *g* = 1 μm. The results showed instead a larger influence of the SVE size on the standard deviation, as qualitatively shown in Figure <sup>3</sup> by the less steep CDF curve for the 2 × <sup>2</sup> <sup>μ</sup>m<sup>2</sup> case with respect to the 3 × <sup>3</sup> <sup>μ</sup>m2 case. Quantitatively, the standard deviation for the uniform strain BCs increased from 3.34 GPa to 5.47 GPa when the SVE size decreased from - = 3 μm to - = 2 μm: the correspondent coefficients of variations (i.e. the ratio between the standard deviation and the mean value) were around 2.2% and 3.5%, respectively. We can interpret this larger scattering as the quantification of the effect of the

grain morphology on the statistical distribution of *E*, as the grain size becomes comparable with the SVE size.

Similar considerations can be made for the apparent shear modulus *G*; see Table 2: the scattering between the mean values was small (on the order of 2%) when the BCs and/or the SVE size were varied; the scattering became more significant if the standard deviation was considered, as evidenced by the values of the coefficients of variations, ranging from 3.9% (3 × <sup>3</sup> <sup>μ</sup>m2) to 6.2% (2 × <sup>2</sup> <sup>μ</sup>m2).

For the apparent Poisson's ratio (see Table 3), the mean values did not change significantly with the BCs, but they varied with the reduction of the SVE size of about 9%. The corresponding standard deviation instead showed a more meaningful change, varying from about 10% for the 3 × <sup>3</sup> <sup>μ</sup>m<sup>2</sup> case to about 17% for the 2 × <sup>2</sup> <sup>μ</sup>m<sup>2</sup> case.

As in any MC analysis, results relevant to a sufficiently high number of samples have to be collected to guarantee the convergence in the statistics of the apparent properties. With the focus on the homogenized Young's modulus of the polycrystalline film, some exemplary results are collected in Figure <sup>7</sup> for the 2 × <sup>2</sup> <sup>μ</sup>m2 SVE (top row), in terms of the convergence of the mean value and of the standard deviation with an increasing number *n* of samples. Besides some small fluctuations, visible only because of the reduced scale of the graphs, outcomes tended to converge and became stable even before *<sup>n</sup>* = 100. For the 3 × <sup>3</sup> <sup>μ</sup>m<sup>2</sup> SVE shown in Figure <sup>7</sup> (bottom row), the convergence was attained even faster, as expected for a larger SVE. Such an evolution of the graphs was not directly related to some specific SVE morphologies, since they remained unchanged (within the limits of the statistical representation) by shuffling the order of the SVEs handled in the MC analysis. Similar results were observed for any elastic property of the film, and for any SVE size; the only slight difference at varying SVE size was the different value at convergence for the standard deviation.

A further interesting feature of the mechanical response of the polysilicon film was its in-plane degree of anisotropy. In Figure 8, box-whisker plots are reported, again at varying BCs and SVE sizes, for the ratio *G*/*G*iso, where *G*iso = *E*/(2 (1 + *ν*)) would be the shear modulus of a perfectly in-plane isotropic material, whose Young's modulus *E* and Poisson's ratio *ν* were those obtained by means of the same homogenization procedure. This measure was already adopted in [17] to assess the material-dependent isotropy of microstructured materials and, to some extent, the representativeness of the SVE size. Since silicon, due to its FCC crystal structure, is a moderately-anisotropic material, the results related to aggregates of randomly-oriented grains are supposed to be characterized by a ratio *G*iso = *E*/(2 (1 + *ν*)) tightly bonded around isotropy, i.e. around *G*/*G*iso = 1. This ratio in the case of single crystal silicon, which has a low degree of anisotropy, turns out to be 1.57; it can be therefore appreciated that, in the case of a polycrystalline material, the ratio tends towards a unitary value, even for small SVEs as the ones considered in this work (see Figure 2).

To deal with the analysis of the imperfection sensitivity at the spring scale, the numerical statistical distributions were fitted with analytical lognormal CDFs (represented by the dashed lines in the graphs of Figures 3–5):

$$\mathcal{L}\left(\xi\right) = \frac{1}{2} + \frac{1}{2}\operatorname{erf}\left(\frac{\ln\xi - \mu}{\sqrt{2}\omega}\right) \tag{6}$$

where *μ* is the mean value and *ω* is the standard deviation of the natural logarithm of the random variable *ξ* (i.e. alternatively *E*, *G*, or *ν*); erf(·) is the normalized Gaussian function (also known as the "error function"). These distributions statistically satisfied the constraint for *E* and *G* to be always positive, with no artificial procedures necessary to avoid handling negative elastic moduli in the calculation of the overall spring stiffness; see also [35]. Moreover, an analytical, but SVE-informed CDF had the advantage of the ease of handling in further theoretical developments, as we will show for the spring stiffness in Section 3.2.

In Tables 1–3, the discussed mean and the standard deviation values obtained with the homogenization procedure are shown, together with the parameters *μ* and *ω* of Equation (6) as identified through a nonlinear fitting. These latter parameters made the analytical lognormal distribution ready to be used for the estimate of the stiffnesses of MEMS springs whose polysilicon morphology showed a grain size on the order of *g* = 1 μm. The quality of the adopted fitting via the lognormal distributions was quantified by the very small standard deviations shown in the tables for the parameters *μ* and for *ω* in all the considered cases: for the mean and standard deviation values corresponding to *E*, *G*, and *ν*, the scattering was even smaller and not affecting the least significant digit there reported.

The choice of the probability distribution for the fitting of the results of the MC analysis was a relevant decision. A lognormal probability distribution was assumed for each elastic modulus, but we assumed independence among *E*, *G*, and *ν*. This latter hypothesis does not prove always correct: in this way, in fact, one neglects the influence of the microstructure on the correlation among them. However, since in this work, we used the outcomes of the stochastic analysis to compute the stiffness of slender beams depending on the Young's modulus only, the hypothesis did not detrimentally affect the results. This would not be the case for more complex inferences, i.e. in the case of the stiffness evaluation for a shorter beam, where also *G* and shear deformations play a role. In these cases, a multivariate probability distribution should be considered, and the correlations between the variables should be evaluated, e.g. through marginal distributions. Assuming that the dependency structure of the multivariate probability distribution would be static, i.e. the stochastic process is stationary, copulae could be used to evaluate the said marginal probability distributions. For time-evolving stochastic processes, other approaches such as spectral analysis would be necessary. Both copulae and spectral analysis are out of the scope of this paper; interested readers can find information on them, e.g. in [36,37], respectively.

**Table 1.** Young's modulus *E*: apparent values of the mean and standard deviation obtained with the MC analysis and parameters for the lognormal CDFs in Equation (6).


**Table 2.** Shear modulus *G*: apparent values of the mean and standard deviation obtained with the MC analysis and parameters for the lognormal CDFs in Equation (6).


**Table 3.** Poisson's ratio *ν*: apparent values of the mean and standard deviation obtained with the MC analysis and parameters for the lognormal CDFs in Equation (6).


**Figure 3.** Effects of SVE size and BCs on the cumulative distribution function of the homogenized Young's modulus: (**a**) 2 <sup>×</sup> <sup>2</sup> <sup>μ</sup>m2 SVE and (**b**) 3 <sup>×</sup> <sup>3</sup> <sup>μ</sup>m<sup>2</sup> SVE. Continuous lines represent the results of MC analyses, whereas dashed lines are the relevant lognormal interpolants.

**Figure 4.** *Cont.*

**Figure 4.** Effects of SVE size and BCs on the cumulative distribution function of the homogenized shear modulus *<sup>G</sup>*: (**a**) 2 <sup>×</sup> <sup>2</sup> <sup>μ</sup>m<sup>2</sup> SVE and (**b**) 3 <sup>×</sup> <sup>3</sup> <sup>μ</sup>m<sup>2</sup> SVE. Continuous lines represent the results of MC analyses, whereas dashed lines are the relevant lognormal interpolants.

**Figure 5.** Effects of SVE size and BCs on the cumulative distribution function of the homogenized Poisson's ratio *<sup>ν</sup>*: (**a**) 2 <sup>×</sup> <sup>2</sup> <sup>μ</sup>m<sup>2</sup> SVE and (**b**) 3 <sup>×</sup> <sup>3</sup> <sup>μ</sup>m<sup>2</sup> SVE. Continuous lines represent the results of MC analyses, whereas dashed lines are the relevant lognormal interpolants.

**Figure 6.** Box-whisker plots of the homogenized elastic properties of the SVE: effects of the BCs (blue: uniform strain BCs; orange: uniform stress BCs) on the scattering of (**a**) Young's modulus, (**b**) shear modulus, and (**c**) Poisson's ratio.

**Figure 7.** *Cont.*

**Figure 7.** Effect of the number *n* of samples in the Monte Carlo analysis on the convergence of the (**a**,**c**) mean value and (**b**,**d**) standard deviation of in-plane Young's modulus *E*, for uniform strain BCs. Top: 2 <sup>×</sup> <sup>2</sup> <sup>μ</sup>m2 SVE. Bottom: 3 <sup>×</sup> <sup>3</sup> <sup>μ</sup>m<sup>2</sup> SVE.

**Figure 8.** Box-whisker plot for the ratio between *G*/*G*iso (blue: uniform strain BCs; orange: uniform stress BCs).

### *3.2. Overall Spring Stiffness Calculation*

For the suspension spring configuration shown in Figure 9, the overall stiffness entering into play in Equation (5) with the relevant scattering can be computed as follows. Results can be easily generalized to the case of folded beams, since the stiffness scales inversely proportional to the number of folds. Stiffness values reported for homogeneous beams (see [38]) cannot be adopted in the present statistical analysis due the local heterogeneity of the material that governs the problem.

We assumed the beam to be subdivided into *N* disjoint subdomains, like in domain decomposition procedures [39–41]: within each subdomain, Young's modulus *E* was assumed constant; as far as the moment of inertia was instead concerned, we assumed that, for the sake of simplicity, its value varied for each subdomain according to a Gaussian distribution accounting for the scattering of the over-etch. The mentioned value of *E* was sampled out of the lognormal fitting of the statistical distributions obtained from the analysis at the polycrystalline material level. In the analysis to follow, we will consider springs characterized by high values of the ratio between length and in-plane width, namely very slender or thin geometries: Young's modulus is then the only parameter affecting the solution, according to the mathematical procedure described in the Appendix. In the case of thick springs, shear stresses gain importance in the solution, and also the distribution of *G* should be accounted for; in such a case, the correlation between the distributions of *E* and *G* has to be also allowed for in the analysis.

**Figure 9.** Scheme adopted for the calculation of the scattered spring stiffness with the spring itself subdivided into *N* subsets, each one with a Young's modulus sampled from the relevant size-dependent CDF and a moment of inertia dependent on the over-etch.

The adopted procedure is schematically reported in Figure 9: a beam has a fixed constraint at one end, where displacements and rotations are prevented, and a slider at the opposite end, where the motion is in the direction perpendicular to the beam axis only, without any rotation. The latter end is supposed to be the link to the proof mass, which moves as a rigid body. By means of the principle of virtual work and allowing for beam slenderness to neglect shear strains, the force *F* and the corresponding displacement *u* can be related to provide the ratio *F*/*u* as the beam stiffness *k*.

The outcome of this reasoning, described in detail in Appendix A, depends on the number of subdivisions *N* along the beam length, and it turns out to provide the spring stiffness as:

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

where *L* is the beam length, *i* is the index running over the subsets handled, *Ii* is the moment of inertia of the *i*th subset, *Ei* is the apparent Young's modulus at the *i*th subset, and *ψ<sup>i</sup>* is a corrective factor dependent on the placement of the *i*th subset and defined in Equation (A10). Since the overall stiffness is a function of the random variables (*E*, *I*), then it is a random variable as well.

When the moment of inertia is assumed constant along the beam, the expression of the stiffness simplifies to:

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

### **4. Evaluation of the Offset at Rest**

### *4.1. Effect of Material Uncertainties Only*

In this section, the CDFs of the offset of the proof mass from its rest position are reported as obtained from Equation (5), in which *u* was taken as unknown, while the stiffnesses were calculated with two alternative choices: (i) from the lognormal CDFs as described in Section 3.2 or (ii) directly from FE simulations of the two-dimensional *<sup>L</sup>* × *<sup>w</sup>* (e.g. 200 × <sup>2</sup> <sup>μ</sup>m2) beam with the geometry shown in Figure 9. In these latter FE analyses, a random silicon grain morphology was considered for the whole domain: it is worthwhile to emphasize that the computational burden of these analyses was very high, since hundreds of grains were involved. To characterize the statistics of the stiffness via FE

analyses of the whole beam, 200 simulations were carried out by varying in each case the whole silicon grain morphology.

As for the CDF of *F* in Equation (5) resulting from the production process, we considered the experimental data shown in [42] measured in a 0.3 μm-thick polysilicon layer: a residual strain constant through the thickness *<sup>r</sup>* = <sup>360</sup> ± <sup>5</sup> · <sup>10</sup>−<sup>6</sup> was reported. In absence of further information, we assumed a Gaussian distribution for *<sup>r</sup>* with the mentioned mean and standard deviation, and we simply computed *F* = *E*˜ *<sup>r</sup> A*, *E*˜ = 150 GPa being a reference Young's modulus and *A* the beam cross-section area. Therefore, the distribution of *F* was Gaussian as well.

In Figure 10, the offset CDFs derived from the stiffnesses obtained through the sampling of the apparent Young's modulus from the (analytical) lognormal distributions are shown with orange lines, while the CDFs derived from the stiffnesses obtained from the FE simulations are instead reported with black lines. An almost zero offset mean value was common to the two cases, so the main difference between the two solutions was related to the variance of *E* and, therefore, of *k*. Looking at the results, we can conclude that the semi-analytical approach, which was far less expensive than the FE one, correctly represented the information obtained from the FE analyses. Around the zero mean value, the offset became positive or negative depending on the difference between the stiffnesses of the right and left springs (see Figure 1). It is noted that the offset was generated by the scattering of the values of the spring stiffnesses around the mean; therefore, it depended on the standard deviations of the CDFs, and not on the mean values. Any stochastic method providing an offset estimate should then address also the quantification of the variance of the random variables associated with the elastic properties, not only their mean values.

**Figure 10.** Offset from the rest position for length <sup>×</sup> SVE size cases: (**a**) 200 <sup>×</sup> <sup>2</sup> <sup>μ</sup>m<sup>2</sup> and (**b**) 200 <sup>×</sup> <sup>3</sup> <sup>μ</sup>m2.

An alternative representation is provided in Figure 11 for the 200 μm × 2 μm case. The histogram represents the percentage of a given realization with a certain value of offset from the rest position. Though not further investigated in this work, it can be seen that the distribution of the offset was not unimodal.

**Figure 11.** Histogram of the offset from the rest position for the 200 <sup>×</sup> <sup>2</sup> <sup>μ</sup>m2 case (uniform strain BCs).
