Stochastic E ﬀ ects on the Dynamics of the Resonant Structure of a Lorentz Force MEMS Magnetometer

: Resonance features of slender mechanical parts of Lorentz force MEMS magnetometers are a ﬀ ected by the (weakly) coupled thermo-electro-magneto-mechanical multi-physics governing their dynamics. We recently showed that reduced-order models for such parts can be written in the form of the Du ﬃ ng equation, whose nonlinear term stems from the mechanical constraint on the vibrations and is a ﬀ ected by the driving voltage. As some device performance indices vary proportionally to the amplitude of oscillations at resonance, an optimization of the operational conditions may lead to extremely slender, imperfection-sensitive movable structures. In this work, we investigate the e ﬀ ects of imperfections on the mechanical response of a single-axis magnetometer. At the microscopic length-scale, imperfections are given in terms of uncertainties in the values of the over-etch depth and of the Young’s modulus of the vibrating polycrystalline silicon ﬁlm. Their e ﬀ ects on the nonlinear structural dynamics are investigated through a Monte Carlo analysis, to show how the output of real devices can be scattered around the reference response trend.


Introduction
The rapid development of semiconductor technologies has enabled the mass production of plenty of micro electrical-mechanical system (MEMS) devices [1][2][3], among which the magnetometers considered in this study are relevant for compass applications. To design and predict the performance of these micro-devices, accurate models are required to account for the interaction among various physical fields, primarily the thermal, electric, magnetic, and mechanical ones. Numerical and analytical methods have been recently proposed to understand the underlying working principles, and increase the capability to match experimental evidences [4][5][6][7].
Standard approaches typically assume that the geometrical and physical properties of the device are known in a deterministic sense. In reality, especially with a view towards the goal of progressively reducing the device footprint, uncertainties look unavoidable and may become dominant due to a variety of factors linked to the micro-fabrication process [8]. For polysilicon MEMS, in [9][10][11][12][13][14] the effects of the crystalline morphology on the reliability of inertial devices subjected to impacts were studied. At the scale of the full device, it turned out that marginal effects can be actually reported regarding the failure mode, if any, as induced by micro-cracking of silicon. Looking instead more closely to the details of crack initiation and final percolation throughout the structural parts, small deviations may discriminate between failure and no-failure outcomes in different devices, which may differ only in terms of the film morphology in the failing region. Such result was validated against real experimental data [10].
Focusing on the operational conditions of micro-devices, an on-chip test procedure was proposed in [15,16]. To avoid accounting for the additional effects of residual stresses possibly induced by micro-fabrication, the movable structure was purposely designed as statically determinate. By slightly varying the geometry of the tested polysilicon samples and the actuation/sensing setup, the superposed effects of the polycrystalline elastic and geometric properties on the structural stiffness and, therefore, on the device response to actuation were assessed. Data reduction also provided a quantitative measure of the stochastic variability of Young's modulus of the polysilicon film, and of the over-etch depth [17,18].
Along the same line and with the aim to extend the preliminary results reported in [19], we deal here with the operational conditions of a single-axis (z-axis) Lorentz force MEMS magnetometer. As anticipated, this type of magnetometer has been developed for low-power electronic compass applications in handheld devices (see, e.g., [4,[20][21][22]). For this purposes, the MEMS must be able to sense the Earth's magnetic field, which amounts at most to about 100 µT. Indeed, sensitivity is not the only performance index to account for in the design, since power consumption is of paramount importance for mobile applications. In this regard, a so-called balanced optimization scheme was proposed in [6,7] to allow also for bandwidth and resolution of a capacitive, single-axis Lorentz force magnetometer similar to that here considered, see also [23][24][25].
To model the dynamics of the magnetometer resonant structure, which consists of a slender beam in a clamped-clamped configuration, the weak interaction among the thermal, electric, magnetic, and mechanical fields must be accounted for. The electric and magnetic fields interact to provide the Lorentz force inducing in-plane beam deformations, or vibrations in the case of a time-varying electric field flowing longitudinally in the beam. The Joule effect leads to a temperature rise in the beam that, being constrained at both ends, has a tendency to buckle if a critical threshold is attained. Though this instability has to be avoided, like pull-in under operational conditions, its proximity can be exploited to enhance the amplitude of oscillations, until a safe limit, and to simultaneously tune the working frequency of the device. A simplified model of this multi-physics frame has been developed by assuming a specific deformation mode for the beam, driven by its geometry and constraints. The differential equation governing the motion of the beam finally results to be of a Duffing type, where both the linear and nonlinear (cubic) stiffness terms embed the effects of the multi-physics, showing softening effects linked to the thermal load and to the electrostatic forces due to sensing. Uncertainties have been embedded by allowing for the scattered values of the Young's modulus and over-etch depth (see [15]) in a batch of statistically representative samples. To assess the effects of those uncertainty sources on the overall performance of the device, a series of Monte Carlo simulations has been run. Granted that the outcomes of this stochastic analysis cannot be automatically extended to any geometry, it is shown that over-etch has a greater impact on the solution, since it affects all the structural characteristics of the vibrating beam. In some cases, it may provide a geometry more prone to buckling, thus approaching the upper bounds on performance indicators whenever they are proportional to the compliance of the structure. Overall, if the compliance is maximized along with its effects on the (capacitive) sensing, the repeatability from part to part can be adversely affected and calibration becomes necessary to avoid drifts and biases [26,27]. Very recently, the proposal to embed machine learning capabilities into real-world devices has emerged as a disrupting trend (see, e.g., [28]); to train the neural network possibly driving a classifier, an approach can rely upon a simple stochastic model, like the one here proposed. Therefore, handling geometrical and mechanical parameters at the micro-scale in a statistical fashion looks necessary, as well as to develop a kind of embedded artificial intelligence for next generation devices.
The remainder of the paper is organized as follows: In Section 2, the multi-physics governing the problem is discussed in details and, moving from a weak form of energy conservation, the reduced-order model is obtained. To feed the stochastic analyses with data relevant to the scattering in the effective Young's modulus of the beam, a numerical homogenization procedure is discussed in Section 3. The results of Monte Carlo simulations, preceded by a sensitivity analysis, are collected in Section 4. Finally, Section 5 is aimed at highlighting the main achievements of the current work, and foreseeing future developments to avoid the approximations introduced here to simplify the analysis.

A Reduced-Order Model of the Resonating Structure
Different geometries were proposed for the resonating structures of Lorentz force MEMS magnetometers, each one specifically tailored for one-, two-, or three-axis compass applications (see, e.g., [4,29]). In this work, we focus on the configuration depicted in Figure 1: the movable parts of this MEMS are represented by a series of mechanically independent slender beams, connected to the substrate at their end cross-sections. Each beam carries a couple of plates for sensing purposes; to prevent as much as possible any global (rigid body-like) tilting mode due to asymmetric imperfections, such plates are connected to the mid-span cross-section of the beam. The considered device is able to sense a magnetic field aligned with the out-of-plane direction (z, according to the adopted reference frame), by exploiting the in-plane motion of the beam and, therefore, of the sensing plates. in Section 4. Finally, Section 5 is aimed at highlighting the main achievements of the current work, and foreseeing future developments to avoid the approximations introduced here to simplify the analysis.

A Reduced-Order Model of the Resonating Structure
Different geometries were proposed for the resonating structures of Lorentz force MEMS magnetometers, each one specifically tailored for one-, two-, or three-axis compass applications (see, e.g., [4,29]). In this work, we focus on the configuration depicted in Figure 1: the movable parts of this MEMS are represented by a series of mechanically independent slender beams, connected to the substrate at their end cross-sections. Each beam carries a couple of plates for sensing purposes; to prevent as much as possible any global (rigid body-like) tilting mode due to asymmetric imperfections, such plates are connected to the mid-span cross-section of the beam. The considered device is able to sense a magnetic field aligned with the out-of-plane direction ( , according to the adopted reference frame), by exploiting the in-plane motion of the beam and, therefore, of the sensing plates.  [19], and (b) scheme of the structure with the notation adopted.
The analysis to follow refers to a single beam, featuring length and a rectangular cross-section of area = ℎ, ℎ being the in-plane width and the out-of-plane thickness. Due to the end constraints, the beam vibrates in a clamped-clamped configuration.
In accordance with the micro-fabrication process, see e.g., [8], the beam is assumed to be made of a polycrystalline silicon film with a columnar structure. For the present analysis the length is assumed to be large enough, in comparison with a characteristic size of the polycrystalline morphology, to avoid considering local or higher-order effects on the dynamic properties of the beam. The beam is, therefore, assumed to be homogeneous, and the elastic properties governing its in-plane vibrations can be obtained through homogenization over a statistical volume element (SVE) of the polysilicon film; additional details are provided in Section 3.
As the beam is constrained at both ends, softening induced by thermal effects and driving voltage modifies the frequency at which the maximum amplitude of the vibrations is sensed. Out of the operational conditions, by increasing, e.g., the current density flowing in the flexure (see below), vibrations can be further amplified and buckling might get incepted; if the device is appropriately designed to enforce a small amplitude of the oscillations, pull-in instability [30] could be instead easily avoided. The elastic response of the beam is thus modeled accounting for second-order effects, see [31][32][33], so that the lateral deflection does affect the equilibrium state. Due to the beam slenderness, as measured by the ratio ℎ ⁄ , energy conservation is enforced in weak form through:  [19], and (b) scheme of the structure with the notation adopted.
The analysis to follow refers to a single beam, featuring length L and a rectangular cross-section of area A = bh, h being the in-plane width and b the out-of-plane thickness. Due to the end constraints, the beam vibrates in a clamped-clamped configuration.
In accordance with the micro-fabrication process, see e.g., [8], the beam is assumed to be made of a polycrystalline silicon film with a columnar structure. For the present analysis the length L is assumed to be large enough, in comparison with a characteristic size of the polycrystalline morphology, to avoid considering local or higher-order effects on the dynamic properties of the beam. The beam is, therefore, assumed to be homogeneous, and the elastic properties governing its in-plane vibrations can be obtained through homogenization over a statistical volume element (SVE) of the polysilicon film; additional details are provided in Section 3.
As the beam is constrained at both ends, softening induced by thermal effects and driving voltage modifies the frequency at which the maximum amplitude of the vibrations is sensed. Out of the operational conditions, by increasing, e.g., the current density flowing in the flexure (see below), vibrations can be further amplified and buckling might get incepted; if the device is appropriately designed to enforce a small amplitude of the oscillations, pull-in instability [30] could be instead easily avoided. The elastic response of the beam is thus modeled accounting for second-order effects, see [31][32][33], so that the lateral deflection does affect the equilibrium state. Due to the beam slenderness, as measured by the ratio L/h, energy conservation is enforced in weak form through: where the integrals respectively denote the first-order and second-order terms of the internal energy, the kinetic energy and the work of external forces. In Equation (1), that must hold for any variation δv of the beam deflection and of the relevant space and time derivatives, an Euler-Bernoulli kinematics has been enforced, so that beam cross-sections are always assumed to maintain their planar geometry and their orthogonality with the deflected longitudinal axis. Terms appearing in the conservation law are: v(x, t), t being time, is the beam deflection; space derivatives v = ∂v/∂x and v = ∂ 2 v/∂x 2 represent the rotation and the curvature of the beam axis; time derivative .. v = ∂ 2 v/∂t 2 represents the lateral acceleration; EI is the flexural rigidity of the beam, E being the effective in-plane polysilicon Young's modulus and I = 1 12 bh 3 the moment of inertia of the beam cross-section; η = A is the mass per unit length, being the mass density; P is the axial load, positive if tensile, caused by thermal effects; and { is the magnitude of the lateral load per unit length, due to the external actions.
The energy conservation equation is solved by imposing that the beam deflection can be written as: thereby assuming in the solution a multiplicative decomposition into space and time functions. The resulting single degree of freedom V(t) represents the time evolution of the amplitude of oscillations at the mid-span cross-section. The differential equation governing the time evolution of V can then be recast in the form of Duffing equation, as follows: where: Terms in Equation (4) respectively represent the effective mass, damping, linear, and cubic stiffness, and external load (Lorentz force) ones, see also [34]. In these equations: η * is the mass per unit length of the attached sensing plates, which might be different from η if beam and plates are designed to have different in-plane widths; µ is the effective viscosity coefficient of the surrounding air; g is the gap between the two surfaces of each parallel-plate capacitor; α is the coefficient of (longitudinal) thermal expansion; i is the density of the current, flowing longitudinally in the beam; ∆T is the mid-span temperature raise due to Joule effect; V is the driving voltage; ε 0 is the permittivity of vacuum; and B is the magnitude of the out-of-plane magnetic field to be sensed.
In this formulation, as proposed in [5], the mass term accounts for the contribution of the two sensing plates too, as they move together with the beam and behave like rigid bodies to interact with the massive stators at top and bottom sides, see Figure 1. The damping term is given by the squeeze of the fluid at the sensing electrodes, see [3,35]. As for the stiffness terms K 1 and K 3 , the aforementioned softening effects are measured by the temperature raise caused by the Joule effect, and by the electrostatic forces at the parallel plate capacitors. An increase of the current density i or of the driving voltage V may, thus, increase the amplitude of oscillations, and also the sensitivity to imperfections in the response of the device. Instability phenomena, like mechanical buckling or pull-in, must be avoided under operational conditions: in [5], a trade-off was automatically set in between increasing the magnetometer performances and decreasing the risk of instabilities, by constraining the resonance frequency of the movable structure within a standard range for this type of applications.
In the Duffing equation (Equation (3)), the term K 3 quantifies the deviation from linearity in the system dynamics. As it is affected by V 2 (see Equation (4)), the relevant dependency on the driving mechanism needs to be accurately assessed in the solution. We then assume that the forcing term varies according to F = F cos ωt, where ω is the corresponding frequency; the maximum amplitude V max of the oscillations is then provided as the solution of the following equation [32]: where ω 1 = √ K 1 /m is the natural frequency relevant to the linearized beam response. The analytical solution of Equation (5) for ω = ω 1 was already discussed in [5]. The main point to stress here is that such solution obviously depends on K 1 and K 3 which, in turn, depend on the flexural and axial rigidities of the beam. Hence, uncertainties in the values of the polysilicon Young's modulus E and of the beam width h induce a scattering in the value of V max , that has to be properly quantified at the design stage. For this purpose, in what follows, we provide a model to account for the aforementioned uncertainties at the polycrystalline level. For what concerns the Young's modulus E of the micro-structured silicon film, results of stochastic homogenization are discussed in Section 3. As the width h of the film depends on the over-etch depth induced by the micro-fabrication process (see [15][16][17]) a discussion pertinent to this effect is provided in Section 4.

Effective Elastic Properties of the Polysilicon Film
In the formulation provided in Section 2, the elastic properties of the polysilicon film come explicitly into play through the Young's modulus E in the longitudinal direction x. Though not accounted for here, the out-of-plane beam thickness b or, effectively, the ratio between this thickness and the in-plane beam width h, may lead to plain strain conditions to hold in the micro-structured film; the bending stiffness should be therefore appropriately modified, to account for the additional effect of the Poisson's ratio on the deformation field.
As the stochastic effects of the polycrystalline morphology on the elastic moduli of the film have to be evaluated, we propose two different strategies to obtain quantitative estimations. A former one is aimed at bilaterally bounding the asymptotic solution for samples of the film featuring a size growing to infinity; this is obtained through an analytical homogenization procedure. A latter one is instead aimed at assessing the morphology-induced scattering in the solution for the finite size beam geometry under consideration, and is obtained via a Monte Carlo-driven, numerical homogenization procedure. The two approaches were already adopted in [36] (see also [14]), and they are briefly discussed in what follows to assure a self-consistent description of the whole method.
In both procedures, we account for the columnar micro-structure of the epitaxially grown polycrystalline silicon film. To simplify the analysis, we assume that a texture is aligned with the out-of-plane direction, namely is perpendicular to the substrate: therefore, homogenization is carried out in the x − y plane only (see Figure 1) and two-dimensional digital representations of the film morphology are adopted in the analysis. As detailed in several studies (see e.g., [8,9]) this assumption looks appropriate for the considered micro-fabrication process, leading to a sub-vertical growth of all the silicon grains.

Analytical, Deterministic Homogenization
As anticipated, bounds obtained with this approach are valid asymptotically, for a ratio between the in-plane dimensions of the polysilicon sample and the characteristic grain size (or radius) s g tending to infinity. The procedure is therefore intended to provide bilateral bounds on the mean values of the Young's modulus of polysilicon, since bounds result to be too tight if the ratio between a characteristic size of the beam, like its width h, and s g does not take large values.
Bounds on the elastic moduli are based on the Voigt and Reuss approaches, which assume that either the strain or stress state is uniform throughout the polycrystalline sample [37]. Moving from the Hill-Mandel macro-homogeneity condition, under the Voigt assumption of a uniform strain field, the effective stiffness matrix is bounded by: Under the Reuss assumption of a uniform stress field, the following bound on the effective compliance matrix is instead obtained: In these equations: Ω is (the measure of) the polysilicon sample volume that is, in the present two-dimensional case, the in-plane area; c l is the in-plane single-crystalline silicon stiffness matrix in a local reference frame aligned with the axes of elastic symmetry, which reads: where c 1111 l = 165.7 GPa, c 1122 l = 63.9 GPa and c 1212 l = 79.6 GPa (see, e.g., [9]); t σ and t ξ are the orthogonal transformation matrices, that allow to respectively transform the stress and strain vectors from the global reference frame of Figure 1 to the local one aligned with the aforementioned axes of elastic symmetry. According to a standard notation, superscripts T and −1 respectively refer to the transpose and to the inverse of the relevant matrices.
To compute the bounds, we assume that the in-plane lattice orientation of each grain is a fluctuating field with a uniform statistical distribution. The integral over Ω in Equations (6) and (7) can, therefore, be transformed into an integral over such random orientation. The resulting overall behavior of the polysilicon film turns out to be in-plane isotropic, with bounds on the Young's modulus E given by: E V = 151.4 GPa, as furnished by the Voigt method; E R = 147.1 GPa, as furnished by the Reuss method.
The assumption of a uniform distribution for the lattice orientation was already proved to fail in case of real film geometries, like those here considered, since the number of grains in the film sample is not large enough [38]. Therefore, the study to follow is aimed at providing a stochastic description of the scattering in the beam stiffness.

Numerical, Stochastic Homogenization
In the current work, the length-scale separation principle [39] adopted to define the properties of a representative volume element of the polycrystalline film does not hold true, due to the small values of the h/s g ratio. Because of the expected growth of the effects on the results of micro-structural features of the polycrystalline morphology, an SVE is adopted to feed a Monte Carlo procedure aiming to estimate the scattering in the Young's modulus E of the film, see Figure 2. Perturbations on the SVE geometry are provided in terms of topology of the grain boundary network, and in terms of in-plane orientation of the symmetry axes of each silicon grain.
Each micro-structured sample is numerically generated using a (regularized) Voronoi tessellation, see [9,14,36], wherein a target grain size s g is assumed. To bound the elastic response of each SVE, either uniform strain of uniform stress boundary conditions (BCs) are adopted. The mentioned BCs are aimed to define the displacement or the traction field along the lateral boundary of the SVE, and do not provide any additional constraint on the local stress and strain fields inside. The overall SVE response is computed by defining the macroscopic stress Σ and strain Ξ vectors as the volume integrals of the relevant microscopic, fluctuating fields, and by linking them linearly through Σ = CΞ, where C is the same stiffness matrix considered before and in need of an estimation. To compute the entries of C, a procedure was proposed in [36] on the basis of a set of finite element simulations, each one featuring only one component of either Σ or Ξ different from zero, depending on the type of BCs handled. Procedural details are omitted here, and interested readers can find them in [36]. The cumulative distribution functions bounding , as obtained with the two types of BCs and for ℎ = 2 − 3 μm and = 0.5 μm, are shown in Figure 3 together with the relevant analytical lognormal interpolations. The lognormal distribution, whose details in the cases under consideration are reported in Table 1, has been selected in place of, e.g., a Gaussian or normal one, since it automatically satisfies the thermodynamic constraint for which must be positive. Accordingly, it is deterministically set that unphysical, negative values of the Young's modulus of the film cannot be sampled from the distributions obtained. To provide a comparative assessment between the current bounds and the former deterministic ones given in Section 3.1, the vertical lines in Figure 3 provide such latter tight ones.
(a) The cumulative distribution functions bounding E, as obtained with the two types of BCs and for h = 2 − 3 µm and s g = 0.5 µm, are shown in Figure 3 together with the relevant analytical lognormal interpolations. The lognormal distribution, whose details in the cases under consideration are reported in Table 1, has been selected in place of, e.g., a Gaussian or normal one, since it automatically satisfies the thermodynamic constraint for which E must be positive. Accordingly, it is deterministically set that unphysical, negative values of the Young's modulus of the film cannot be sampled from the distributions obtained. To provide a comparative assessment between the current bounds and the former deterministic ones given in Section 3.1, the vertical lines in Figure 3 provide such latter tight ones. The cumulative distribution functions bounding , as obtained with the two types of BCs and for ℎ = 2 − 3 μm and = 0.5 μm, are shown in Figure 3 together with the relevant analytical lognormal interpolations. The lognormal distribution, whose details in the cases under consideration are reported in Table 1, has been selected in place of, e.g., a Gaussian or normal one, since it automatically satisfies the thermodynamic constraint for which must be positive. Accordingly, it is deterministically set that unphysical, negative values of the Young's modulus of the film cannot be sampled from the distributions obtained. To provide a comparative assessment between the current bounds and the former deterministic ones given in Section 3.1, the vertical lines in Figure 3 provide such latter tight ones.
(a)  An interesting feature of the results is that the mean values of the computed distributions are very close to the average between the Reuss and Voigt bounds. What is instead majorly affected in the solution is the scattering of the values around the mean, and the dependence of the standard deviation on ℎ ⁄ (see also [36]).

Monte Carlo Analysis: Sensitivity to Imperfections
Stochastic effects on the device sensitivity, which results to be proportional to the amplitude of the oscillations of the resonant structure, are now assessed. These effects play a role in defining confidence intervals for additional performance indices like, e.g., power consumption, bandwidth, and resolution, see [5][6][7]. The focus of this study is anyway on micro-scale imperfections, as measured by the scattering in the beam width ℎ and in the effective Young's modulus of the polysilicon film around the target values.
Micro-structural effects on the value of have been already investigated in Section 3; microfabrication effects on ℎ are instead handled in terms of an over-etch depth . By means of purposely designed on-chip testing devices featuring very slender movable structures, we already quantified the impact of on the aforementioned scattering. In fact, changes the in-plane film width, whose value turns out to be ℎ − 2 while ℎ represents only a target size; accordingly, the cross-sectional area of the beam is affected linearly by , whereas the moment of inertia is affected cubically. The over-etch depth also has primary indirect effects on : the gap at the capacitors, if the depth is assumed to be isotropically distributed within the wafer plane, is given by + 2 , where again represents only a target value. The temperature raise Δ at mid-span, given by: where, besides the variables already introduced in Section 2, and respectively represent the polysilicon resistivity and thermal conductivity, is affected as well due to its dependence on the crosssectional area . An interesting feature of the results is that the mean values of the computed distributions are very close to the average between the Reuss and Voigt bounds. What is instead majorly affected in the solution is the scattering of the values around the mean, and the dependence of the standard deviation E s on h/s g (see also [36]).

Monte Carlo Analysis: Sensitivity to Imperfections
Stochastic effects on the device sensitivity, which results to be proportional to the amplitude V max of the oscillations of the resonant structure, are now assessed. These effects play a role in defining confidence intervals for additional performance indices like, e.g., power consumption, bandwidth, and resolution, see [5][6][7]. The focus of this study is anyway on micro-scale imperfections, as measured by the scattering in the beam width h and in the effective Young's modulus E of the polysilicon film around the target values.
Micro-structural effects on the value of E have been already investigated in Section 3; micro-fabrication effects on h are instead handled in terms of an over-etch depth o. By means of purposely designed on-chip testing devices featuring very slender movable structures, we already quantified the impact of o on the aforementioned scattering. In fact, o changes the in-plane film width, whose value turns out to be h − 2o while h represents only a target size; accordingly, the cross-sectional area A of the beam is affected linearly by o, whereas the moment of inertia I is affected cubically. The over-etch depth also has primary indirect effects on V max : the gap at the capacitors, if the depth o is assumed to be isotropically distributed within the wafer plane, is given by g + 2o, where g again represents only a target value. The temperature raise ∆T at mid-span, given by: where, besides the variables already introduced in Section 2, ψ and K H respectively represent the polysilicon resistivity and thermal conductivity, is affected as well due to its dependence on the cross-sectional area A. Before moving to the results of the Monte Carlo simulation, a sensitivity analysis is provided to assess the dependence of the dynamic response of the beam on E and o. In this analysis, the values of all the parameters have been set as gathered in Table 2, see also [5]. The value of the magnetic field intensity, with the aim discussed in [29] of developing devices for compass applications, has been assumed to be B = 50 µT. Table 2. Values of the geometrical, physical, actuation, and sensing parameters adopted in the analysis.

Property
Value 2 mass density (kg/m 3 ), [40] 2330 thermal expansion coefficient α (K −1 ), [40] 2.5 × 10 −6 thermal conductivity K H (W/mK), [41] 25 resistivity ψ (Ωm), [41] 10 −5 permittivity ε 0 (F/m) 8.85 × 10 −12 air viscosity µ (Ns/m 2 ), [35] 5 × 10 −8 actuation current i (mA) 1 bias voltage V (V) 2 magnetic field intensity B (µT), [4] 50 The impact of E and o on the response of the movable structure is shown in Figure 4, in terms of V max for both the target values of the beam thickness considered in the stochastic homogenization procedure of Section 3. In the plots, results are reported in a dimensionless form, to focus more on the trends at varying parameter values, rather than on the amplitude itself; for both the considered geometries, the amplitude of oscillations is therefore divided by a reference value V * max linked to E * = (E R + E V )/2 =149.25 GPa and o = 0. In the considered cases, this reference amplitude amounts to: V * max (h = 2 µm) =0.028 µm; V * max (h = 3 µm) =0.016 µm. The effect of the increased beam stiffness for h = 3 µm is clear in these values; keeping aside such difference between the reference values, the dimensionless plots of Figure 4 allow to catch the trends in the system response. E affects the linear stiffness K 1 and, therefore, the natural frequency f 1 = ω 1 /2π at which V max is computed; the solution in terms of the oscillations amplitude accordingly turns out to be a weakly nonlinear function of the beam elastic properties, as shown in Figure 4a. The impact of E on the solution is rather small; in order to report a variation still below 10% around the reference one, the range of values of the Young's modulus has been extended to include all those possible for an ideal single-crystalline silicon sample, that is 130 GPa ≤ E ≤ 169 GPa [40]. To better assess the smaller variability for a polycrystalline aggregate, independently of the target value h of the film width, the analytical bounds E R and E V are also displayed with the vertical blue dotted lines. A much larger effect is due to the depth o of the over-etch: Figure 4b shows that a variation of up to about 50% around V * max can be attained in the range -0.15 µm ≤ o ≤ 0.15 µm, as selected for the considered micro-fabrication process. Terms m, d, K 1 and K 3 are all affected by o, and the solution thus shows the nonlinear dependence depicted in the graph, where the vertical blue dotted line is reported just to identify the condition of null o leading to the reference value V * max . In the adopted nondimensional format, the solution relevant to h = 2 µm, though characterized by a smaller amplification for positive o values, shows a remarkable nonlinearity induced by the Joule effect. This is reported also in Figure 5, where the mid-span temperature raise ∆T is plotted against the over-etch depth; once again, results are provided in a dimensionless form to better assess how close the actuation parameters i and V in Table 2 are to those leading to buckling. For this configuration, the critical value ∆T cr of the temperature raise is given by: ∆T cr = 6π 2 I αAL 2 (10) and, due to its dependence on I and A, it is affected by o as well. A reduction of the film width, given by positive values of o, leads to an increase of ∆T and a simultaneous decrease of ∆T cr , overall causing the nonlinear variation of ∆T/∆T cr shown in Figure 5 for both the target beam widths. The dependence is obviously more pronounced for h = 2 µm, so that an increase of the excitation current i might induce buckling by exceeding the threshold ∆T/∆T cr = 1, especially in the stochastic environment to be considered next.
Additional results of the parametric analysis are reported in Figure 6, in terms of the variation of the fundamental resonance frequency f 1 of the beam due to E and o. Results are obtained by dividing the frequency f 1 by the reference one f * 1 linked to E * and o = 0, which amounts to: f * 1 (h = 2 µm) = 15.5 kHz; f * 1 (h = 3 µm) =28.9 kHz. Alike the results of Figure 4, the solution is slightly affected by E and varies in a highly nonlinear manner due to o for h = 2 µm. In this case, the reduction of f 1 for values approaching o = 0.15 µm is another indicator of the frequency tuning provided by system nonlinearities, with a larger reduction being an indicator of the proximity to buckling instability. This condition has been purposely selected, to show next how all the nonlinear effects have a role in setting the scattering of V max in a stochastic environment. and, due to its dependence on and , it is affected by as well. A reduction of the film width, given by positive values of , leads to an increase of Δ and a simultaneous decrease of Δ , overall causing the nonlinear variation of Δ /Δ shown in Figure 5 for both the target beam widths. The dependence is obviously more pronounced for ℎ = 2 μm, so that an increase of the excitation current might induce buckling by exceeding the threshold Δ /Δ = 1 , especially in the stochastic environment to be considered next.
Additional results of the parametric analysis are reported in Figure 6, in terms of the variation of the fundamental resonance frequency of the beam due to and . Results are obtained by dividing the frequency by the reference one * linked to * and = 0, which amounts to: * (ℎ = 2 μm) =15.5 kHz; * (ℎ = 3 μm) =28.9 kHz. Alike the results of Figure 4, the solution is slightly affected by and varies in a highly nonlinear manner due to for ℎ = 2 μm. In this case, the reduction of for values approaching = 0.15 μm is another indicator of the frequency tuning provided by system nonlinearities, with a larger reduction being an indicator of the proximity to buckling instability. This condition has been purposely selected, to show next how all the nonlinear effects have a role in setting the scattering of in a stochastic environment.     Results of the Monte Carlo analysis are now discussed. As far as is concerned, for each value of ℎ the lognormal distribution best fitting the results of the numerical homogenization of Section 3 has been adopted; we recall that, depending on the type of BCs, relevant values of mean and standard deviation are reported in Table 1. As for the over-etch, in accordance with former studies [15,16] its distribution has been assumed to be Gaussian, centered around the reference condition = 0; the corresponding standard deviation has been set to = 0.05 μm, so that the range −0.15 μm ≤ ≤ 0.15 μm inspected in the previous parametric analysis represents the 99.7% confidence interval. The output of the stochastic analysis is shown in Figures 7 and 8 for ℎ = 2 μm and ℎ = 3 μm, respectively, in terms of the nondimensional values of . In these graphs, each dot represents a  Results of the Monte Carlo analysis are now discussed. As far as is concerned, for each value of ℎ the lognormal distribution best fitting the results of the numerical homogenization of Section 3 has been adopted; we recall that, depending on the type of BCs, relevant values of mean and standard deviation are reported in Table 1. As for the over-etch, in accordance with former studies [15,16] its distribution has been assumed to be Gaussian, centered around the reference condition = 0; the corresponding standard deviation has been set to = 0.05 μm, so that the range −0.15 μm ≤ ≤ 0.15 μm inspected in the previous parametric analysis represents the 99.7% confidence interval. The output of the stochastic analysis is shown in Figures 7 and 8 for ℎ = 2 μm and ℎ = 3 μm, respectively, in terms of the nondimensional values of . In these graphs, each dot represents a Results of the Monte Carlo analysis are now discussed. As far as E is concerned, for each value of h the lognormal distribution best fitting the results of the numerical homogenization of Section 3 has been adopted; we recall that, depending on the type of BCs, relevant values of mean and standard deviation are reported in Table 1. As for the over-etch, in accordance with former studies [15,16] its distribution has been assumed to be Gaussian, centered around the reference condition o = 0; the corresponding standard deviation has been set to o s = 0.05 µm, so that the range −0.15 µm ≤ o ≤ 0.15 µm inspected in the previous parametric analysis represents the 99.7% confidence interval.
The output of the stochastic analysis is shown in Figures 7 and 8 for h = 2 µm and h = 3 µm, respectively, in terms of the nondimensional values of V max . In these graphs, each dot represents a solution of the Monte Carlo simulation handling 10,000 samples. In the charts, the horizontal black dashed lines stand for the limits of the aforementioned confidence interval; the outliers falling outside such interval are clearly reported in the graphs versus E, while they are partially dropped in the graphs versus o to avoid widening too much the displayed interval beyond the previously discussed limits. The vertical blue dotted lines still represent the size-independent bounds E R and E V on the polysilicon Young's modulus, and the reference no-over-etch case. The orange lines are, as in Figure 4 To focus more on the quantitative analysis of the results, Figure 9 provides the cumulative distribution functions of V max for the two film widths considered. The relevant values of mean and standard deviation of V max , as affected by film width and type of BCs, are all collected in Table 3. As shown also by the graphs, the effect of the type of BCs is marginal: the curves relevant to uniform stress and uniform strain BCs can be hardly distinguished at the reported global scale. Once more, this is a proof of the much smaller impact of E on the solution and, therefore, on the device performance, in comparison to that of o. solution of the Monte Carlo simulation handling 10,000 samples. In the charts, the horizontal black dashed lines stand for the limits of the aforementioned confidence interval; the outliers falling outside such interval are clearly reported in the graphs versus , while they are partially dropped in the graphs versus to avoid widening too much the displayed interval beyond the previously discussed limits. The vertical blue dotted lines still represent the size-independent bounds and on the polysilicon Young's modulus, and the reference no-over-etch case. The orange lines are, as in Figure  4, the trends provided by a varying at = 0, and by a varying at = * . They are depicted to catch the relative importance of and on the scattering in the solution: while all the sampled outcome values lie close to the trend given by alone, they are much spread around the trend induced by . This is a clear qualitative measure of the higher importance of the role played by , as compared to that of .
To focus more on the quantitative analysis of the results, Figure 9 provides the cumulative distribution functions of for the two film widths considered. The relevant values of mean and standard deviation of , as affected by film width and type of BCs, are all collected in Table 3. As shown also by the graphs, the effect of the type of BCs is marginal: the curves relevant to uniform stress and uniform strain BCs can be hardly distinguished at the reported global scale. Once more, this is a proof of the much smaller impact of on the solution and, therefore, on the device performance, in comparison to that of .   To quantify next the sensitivity of the magnetometer, we refer to the differential sensing scheme provided by the two parallel plate capacitors depicted in Figure 1, see [4,7]. The magnitude ∆C of the capacitance change due to the deformation of the beam, or flexure is given by: where ∆g represents variation of the gap at the capacitors, or the amplitude of the lateral oscillations of the beam. By driving the device at frequency f 1 , at which V max represents the said variation of the gap, and by normalizing the results with respect to the intensity of the external magnetic field to sense, the following measure of the mechanical sensitivity is arrived at: In accordance with the model discussed in Section 2, in Equation (12) it has been assumed that the gap at capacitors is homogeneously varied all over their surfaces. Hence, the plates attached to the beam get displaced as rigid bodies and, due to the geometrical symmetry, tilting of the plates is disregarded. Moreover, even though Figure 1 shows that the in-plane length of the stators is smaller than the length L of the beam, the two measures have been assumed the same; results can be anyway generalized since, according to Equation (12), φ s is proportional to L.
For the two target values of the beam width h, reference values of φ s turn out to be: φ * s (h = 2 µm) =31.35 pF/T; φ * s (h = 3 µm) =18.47 pF/T. The scattering in the sensitivity provided by the Monte Carlo analysis is shown in Figure 10, in terms of the cumulative distribution functions of φ s . In these graphs, the vertical blue dotted lines represent the above mentioned h-dependent reference values of the sensitivity. Statistics are also collected in Table 4, in terms of the values of mean and standard deviation of φ s , depending on film width and type of BCs. The results testify that, due to the nonlinearities in the definition of the sensitivity and in the system response to the magnetic field, and to the scattering in the values of V max and g, the cumulative distribution functions of the output result to be distorted in comparison with those relevant to the input parameters. Due to the compliance of the beam in the h = 2 µm case, the curves relevant to the two types of BCs are not superposed, and can be distinguished at the scale of graph. At any rate, the plot clearly shows that, also for this performance index, the variation in the statistics linked to the BCs plays a marginal role, being largely superseded by the scattering induced by the two main parameters E and o handled in the Monte Carlo analysis.

Conclusions
In this work, we have discussed a single degree of freedom reduced-order model for the dynamics of the resonant structure of a single-axis Lorentz force MEMS magnetometer. By accounting for the (weakly) coupled thermo-electro-magneto-mechanical multi-physics governing the vibrations of the slender clamped-clamped resonant beam, its motion turned out to be described by the Duffing equation, where nonlinearities are enhanced by the end constraints on the motion itself.
Since the Joule effect, caused by the current flowing along the longitudinal axis of the beam, induces a temperature rise that may get close to values leading to buckling instability (while pull-in is avoided due to the small amplitude of the vibrations of the flexure), the characteristics of the structural vibrations turn out to be highly affected by micromechanical imperfection-like fluctuations of stiffness and geometry of the moving polysilicon film.
A statistical analysis has been developed, in order to assess the sensitivity of the sensor response to uncertainties in the values of the over-etch depth and of the morphology-driven overall Young's modulus of the polycrystalline silicon. A Monte Carlo investigation has been carried out, wherein the over-etch depth has been sampled out of a micro-fabrication tailored normal distribution, and the Young's modulus has been sampled out of a size-dependent lognormal distribution obtained by best fitting the outcomes of a stochastic homogenization on samples (or statistical volume elements) of the polycrystalline film.
A sensitivity analysis has shown that the amplitude of oscillations of the resonant beam of this Lorentz force MEMS magnetometer is marginally affected by the elasticity of polysilicon, while it strongly depends on the over-etch depth. Similar results have been arrived at with the Monte Carlo simulations, which have highlighted that the scattering in the beam response and in the mechanical

Conclusions
In this work, we have discussed a single degree of freedom reduced-order model for the dynamics of the resonant structure of a single-axis Lorentz force MEMS magnetometer. By accounting for the (weakly) coupled thermo-electro-magneto-mechanical multi-physics governing the vibrations of the slender clamped-clamped resonant beam, its motion turned out to be described by the Duffing equation, where nonlinearities are enhanced by the end constraints on the motion itself.
Since the Joule effect, caused by the current flowing along the longitudinal axis of the beam, induces a temperature rise that may get close to values leading to buckling instability (while pull-in is avoided due to the small amplitude of the vibrations of the flexure), the characteristics of the structural vibrations turn out to be highly affected by micromechanical imperfection-like fluctuations of stiffness and geometry of the moving polysilicon film.
A statistical analysis has been developed, in order to assess the sensitivity of the sensor response to uncertainties in the values of the over-etch depth and of the morphology-driven overall Young's modulus of the polycrystalline silicon. A Monte Carlo investigation has been carried out, wherein the over-etch depth has been sampled out of a micro-fabrication tailored normal distribution, and the Young's modulus has been sampled out of a size-dependent lognormal distribution obtained by best fitting the outcomes of a stochastic homogenization on samples (or statistical volume elements) of the polycrystalline film.
A sensitivity analysis has shown that the amplitude of oscillations of the resonant beam of this Lorentz force MEMS magnetometer is marginally affected by the elasticity of polysilicon, while it strongly depends on the over-etch depth. Similar results have been arrived at with the Monte Carlo simulations, which have highlighted that the scattering in the beam response and in the mechanical sensitivity of the device basically follow the trend induced by the over-etch and, therefore, by its geometry rather than by its elastic properties.
In the present analysis some simplifying assumptions have been introduced, in order to obtain the handled single degree of freedom reduced-order model. First, as motivated in the end by the far different effects of elastic and geometrical features on the vibrations amplitude, the statistical distributions of Young's modulus and over-etch depth have been assumed independent. This is clearly an approximation, as it has been remarked here above that the lognormal function describing the statistics of the polysilicon Young's modulus is size-dependent and, therefore, etch-dependent. Second, it has been assumed that the polysilicon Young's modulus and the over-etch depth do not vary along the beam axis, hence they are homogeneous in space. This is a further approximation since, in reality, the morphology of the film stochastically varies from point to point, and the over-etch depth is not uniform, even inside a single die. Third, the purposely developed reduced-order model has been based on one deformation mode for the structure, so that internal resonances and additional higher-order effects turned out to be all disregarded. Fourth, the values of the silicon parameters handled in the analysis have been assumed constant, though they can depend on the concentration of dopants like phosphorous and boron. All these effects are going to be taken into account in future developments of the present work.
Author Contributions: The authors contributed equally to this work.
Funding: Financial support provided by STMicroelectronics through project MaRe (Material Reliability) is gratefully acknowledged.