Nonlinear stochastic dynamics of an array of coupled micromechanical oscillators

The stochastic response of a multi‐degree‐of‐freedom nonlinear dynamical system is determined based on the recently developed Wiener path integral (WPI) technique. The system can be construed as a representative model of electrostatically coupled arrays of micromechanical oscillators, and relates to an experiment performed by Buks and Roukes. Compared to alternative modeling and solution treatments in the literature, the paper exhibits the following novelties. First, typically adopted linear, or higher‐order polynomial, approximations of the nonlinear electrostatic forces are circumvented. Second, for the first time, stochastic modeling is employed by considering a random excitation component representing the effect of diverse noise sources on the system dynamics. Third, the resulting high‐dimensional, nonlinear system of coupled stochastic differential equations governing the dynamics of the micromechanical array is solved based on the WPI technique for determining the response joint probability density function. Comparisons with pertinent Monte Carlo simulation data demonstrate a quite high degree of accuracy and computational efficiency exhibited by the WPI technique. Further, it is shown that the proposed model can capture, at least in a qualitative manner, the salient aspects of the frequency domain response of the associated experimental setup.

mechanisms relates to potential enhancement of the nanomechanical system detection sensitivity; see, for instance, the pioneering experimental work by Buks and Roukes. 8 Overall, it becomes clear that optimization of the design and enhancement of the detection capabilities of nanomechanical systems and devices dictate, first, modeling the micro/nanooscillator as a nonlinear multi-degree-of-freedom (multi-DOF) dynamical system subject to stochastic excitation. Second, potent uncertainty propagation methodologies are required for solving the governing equations of motion and for determining the system stochastic response. In other words, a high-dimensional system of coupled nonlinear stochastic differential equations needs to be solved accurately and in a computationally efficient manner for determining response statistics. In this regard, the state-of-the-art solution techniques in stochastic engineering dynamics can be broadly divided into two categories: first, those that exhibit a high degree of accuracy, but the associated computational cost becomes prohibitive with an increasing number of stochastic dimensions, and second, those that can readily treat high-dimensional systems, but provide reliable estimates for low-order response statistics only. The interested reader is directed to indicative standard books 9-11 for a broad perspective. Clearly, the development of versatile solution techniques that exhibit both high accuracy and low computational cost is critical for efficient stochastic response analysis of highdimensional systems of coupled nonlinear micro/nano-oscillators.
One of the promising solution techniques, recently pioneered in the field of engineering mechanics by Kougioumtzoglou and coworkers, 12,13 relates to the concept of the Wiener path integral (WPI). 14,15 According to the WPI technique, 16 the system response joint transition probability density function (PDF) is expressed as a functional integral over the space of all possible paths satisfying the initial and final conditions in time. Next, employing a functional integral series expansion, the contribution only of the first term is typically considered, pertaining to the path with the maximum probability of occurrence. This is referred to in the literature as the most probable path and corresponds to an extremum of the functional integrand. In this regard, the most probable path, which is used for determining the system response joint transition PDF approximately, is computed by solving a functional minimization problem that takes the form of a deterministic boundary value problem. 17 It is remarked that the WPI technique is capable of treating systems exhibiting diverse nonlinear/ hysteretic behaviors 18,19 and subjected to non-white and non-Gaussian stochastic excitations. 20 Further, it was shown by Psaros et al. 21 and by Katsidoniotaki et al. 22 that the associated computational cost can be reduced drastically by using sparse representations for the system response PDF in conjunction with compressive sampling concepts and tools. 23 Remarkably, there are only a few papers in the literature pertaining to stochastic modeling and analysis of nonlinear micro/ nano-oscillators. The vast majority of these research efforts relate to low-dimensional (typically single-DOF) systems, for which an analytical or numerical solution treatment is tractable. [24][25][26][27][28] The few papers referring to large arrays of coupled micro/nano-beams modeled as high-dimensional multi-DOF systems rely on significant simplifications and approximations that reduce, unavoidably, the accuracy degree of the stochastic response estimates. Indicatively, a moments equations solution approach was used by Ramakrishnan and Balachandran 29 for determining the stochastic response of an array of microcantilevers under the assumption of weak coupling.
Note, however, that the approximations inherent in the standard moments equations solution scheme 9 render it capable of yielding relatively accurate estimates only for the system response first-and second-order statistics (i.e., mean vector and covariance matrix). In this context, it can be argued that developing efficient and robust procedures for designing and optimizing nanomechanical systems requires the determination of the complete joint response PDF, or at least, of lower-dimensional joint PDFs corresponding to selected response coordinates of interest. In this regard, additional information relating, for example, to low-probability events (e.g., failures), and obtained by accurate estimation of the response PDF tails, can be integrated into the optimization methodology, leading to enhanced design. To this aim, the response PDF of a 100-DOF micromechanical system 30 was determined accurately and in a computationally efficient manner by Petromichelakis and Kougioumtzoglou. 31 This was done by developing a variational formulation of the WPI technique with mixed fixed/free boundary conditions that renders the computational cost independent from the total number of stochastic dimensions.
In this paper, attention is directed to an experiment performed

| Micromechanical system equations of motion
In the ensuing analysis, the focus is on the experimental setup by Buks and Roukes, 8 which can be construed as a representative case of electrostatically coupled micro/nano-resonators. In particular, the experiment pertained to a fabricated array of 67 doubly clamped microbeams that were parametrically excited by applying a timevarying voltage. In the same paper, the authors proposed a linear model to represent the dynamics of the coupled system of microbeams, that is, To address the inadequacy of the linear model of Equation (1) to capture satisfactorily the rich frequency content of the system response, Lifshitz and Cross 32 proposed a model comprising stiffness and damping nonlinearities of the polynomial kind, that is, and In Equation (2), ϵ 1 and ϵ 2 are parameters controlling the magnitudes of the stiffness and damping nonlinearities, respectively, and c 0 is a linear damping coefficient. Applying a perturbation solution treatment to the nonlinear Equation (2), it was shown 32 that the model can predict, qualitatively, the salient aspects of the system response in the frequency domain. Note, however, that the electrostatic forces were approximated as linear, in a similar manner as in Equation (1). In fact, this simplification was justified 32 by arguing that the effect of the elastic cubic nonlinearities g x ( ) 1 is significantly stronger relative to the electrostatic force nonlinearities, particularly for the case of small beam thickness compared to the gap between adjacent beams.
To enhance the accuracy of the model of Equation (2) and to also account for more general cases of microbeam arrays, where the gap between adjacent beams is comparable to their thickness (see Figure 1 for an indicative schematic representation), Equation (2) becomes Clearly, the electrostatic forces in Equations (5) and (6) are modeled as nonlinear, exhibiting a decaying behavior with increasing relative distance between the beams. In fact, Zhu et al. 33,34 considered an approximation of Equation (6) by utilizing the first two terms of a Taylor expansion of the nonlinear electrostatic forces.
It was shown, based on parametric resonance analysis, that the system response showed a more complex frequency content compared to the model of Equation (2). In passing, note that the linear assumption for the electrostatic forces in the models of Equation (1) and Equation (2) is equivalent to considering the first term only of a Taylor expansion of the nonlinear electrostatic forces.
In all previous models, the impact of stochastic noise sources 6,7 on the system modeling was ignored and the response analyses were purely deterministic. In this paper, a stochastic version of the model of Equation (5) where and . n n n n n n n n n n n n n N Further, the excitation vector w t ( ) represents a Gaussian , where t t , l l+1 are two arbitrary time instants and S w denotes an N N × power spectrum matrix given by where It is remarked that analytical evaluation of the functional integral of Equation (13) is, in general, an impossible task. Thus, researchers resort, routinely, to the most probable path approximation; that is, only the path x t ( ) c is used for the approximate evaluation of Equation (13) for which the stochastic action According to calculus of variations, 35 an extremal of can be determined by enforcing the necessary condition that the first variation equals zero, that is, δ = 0. Next, using a Taylor-type expansion for and integrating by parts, the condition δ = 0 for computing the most probable path x t ( ) c . Furthermore, substitut- into Equation (13), a specific point of the joint response transition PDF is determined approximately as where C is a normalization constant.
Clearly, in general, the boundary value problem of Equations (16) and (17) is not amenable to an analytical solution treatment, and therefore, numerical schemes are required. In this regard, adopting a brute-force solution approach, for a specific time instant t f , the values of the joint response PDF are computed based on Equation (18)  In this regard, the path integral representation of Equation (13) becomes Finally, solving Equations (16) and (20) numerically yields the most probable path  t x( ), and a specific point of the (lowerdimensional) joint response PDF is obtained as It is clearly seen that the WPI technique variational formulation with mixed fixed/free boundaries can reduce the associated computational cost drastically by determining directly any lowerdimensional joint response PDF. This capability is particularly advantageous for problems where the interest lies in few p q ( + ) specific DOF whose stochastic response is critical for the design and optimization of the N 2 -DOF system. In this regard, the L N 2 boundary value problems required to be solved by the standard formulation of the technique decrease to L p q + problems only, where L 30 < < 50 is a reasonable range of values for various diverse engineering applications. [16][17][18][19] Note that for small values of p q + relative to N 2 , it has been shown 31 that the technique becomes orders of magnitude more efficient than an alternative MCS solution treatment.

| NUMERICAL EXAMPLE
In agreement with the experiment carried out by Buks and Roukes, 8 consider a 67-DOF version of the system of Equation (7) 36 Second, a standard Runge-Kutta numerical integration scheme is used for solving Equation (7) and for determining response realizations. Finally, the system response PDF is estimated by performing a statistical analysis on the ensemble of the response realizations. Note that due to the specific type of nonlinearities in Equation (11), and in particular because of the form of the electrostatic forces, it is possible for the effective stiffness of the system to become negative. 37 This is followed, typically, by an unbounded response behavior. Nevertheless, in the numerical example, such an event has practically zero probability of occurrence for the selected parameters. In fact, for relatively small values of the excitation power spectrum magnitude S 0 , the restoring force in Equation (7) acts in the direction opposite to the displacement, and the system response exhibits a bounded behavior.
Next, to demonstrate the efficiency and reliability of the WPI technique presented in Section 2.2, only the final displacement x 10 and velocity ẋ 10 corresponding to the arbitrarily chosen 10th DOF are considered fixed; thus, p q + = 2, and the value L = 31 is used. This yields 31 2 boundary value problems to be solved for evaluating the joint PDF p x x ( ,˙) 10 10 at a given time instant. The WPI-based joint PDF p x x ( ,˙) 10 10 is plotted in Figure 2A as it evolves with time, and is compared in Figure 2B with MCS-based estimates (30 000 realizations). Obviously, the WPI technique exhibits a quite high degree of accuracy. This is further corroborated in Figure 3, where the WPIbased joint PDF p x x ( ,˙) 10 10 is compared with MCS-based estimates at indicative time instants. In Figure 4, the marginal displacement and velocity PDFs, that is, p x ( ) 10 and p x (˙) 10 , are plotted for time instants Lastly, the capability of the model of Equation (7) to reproduce, at least in a qualitative manner, the frequency domain response of the experimental setup by Buks and Roukes 8 is demonstrated next.
Specifically, for a given value of V dc , the frequency content of the response is estimated by applying Fourier transform to an arbitrary KATSIDONIOTAKI ET AL. | 7 response displacement realization corresponding to x 10 . This is produced by solving the system of Equation (7) excited by a randomly generated white noise realization. The results are plotted in Figure 5, which shows quite similar characteristics to Figure 5 in Buks and Roukes. 8 Clearly, the stochastic nonlinear model of Equation (7) proposed herein is capable of capturing, to a large extent, the salient aspects of the rich frequency content of the system response.

| CONCLUDING REMARKS
In this paper, an experimental setup by Buks