Some singular equations modeling MEMS

TL;DR

This survey reviews recent mathematical modeling efforts for MEMS, focusing on complex models with moving boundaries and their qualitative properties, highlighting recent advances in analysis of these sophisticated systems.

Contribution

It provides an overview of recent progress in analyzing complex MEMS models with moving boundaries, emphasizing the transition from simplified to more realistic models.

Findings

Recent analysis of complex MEMS models with moving boundaries

Progress in understanding qualitative properties of MEMS equations

Survey of recent mathematical techniques applied to MEMS modeling

Abstract

In the past fifteen years mathematical models for microelectromechanical systems (MEMS) have been the subject of several studies, in particular due to the interesting qualitative properties they feature. Still most research is devoted to an illustrative, but simplified model which is deduced from a more complex model when the aspect ratio of the device vanishes, the so-called vanishing (or small) aspect ratio model. The analysis of the aforementioned complex model involving a moving boundary has started only recently and an outlook of the results obtained so far in this direction is provided in this survey.

Equations

$${\epsilon }^2{\partial }_x^2\psi +{\epsilon }^2{\partial }_y^2\psi +{\partial }_z^2\psi =0,(x,y,z)\in \Omega (u(t)),t>0,$$

$$\Omega (u(t)):=\left\{(x,y,z)\in D\times (-1,\infty ):-1<z<u(t,x,y)\right\}$$

$$\psi (t,x,y,z)=\frac{1+z}{1+u(t,x,y)},(x,y,z)\in \partial \Omega (u(t)),t>0,$$

$$\begin{array}{rl}{\gamma }^2{\partial }_t^{2}u& +{\partial }_{t}u+\beta {\Delta }^{2}u-\tau \Delta u\\ & =-\lambda \left({\epsilon }^2|{\nabla }^{\prime }\psi (t,x,y,u(t,x,y)){|}^2+|{\partial }_z\psi (t,x,y,u(t,x,y)){|}^2\right)\end{array}$$

$$u=\beta {\partial }_{\nu }u=0\text{ on }\partial D,t>0,$$

$$u(0,x,y)=\gamma {\partial }_{t}u(0,x,y)=0,(x,y)\in D.$$

$$u=\beta \Delta u-\beta (1-\sigma )\kappa {\partial }_{\nu }u=0\text{ on }\partial D,t>0,$$

$$u=\beta \Delta u=0\text{ on }\partial D,t>0.$$

$${\psi }_0(t,x,y,z)=\frac{1+z}{1+u_0(t,x,y)}$$

$$\begin{array}{r}{\gamma }^2{\partial }_t^2{u}_0+{\partial }_t{u}_0+\beta {\Delta }^2{u}_0-\tau \Delta u_0=-\frac{\lambda }{(1+u_0{)}^2}\end{array}$$

$${\lambda }_0^{stat}={\lambda }_{0,0}^{dyn},$$

$${\lambda }_{0,\gamma }^{dyn}<{\lambda }_0^{stat},\gamma >0.$$

$${\lambda }_{\epsilon ,\gamma }^{dyn}<{\lambda }_{\epsilon ,0}^{dyn}={\lambda }_{\epsilon }^{stat},\epsilon >0,\gamma >0.$$

$$\hat{\Omega }(\hat{u}):=\left\{(\hat{x},\hat{y},\hat{z});(\hat{x},\hat{y})\in \hat{D},-H<\hat{z}<\hat{u}(\hat{x},\hat{y})\right\}$$

$$\Delta \hat{\psi }=0\text{in}\hat{\Omega }(\hat{u})$$

$$\hat{\psi }(\hat{x},\hat{y},-H)=0,\hat{\psi }(\hat{x},\hat{y},\hat{u}(\hat{x},\hat{y}))=V,(\hat{x},\hat{y})\in \hat{D},$$

$$\hat{\psi }(\hat{x},\hat{y},\hat{z})=V(1+\hat{z}),(\hat{x},\hat{y})\in \partial \hat{D}.$$

$$\hat{\mathcal{E}}(\hat{u})={\hat{\mathcal{E}}}_b(\hat{u})+{\hat{\mathcal{E}}}_s(\hat{u})+{\hat{\mathcal{E}}}_e(\hat{u}).$$

$${\hat{\mathcal{E}}}_b(\hat{u})=\frac{Yb{h}^3}{24}{\int }_{\hat{D}}{\left|\Delta \hat{u}(\hat{x},\hat{y})\right|}^2\mathrm{d}(\hat{x},\hat{y}),$$

$${\hat{\mathcal{E}}}_s(\hat{u})=\frac{T}{2}{\int }_{\hat{D}}{\left|\nabla \hat{u}(\hat{x},\hat{y})\right|}^2\mathrm{d}(\hat{x},\hat{y})+\frac{bhY}{8L}{\left({\int }_{\hat{D}}{\left|\nabla \hat{u}(\hat{x},\hat{y})\right|}^2\mathrm{d}(\hat{x},\hat{y})\right)}^2.$$

$${\hat{\mathcal{E}}}_e(\hat{u})=-\frac{{ϵ}_0}{2}{\int }_{\hat{\Omega }(\hat{u})}|\nabla \hat{\psi }(\hat{x},\hat{y},\hat{z}){|}^2\mathrm{d}(\hat{x},\hat{y},\hat{z})$$

$$x=\frac{\hat{x}}{L},y=\frac{\hat{y}}{L},z=\frac{\hat{z}}{H},u=\frac{\hat{u}}{H},\psi =\frac{\hat{\psi }}{V},$$

$$\begin{split}\mathcal{E}(u)&=\frac{Ybh^{3}\varepsilon^{2}}{24}\int_{D}\left|\Delta u(x,y)\right|^{2}\,\mathrm{d}(x,y)+\frac{TL^{2}\varepsilon^{2}}{2}\int_{D}\left|\nabla u(x,y)\right|^{2}\,\mathrm{d}(x,y)\\ &\quad+\frac{bhYL^{3}\varepsilon^{4}}{8}\left(\int_{D}\left|\nabla u(x,y)\right|^{2}\,\mathrm{d}(x,y)\right)^{2}\\ &\quad-\frac{\epsilon_{0}V^{2}L}{2\varepsilon}\int_{\Omega(u)}\big{(}\varepsilon^{2}|\partial_{x}\psi(x,y,z)|^{2}+\varepsilon^{2}|\partial_{y}\psi(x,y,z)|^{2}+|\partial_{z}\psi(x,y,z)|^{2}\big{)}\,\mathrm{d}(x,y,z)\end{split}$$

$$\Omega (u):=\left\{(x,y,z)\in D\times (-1,\infty ):-1<z<u(x,y)\right\}.$$

$${\epsilon }^2{\partial }_x^2\psi +{\epsilon }^2{\partial }_y^2\psi +{\partial }_z^2\psi =0,(x,y,z)\in \Omega (u),$$

$$\psi (x,y,z)=\frac{1+z}{1+u(x,y)},(x,y,z)\in \partial \Omega (u),$$

$$\begin{split}0&=-\varepsilon^{2}\,\beta\,\Delta^{2}u+\varepsilon^{2}\left(\tau+a\|\nabla u\|_{2}^{2}\right)\Delta u\\ &\quad-\varepsilon^{2}\lambda\big{(}\varepsilon^{2}|\nabla^{\prime}\psi(x,y,u(x,y))|^{2}+|\partial_{z}\psi(x,y,u(x,y))|^{2}\big{)}\end{split}$$

$$\beta :=\frac{Yb{h}^3}{12},\tau :=T{L}^2,a=a(\epsilon ):=\frac{bhY{L}^3{\epsilon }^2}{2},\lambda =\lambda (\epsilon ):=\frac{{ϵ}_0{V}^{2}L}{2{\epsilon }^3}.$$

$$\begin{split}\gamma^{2}\,\partial_{t}^{2}u+\partial_{t}u\,&=\,-\beta\,\Delta^{2}u+\left(\tau+a\|\nabla u\|_{2}^{2}\right)\Delta u\\ &\quad-\lambda\big{(}\varepsilon^{2}|\nabla^{\prime}\psi(t,x,y,u(t,x,y))|^{2}+|\partial_{z}\psi(t,x,y,u(t,x,y))|^{2}\big{)}\end{split}$$

$$D=(-1,1)\times (-b/L,b/L)\subset {\mathbb{R}}^2,$$

Full text

Partially supported by the French-German PROCOPE project 30718ZG

1. 1 Introduction

2. 1.1 The Free Boundary Problem for an Idealized MEMS

3. 1.2 The Vanishing Aspect Ratio Model for an Idealized MEMS

4. 1.3 A Bird’s Eye View on MEMS Dynamics

5. 1.4 Outline

6. 2 Derivation of the Governing Equations

7. 2.1 Linear Elasticity Model

8. 2.2 A One-Dimensional Linear Elasticity Model

9. 2.3 Nonlinear Elasticity Model

10. 2.4 Other Extensions

11. 3 Tool Box

12. 3.1 The Electrostatic Potential

13. 3.2 The Electrostatic Energy

14. 3.3 The Total Energy

15. 3.4 Maximum Principles

16. 4 Stationary Problem

17. 4.1 Stationary Solutions to the Free Boundary Problem

18. 4.1.1 Non-Existence

19. 4.1.2 Existence

20. 4.1.3 Multiplicity

21. 4.1.4 Vanishing Aspect Ratio Limit

22. 4.2 Stationary Solutions to the Vanishing Aspect Ratio Model

23. 5 Evolution Problem

24. 5.1 Local Existence

25. 5.2 Refined Criterion for Global Existence

26. 5.3 Global Existence

27. 5.4 Finite Time Singularity

28. 5.5 The Damping Dominated Limit $\gamma\rightarrow 0$

29. 5.6 The Vanishing Aspect Ratio Limit $\varepsilon\rightarrow 0$

30. 5.7 Dynamics in the Vanishing Aspect Ratio Model

31. 6 The Three-Dimensional Problem

32. 7 Summary and Open Problems

1. Introduction

Feynman anticipated and popularized the growing need for micro- and nanostructures with his famous lecture entitled “There is plenty of room at the bottom” [39] on the occasion of the American Physical Society’s annual meeting in 1959. Less than a decade later, Nathanson et al. [119] produced in 1967 the resonant gate transistor: the first batch-fabricated microelectromechanical system (MEMS) converting an input electrical signal into a mechanical force by utilizing electrostatic attraction. Since then, integrated devices on the scale of $1-100\mu m$, combining mechanical and electrical components, have become a well-established technology for everyday electronic products with applications in various domains, including automotive domain (e.g. accelerometers in airbag deployment systems, fuel or air pressure sensors, rollover detection, headlight leveling), commercial domain (e.g. inertia sensors in smartphones, inkjet printer heads), biotechnology (e.g. sensors for infusion pumps), industrial domain (e.g. optical switches, microvalves, micropumps, capacitors, gyroscopes), and many more [36, 139, 128, 75, 74]. With the increasing demand for miniaturized sensors and actuators of sophisticated functionalities, MEMS is regarded as one of the most promising technologies of this century. Needless to say that there are many challenges associated with this advancing miniaturization, both from a technological and theoretical point of view.

Electrostatically actuated MEMS devices suffer from an ubiquitous pull-in instability which limits their effectiveness. This phenomenon has been identified and described already in the seminal work of Nathanson et al. [119] for the simple situation of a lumped mass and spring system with two parallel plates at different potentials as depicted in Figure 1.

In this instability, when applied voltages are increased beyond a certain critical value, the top plate “pulls in” in the sense that the plates do no longer remain separate, restricting the range of stable operating of the device. Whether or not the touchdown of the upper plate on the lower one is a sought-for effect, key issues are not only to determine the range of applied voltages for which this phenomenon takes place but also the state of the device when it does not (for instance, it might be of interest to optimize the gap between the two plates). For the design of reliable high-performance devices, a detailed understanding of stable operating regimes and the occurrence of pull-in instabilities is thus of utmost importance.

Even in greatly simplified and idealized physical situations, the governing equations are mathematically challenging. So not only the above mentioned applications are of true interest, but also the mathematical investigation of the equations which involves various tools from many different fields from nonlinear analysis. The focus of this article is on a mathematical model for an idealized electrostatically actuated device in its “original” form of a free boundary problem [128, Section 7.4]. Maybe due to the many mathematical difficulties that come along with it, this free boundary problem was merely the source of simplified models and the research has focused so far mostly on a version thereof where the ratio height/length of the device is taken to vanish, thereby allowing one to compute the moving domain and to reduce the original free boundary problem to a single partial differential equation. While we only provide a short overview on this simplified model, we shall introduce the reader to some recent advances in the mathematical analysis of the free boundary problem.

1.1. The Free Boundary Problem for an Idealized MEMS

An idealized electrostatically actuated MEMS device which features the above mentioned instability is schematically drawn in Figure 2.

It consists of a rigid conducting ground plate above which an elastic plate (or membrane), coated with a thin layer of dielectric material and clamped on its boundary, is suspended. Holding the ground plate and the deformable plate at different electric potentials induces a Coulomb force across the device resulting in a deformation of the elastic plate and thus in a change in geometry, while electrostatic energy is converted into mechanical energy. There is thus a competition between attractive electrostatic forces and restoring elastic forces acting on the upper plate. After a suitable rescaling, the elastic plate’s position at rest is a two-dimensional domain $D$ located at $z=0$ while the fixed ground plate of the same shape is positioned at $z=-1$. The electrostatic potential on the ground plate is assumed to be zero while it takes a positive value $V>0$ on the elastic plate. The state of the device is fully described by the deformation $u=u(t,x,y)$ of the elastic plate in $z$-direction and the (dimensionless) electrostatic potential $\psi=\psi(t,x,y,z)$ between the plates. Introducing the aspect ratio $\varepsilon$ of the device (i.e. the ratio height/length), the electrostatic potential satisfies the rescaled Laplace equation

[TABLE]

in the time-varying cylinder

[TABLE]

between the rigid ground plate and the deflected elastic plate. Scaling the potential difference to $V=1$, the boundary conditions for $\psi$ are taken to be

[TABLE]

that is, $\psi=1$ on the elastic plate and $\psi=0$ on the ground plate. The dynamics of the elastic plate deformation is governed by a damped wave/plate equation

[TABLE]

for time $t>0$ and position $(x,y)\in D$ with clamped boundary conditions

[TABLE]

see, e.g., [128, 75, 36, 139] and Section 2.1. In (1.3) we use the notation $\nabla^{\prime}:=(\partial_{x},\partial_{y})$ for the gradient with respect to the horizontal directions. A commonly used initial state is that the elastic plate is at rest, that is,

[TABLE]

The fourth-order term $\beta\Delta^{2}u$ in equation (1.3) reflects plate bending while the linear second-order term $\tau\Delta u$ with $\tau\geq 0$ accounts for external stretching. Taking $\beta>0$ (with $\tau\geq 0$) reflects the case of an elastic plate while $\beta=0$ (with $\tau>0$) rather refers to an elastic membrane. The right-hand side of (1.3) is due to the electrostatic forces exerted on the elastic plate. The parameter $\lambda>0$, being proportional to the square of the applied voltage difference and the device’s aspect ratio $\varepsilon>0$ are key quantities in the model and its analysis. In the sequel we focus on the situation of an elastic plate clamped at its boundary as modeled in (1.4) though other boundary conditions are possible. For instance, assuming that the domain $D$ is regular enough and denoting the curvature of $\partial D$ (taken to be positive for convex domains) by $\kappa$, the hinged boundary conditions (also referred to as Steklov boundary conditions)

[TABLE]

are derived from the description of the deformation of a device which is ideally hinged along all its boundary so that it is free to rotate and cannot support a bending moment [135, 46], the parameter $\sigma$ being the Poisson ratio which in general ranges in $[0,1/2)$. A variation hereof, often met in the literature, are the pinned boundary conditions

[TABLE]

From a mathematical point of view, the aforementioned pull-in phenomenon occurring when the elastic plate contacts the ground plate corresponds to the deflection $u$ attaining the value $-1$ somewhere in $D$ at some time $t>0$ and is also referred to as the touchdown phenomenon. Equations (1.1)-(1.2) are then no longer meaningful since the set $\Omega(u(t))$ becomes disconnected and solutions cease to exist. Moreover, if $(x_{0},y_{0})$ denotes a touchdown point in $D$, the vertical derivative $\partial_{z}\psi(t,x_{0},y_{0},u(t,x_{0},y_{0}))$ appearing on the right-hand side of (1.3) develops a singularity as $\psi=1$ along $z=u(t,x,y)$ while $\psi=0$ along $z=-1$. In some sense, the parameter $\lambda$ tunes this singularity and consequently plays a crucial rôle in the pull-in instability which is closely related to the existence of stationary (i.e. time-independent) solutions and the question of global existence of solutions to (1.1)-(1.5).

The elliptic equation (1.1) is a free boundary problem posed in the non-smooth cylinder $\Omega(u(t))$ which changes its geometry as its upper surface – that is, the elastic plate being described at instant $t$ by $u(t)$ – evolves in time. The inherently nonlinear character of this free boundary problem induces another mathematical challenge to the analysis of (1.1)-(1.5) in the form of an intricate coupling of the two unknowns $\psi$ and $u$ : Any solution $\psi$ to (1.1)-(1.2) and accordingly the right-hand side of the wave equation (1.3) are nonlocal and nonlinear functions of $u$ and likely to enjoy restricted regularity. It is thus certainly not obvious how the function $\psi$ reacts on changes with respect to $u$.

In this connection, let us point out that the system (1.1)-(1.4) has a gradient structure [5, 9, 11, 36, 95]. More precisely, there is an energy functional $\mathcal{E}$ defined for suitable deflections $u$ such that stationary solutions to (1.1)-(1.4) correspond to critical points of $\mathcal{E}$ while it decays along any trajectory of the evolution problem. However, as we shall see below, the use of this variational structure is delicate due to the non-coercivity of the functional. In fact, the energy $\mathcal{E}(u)$ diverges to $-\infty$ if the minimum of $u$ approaches $-1$. Nonetheless, this structure provides a simple way of deriving the model (Section 2) and is also helpful in excluding the occurrence of some singularities.

1.2. The Vanishing Aspect Ratio Model for an Idealized MEMS

The geometrical difficulties in the governing equations can be hidden when the aspect ratio $\varepsilon$ is negligibly small. Indeed, setting formally $\varepsilon=0$ in (1.1), we readily infer from (1.1)-(1.2) that the electrostatic potential $\psi_{0}$ is given by the explicit formula

[TABLE]

as a function of the yet to be determined deflection $u_{0}$ according to

[TABLE]

for $t>0$, $(x,y)\in D$, supplemented with the boundary condition (1.4) and initial conditions (1.5). The original free boundary problem is thus explicitly solved and equation (1.3) is reduced to a single semilinear equation (1.9) known as the vanishing aspect ratio equation. The nonlinear reaction term on its right-hand side is no longer nonlocal but still singular when $u_{0}=-1$, which corresponds to the touchdown phenomenon already alluded to. This singularity is also referred to as quenching in the literature concerning parabolic equations with singular right-hand sides, see [40, 99, 100] and the references therein. In addition, from the explicit form of the nonlinearity on the right-hand side of (1.9) utmost favorable properties can be directly read off, important and well exploited in the investigation of qualitative aspects of solutions: the reaction term of equation (1.9) is monotone, concave, and of zero order. These properties are not known and by no means obvious for the equation (1.3) involving the nonlocal reaction term on its right-hand side. Still, the singular equation (1.9) depending on the parameter $\lambda$ presents various mathematical challenges and exhibits many interesting phenomena which is the reason why the vast literature dedicated to vanishing aspect ratio models for MEMS is continuously growing.

1.3. A Bird’s Eye View on MEMS Dynamics

As we shall see more precisely later on, the qualitative behavior of solutions to both (1.3) and (1.9) significantly depend on whether or not some of the coefficients $\beta$, $\gamma$, and $\tau$ are assumed to be zero as well as on the shape and the dimension of the plate $D$. In particular, the order of the partial differential equation (fourth-order when $\beta>0$ versus second-order when $\beta=0$) and its character (hyperbolic when $\gamma>0$ versus parabolic when $\gamma=0$) play important rôles.

As already mentioned a common peculiarity of the dynamics of MEMS is the pull-in instability which is widely observed in real world applications and extensively discussed in engineering literature (see, e.g., [139, 36] for references). Mathematically it manifests in different ways:

• $\star$

In the stationary regime it is (presumably) revealed as a sharp stationary threshold $\lambda_{\varepsilon}^{stat}>0$ (depending on $\varepsilon$) for the voltage parameter $\lambda$ above which there is no stationary solution while below there is at least one, the exact number of steady-state solutions for a specific value of $\lambda\in(0,\lambda_{\varepsilon}^{stat}]$ depending heavily on the space dimension. This threshold value of the voltage parameter is often referred to as the (stationary) pull-in voltage in the literature.

• $\star$

For the evolution problem the pull-in instability corresponds to the occurrence of a finite time singularity, the touchdown phenomenon, in which the elastic plate, initially at rest, collapses onto the ground plate or, equivalently, the plate deflection $u$ reaches the value $-1$. Similarly as in the stationary case, it is expected that a dynamic threshold value $\lambda_{\varepsilon,\gamma}^{dyn}>0$ (depending on $\varepsilon$ and $\gamma$) exists for the dynamical regime such that touchdown takes place when $\lambda>\lambda_{\varepsilon,\gamma}^{dyn}$ while the deflection $u$ stays above $-1$ for all times when $\lambda<\lambda_{\varepsilon,\gamma}^{dyn}$. It is worth pointing out that, as observed in [80], this threshold value is likely to depend on the initial deformation of the elastic plate. We emphasize that throughout this paper “pull-in” in the dynamic case always refers to a touchdown of the elastic plate when it is initially at rest as in (1.5).

Neither the existence of sharp pull-in voltages $\lambda_{\varepsilon}^{stat}$ and $\lambda_{\varepsilon,\gamma}^{dyn}$ nor the relation between them are fully established in all cases, not even for the vanishing aspect ratio model (1.9) where $\varepsilon=0$. In fact, for the fourth-order case $\beta>0$ with clamped boundary conditions (1.4), the only available result appears to be that $\lambda_{0,\gamma}^{dyn}\leq\lambda_{0}^{stat}$ as shown in [94, 109]. Still when $\varepsilon=0$, more can be said regarding the second-order setting with $\beta=0$ which we describe now, see also Section 4.2. It is known [49] that

[TABLE]

while numerical results [41] and experimental investigations [134] indicate that

[TABLE]

This implies that the hyperbolic evolution can drive the device to a touchdown even though there is a stationary configuration. Thus, precise details and the exact values of pull-in voltages are of utmost importance from the viewpoint of applications.

We also point out that the order of the differential equation impacts not only on the structure of the stationary solutions but also on the dynamics of the touchdown behavior, at least for the vanishing aspect ratio model (1.9). Indeed, on the one hand, the fourth-order ($\beta>0$) stationary problem in a one- or two-dimensional domain $D$ is in many respects qualitatively the same as that of the second-order equation ($\beta=0$) in a one-dimensional interval $D$ (see [35] and Section 2.2). On the other hand, there is numerical evidence that touchdown occurs in one space dimension at a single point when $\beta=0$ and $\lambda>\lambda_{0}^{stat}=\lambda_{0,0}^{dyn}$, while it takes place at two different points when $\beta>0$ and $\lambda\gg\lambda_{0,0}^{dyn}$ [109].

A similar behavior as observed in the vanishing aspect ratio equation is expected for the free boundary problem (1.1)-(1.5) with $\varepsilon>0$, that is, that a stationary and dynamical pull-in voltage values $\lambda_{\varepsilon}^{stat}$ and $\lambda_{\varepsilon,\gamma}^{dyn}$ exist which may be different. Numerical evidence [43] implies that the values for the pull-in voltages indeed depend on both $\varepsilon$ and $\gamma$ and, furthermore, on the hyperbolic or parabolic character since

[TABLE]

1.4. Outline

In the following, we shall give more details regarding analytical results on the free boundary problem (1.1)-(1.5). In Section 2 we formally derive the equations as Euler-Lagrange equations of the underlying energy balance within linear (Section 2.1) and nonlinear (Section 2.3) elasticity theory. We also present the equations in reduced dimensions (Section 2.2) – which will be the main focus of this article – and briefly discuss various extensions of the model (Section 2.4), which, however, will not be pursued further in this article. In Section 3 we collect some useful tools needed for the analysis of (1.1)-(1.5). We first investigate the subproblem (1.1)-(1.2) for the electrostatic potential and show how to formulate the free boundary problem as a single equation for the elastic plate deflection (Section 3.1). We then discuss the rigorous proof that the free boundary problem has indeed a variational structure (Section 3.2) and finally present some results related to the fourth-order maximum principle (Section 3.4), which – when available – is a central tool for the analysis. Section 4 is dedicated to existence and non-existence results for stationary solutions and the related pull-in instability (Section 4.1). We also summarize some of the research on stationary solutions to the vanishing aspect ratio model (Section 4.2). The evolution problem is considered in Section 5 which includes local (Section 5.1) and (non-) global existence results (Section 5.3) with regard to the pull-in phenomenon. This section also contains a brief discussion on the damping dominated limit (Section 5.5) and the vanishing aspect ratio limit (Section 5.6). In Section 6 we give a short summary on the (original) three-dimensional case. Finally, in Section 7 we mention some open problems that could give raise to further investigations.

2. Derivation of the Governing Equations

The equations governing the dynamics of the elastic plate in an electrostatically actuated MEMS may be derived from its energy balance after a suitable rescaling as Euler-Lagrange equations [5, 9, 11, 36, 96, 128]. We shall first derive the equations presented in the introduction using linear elasticity theory and then introduce variants thereof under different modeling assumptions.

2.1. Linear Elasticity Model

In Figure 2, an idealized MEMS device is schematically drawn. We assume that the elastic plate, held at potential $V$ and clamped along its boundary, is of shape $\hat{D}\subset\mathbb{R}^{2}$ with characteristic dimension $2L$ and rest position at $\hat{z}=0$. The grounded fixed plate is also of shape $\hat{D}$ and located at $\hat{z}=-H$. We denote the vertical deflection111See e.g. [121] for a more sophisticated model including non-vertical large deflections based on the von Kármán equations. of the elastic plate at position $(\hat{x},\hat{y})\in\hat{D}$ and time $\hat{t}$ by $\hat{u}=\hat{u}(\hat{t},\hat{x},\hat{y})>-H$ and the electrostatic potential between the plates at position $(\hat{x},\hat{y},\hat{z})$ and time $\hat{t}$ by $\hat{\psi}=\hat{\psi}(\hat{t},\hat{x},\hat{y},\hat{z})$. We suppress the time variable $\hat{t}$ for the moment. The electrostatic potential $\hat{\psi}$ is harmonic in the region

[TABLE]

between the ground plate and the elastic plate, that is,

[TABLE]

and with given values on the plates

[TABLE]

extended continuously on the vertical sides of $\hat{\Omega}(\hat{u})$ by

[TABLE]

This choice is motivated by formula (1.8) for $\psi$ in the vanishing aspect ratio limit. The total energy generated by the deformation of the elastic plate is [5, 9, 11, 36, 128]

[TABLE]

It involves the bending energy $\hat{\mathcal{E}}_{b}$, for simplicity taking into account only the linearized curvature of the elastic plate, that is,

[TABLE]

where $Y$ is Young’s modulus, $h$ is the height of the plate and $b$ its width, the latter being assumed to be of the same order as the plate’s length $2L$. The stretching energy $\hat{\mathcal{E}}_{s}$ retains the linearized contribution due to the axial tension force with coefficient $T$ together with a self-stretching force coming from elongation of the plate in the case of (moderately) large oscillations. It is of the form [36]

[TABLE]

The electrostatic energy $\hat{\mathcal{E}}_{e}$ is given by

[TABLE]

with $\epsilon_{0}$ being the permittivity of free space. Note that $\hat{\psi}$ is the maximizer of the (negative) Dirichlet energy $-(\epsilon_{0}/2)\|\nabla\hat{\xi}\|_{L_{2}(\hat{\Omega}(\hat{u}))}^{2}$ among functions $\hat{\xi}\in H^{1}(\hat{\Omega}(\hat{u}))$ satisfying the prescribed boundary conditions on $\partial\hat{\Omega}(\hat{u})$. Introducing dimensionless variables

[TABLE]

and denoting the aspect ratio of the device by $\varepsilon=H/L$, the total energy in these variables is

[TABLE]

with $D:=\{(x,y)\,;\,(Lx,Ly)\in\hat{D}\}$ and

[TABLE]

The equilibrium configurations of the device are the critical points of the total energy, i.e. satisfy $\partial_{u}\mathcal{E}(u)=0$. The computation of this Fréchet derivative with respect to $u$ is quite involved since $\psi$ solves the rescaled equation

[TABLE]

in the $u$-dependent domain $\Omega(u)$ along with the boundary condition

[TABLE]

and hence $\psi=\psi_{u}$ depends in a nonlocal way on the deformation $u$. Interpreting the derivative of the electrostatic energy as the shape derivative of the Dirichlet integral of $\psi=\psi_{u}$ and using shape optimization arguments [69], the Euler-Lagrange equation is (see Section 3.2 and [95])

[TABLE]

for $(x,y)\in D$, where $\nabla^{\prime}=(\partial_{x},\partial_{y})$ and

[TABLE]

The dynamics of the time-dependent deflection $u=u(\hat{t},x,y)$ is derived by means of Newton’s second law. Letting $\rho$ and $\delta$ denote the mass density per unit volume of the elastic plate respectively its thickness, the sum over all forces equals $\rho\delta\partial_{\hat{t}}^{2}u$. The elastic and electrostatic forces are combined with a damping force of the form $-r\partial_{\hat{t}}u$ being linearly proportional to the velocity. Scaling time based on the strength of damping according to $t=\hat{t}\varepsilon^{2}/r$ and setting $\gamma:=\varepsilon\sqrt{\rho\delta}/r$, we derive for the dimensionless deflection $u$ the evolution equation

[TABLE]

for $t>0$ and $(x,y)\in D$, which is (1.3) when $a=0$, that is, when self-stretching forces due to elongation of the plate are neglected.

Remark 2.1**.**

An alternative, though formal, way to obtain the electrostatic force exerted on the elastic plate, which avoids the use of shape derivatives, is to treat the device as a parallel plate capacitor [128, 127, 35, 107] so that classical electrostatics can be applied, see e.g. [73].

2.2. A One-Dimensional Linear Elasticity Model

A common building block in MEMS devices are rectangular plates, i.e., after rescaling,

[TABLE]

where the elastic plate is only held fixed along the edges in $y$-direction while the edges in $x$-direction are free. This structure is depicted in Figure 3 and yields a one-dimensional equation for the plate deformation when assuming homogeneity in $y$-direction.

In this case, $u=u(t,x)$ and $\psi=\psi(t,x,z)$ are independent of the variable $y$ so that the equations (1.1)-(1.5) reduce to

[TABLE]

in the two-dimensional time-varying domain

[TABLE]

with $I:=(-1,1)$, coupled to the one-dimensional equation

[TABLE]

for $t>0$ and $x\in I$ in which the self-stretching force has been neglected (i.e. $a=0$), supplemented with initial and boundary conditions

[TABLE]

As we shall see later, the mathematical analysis in this simplified geometry is already quite involved and we shall restrict in this survey to the particular case of equations (2.1)-(2.5). In Section 6 we present extensions for the two-dimensional case (1.1)-(1.5). In principle, the same approach can be used in both cases. However, much weaker regularity is obtained in the two-dimensional setting and consequently, the analysis is less complete so far.

2.3. Nonlinear Elasticity Model

The previous linear elasticity model is restricted to small deformations neglecting curvature effects from the outset. The governing equation (1.3) (respectively (2.3)) for the plate deflection is thus semilinear as its highest order term $\beta\Delta^{2}u$ is linear. As pointed out in [9] it may well be, however, that retaining gradient terms affects quantitatively and qualitatively the occurrence of the pull-in instability and is thus important in applications. To account for curvature effects in the simpler one-dimensional setting without self-stretching forces (i.e. $a=0$), the bending energy $\mathcal{E}_{b}$ is replaced after rescaling by

[TABLE]

involving the curvature of the graph of $u$ and its arc length element, while the stretching energy $\mathcal{E}_{s}$ is proportional to the change of arc length and replaced by

[TABLE]

Analogously as before, the Euler-Lagrange equation of the total energy gives rise to the governing equation for the plate deflection. Consequently, taking curvature effects into account, equation (2.3) is replaced by

[TABLE]

with the quasilinear fourth-order operator

[TABLE]

The resulting problem (2.1), (2.2), (2.6), (2.4), (2.5) is studied in [96] and in the second-order case (i.e. $\beta=0$) in [31, 29].

2.4. Other Extensions

There are various variants and extensions of the governing equations for a MEMS device in order to incorporate additional physical effects. Except for [32, 105, 106], where the case of a non-uniform potential applied along the elastic plate is studied, and [83, 84, 85], where a model involving two elastic plates suspended one above each other is considered, the analytical investigations regarding extended models have been concerned so far almost exclusively with vanishing aspect ratio models like (1.9). For instance, spatial variations in the dielectric properties of the elastic plate which modify the strength of the Coulomb force may be included [127] or the device may be embedded in a circuit in an attempt to control the pull-in instability [131]. If a simple capacitive control scheme is used, the right-hand side of equation (1.9) is then replaced by a term

[TABLE]

where $\chi$ is a ratio of capacitances in the system and $1/f(x,y)$ is proportional to the local dielectric permittivity. Observe that, when $\chi>0$, the additional nonlocal term does not depend on the spatial variables and may in some way be included in the parameter $\lambda$.

Actually, the common vanishing aspect ratio approximation leading to (1.9) is not uniformly valid in the spatial variables as fringing fields – unaccounted for in the model (also with free boundary) – may cause corrections in regions close to the boundary of the domain. Corner-corrected models are derived in [126, 128, 130] and the right-hand side of equation (1.9) is replaced by

[TABLE]

On a nanoscale, when the device is below a critical size, van der Waals and Casimir forces may become important which represent attraction of two uncharged bodies and may lead to a spontaneous collapse of the system with zero applied voltage. Such a phenomenon is considered in [3, 53] for instance and leads to the addition of a term of the form

[TABLE]

on the right-hand side of (1.9) with $m=3$ (van der Waals forces) or $m=4$ (Casimir forces).

Somewhat in the opposite direction, the addition of a term $\lambda\eta^{m-2}(1+u_{0})^{-m}$, $m>2$, on the right-hand side of (1.9) is suggested in [110, 111, 108] with small parameter $\eta>0$ mimicking the effect of a thin insulating layer coating the ground plate and thereby preventing the touchdown phenomenon.

Finally, if a MEMS device is not fabricated to operate in a vacuumed capsule, the dynamics of the moving plate can be affected by the gas (or liquid) between the plates. This squeeze film phenomenon is taken into account by adding a pressure term on the right-hand side of (1.9) which itself is governed by the nonlinear Reynolds equation and depends on the deflection $u_{0}$. We refer to [36, 120, 2, 121, 75] for a detailed description of the model which has been investigated only numerically so far.

3. Tool Box

A promising solving strategy for (2.1)-(2.5) is a decoupling of the equations for the electrostatic potential $\psi$ and the plate deflection $u$. Focusing first on the subproblem (2.1)-(2.2) for the electrostatic potential with an (almost) arbitrarily given, but fixed deflection $u$ and investigating the dependence of its solution $\psi=\psi_{u}$ on $u$, is the basis to define

[TABLE]

which also reads

[TABLE]

owing to the constant Dirichlet boundary condition (2.2). The time dependence is omitted in (3.1) and (3.2) since time acts only as a parameter in the elliptic subproblem (2.1)-(2.2). One can then reformulate (2.1)-(2.5) as a single nonlocal problem

[TABLE]

only involving the plate deflection $u$. In the following we shall collect some auxiliary results needed for this approach. In Section 3.1 we first provide the reduction to (3.3) and state important properties of the function $g_{\varepsilon}$. In Section 3.2 we then set the stage for a variational formulation of (3.3). Existence of solutions to (3.3) and their qualitative aspects shall be studied later in Section 4 and Section 5.

3.1. The Electrostatic Potential

For an arbitrarily fixed deflection $u$, the subproblem (2.1)-(2.2) is an elliptic equation with Dirichlet boundary data. Thus only minimal regularity assumptions on $u$ – including its continuity along with the constraint $u(\pm 1)=0$ guaranteeing an open domain $\Omega(u)$ – are needed in order to ensure the existence of a unique weak solution $\psi_{u}$ to (2.1)-(2.2) in the Sobolev space $H^{1}(\Omega(u))$ by the well-known Lax-Milgram Theorem. However, the provided regularity for $\psi_{u}$ is far from being sufficient to give a meaning to the square of its gradient trace $|\nabla\psi_{u}(\cdot,u)|^{2}$ needed later to define $g_{\varepsilon}(u)$ in (3.1). Moreover, no information is so obtained on the dependence of this term on the given deflection $u$. For this reason, a beneficial approach is to restate the Laplace equation (2.1)-(2.2) for $\psi_{u}$ as a boundary value problem with zero Dirichlet data of the form

[TABLE]

for the transformed electrostatic potential

[TABLE]

in the fixed rectangle $\Omega:=I\times(0,1)$. The differential operator $-\mathcal{L}_{u}$ has coefficients depending on $u$ and its $x$-derivatives up to second-order and being singular at touchdown points $x\in I$ where $u(x)=-1$. However, as long as $u(x)>-1$ for all $x\in I$, the operator $-\mathcal{L}_{u}$ is elliptic. We thus restrict to deflections $u$ belonging to

[TABLE]

for some $\alpha>1/q$ and $\kappa\in(0,1)$, where $q\in[1,\infty]$. Here

[TABLE]

We stress here that, when $\beta=0$ and $\alpha>1/q$, the first-order boundary conditions are irrelevant and $W_{q,D}^{\alpha}(I)$ is simply the subspace of the classical Sobolev space $W_{q}^{\alpha}(I)$ consisting of those functions $v$ satisfying only $v(\pm 1)=0$. Observe that if $\alpha>1/q$, then $W_{q}^{\alpha}(I)$ is continuously embedded in $C(\bar{I})$ so that $\Omega(u)$ is a well-defined open set of $\mathbb{R}^{2}$. For simplicity we put

[TABLE]

Depending on the smoothness of the deflection $u$, we derive different regularity properties of the electrostatic potential $\psi_{u}$. Introducing the Sobolev spaces

[TABLE]

we can use elliptic regularity theory to show that, for sufficiently smooth deflections $u$, the operator $-\mathcal{L}_{u}$ is bounded and invertible from $H_{D}^{2}(\Omega)$ onto $L_{2}(\Omega)$ so that the transformed potential $\Phi_{u}$ solving (3.4) belongs to $H_{D}^{2}(\Omega)$ and depends smoothly on $u$. For less regular deflections $u$, the operator $-\mathcal{L}_{u}$ no longer maps $H_{D}^{2}(\Omega)$ into $L_{2}(\Omega)$. Still it is a bounded and invertible operator from $H_{D}^{1}(\Omega)$ onto $(H_{D}^{1}(\Omega))^{\prime}$. A careful analysis reveals that its inverse maps $L_{2}(\Omega)$ into the anisotropic space

[TABLE]

and $\Phi_{u}\in X(\Omega)$ is weakly sequentially continuous with respect to $u$. The so-derived key properties in terms of the electrostatic potential $\psi_{u}$ are collected in the following result.

Proposition 3.1** (Nonlocal Nonlinearity [30, 91]).**

1. (i)

For each $u\in S_{2}^{1}$ there is a unique weak solution $\psi_{u}\in H_{D}^{1}(\Omega(u))$ to (2.1)-(2.2). 2. (ii)

Assume that $\beta>0$ and consider $\alpha\in(3/2,2]$. If $u\in S_{2}^{\alpha}$, then $\psi_{u}\in H^{\alpha}(\Omega(u))$, and the mapping

[TABLE]

is continuous for each $\sigma\in[0,1/2)$. 3. (iii)

Let $\kappa\in(0,1)$, $q\in(2,\infty]$, and $\sigma\in[0,1/2)$. Then $g_{\varepsilon}:S_{q}^{2}(\kappa)\rightarrow H^{\sigma}(I)$ is globally Lipschitz continuous, bounded, and analytic.

The proofs of the statements (ii) and (iii) in Proposition 3.1 are of a completely different nature. While (iii) is a consequence of the regularity theory for elliptic equations with sufficiently smooth coefficients [88], the proof of (ii) is more involved and requires to exploit carefully the specific structure of the right-hand side $f_{u}$ of (3.4) and the vanishing of $\partial_{x}u(\pm 1)$ to derive estimates in $X(\Omega)$ for $\Phi_{u}$. In particular, it is restricted to the fourth-order case $\beta>0$. Note that $g_{\varepsilon}$ acts in statement (ii) as a mapping of (roughly) order one while in (iii) it is of order $3/2$, but with significant additional properties.

The particular geometry of $\Omega(u)$ with its corners apparently precludes a higher regularity than $H^{2}(\Omega(u))$ for the electrostatic potential $\psi_{u}$. This is a far-reaching technical obstruction for the analysis. Also, as pointed out in the introduction, convexity and monotonicity properties of the electrostatic potential or of the function $g_{\varepsilon}$ – obvious and important for the study of the vanishing ratio equation (1.9) – are not known in the present situation.

3.2. The Electrostatic Energy

A useful tool in the investigation of qualitative aspects for solutions to equation (3.3) is its gradient flow structure. In Section 2 we have seen that (the stationary) equation (3.3) is formally the Euler-Lagrange equation of the (rescaled) total energy

[TABLE]

involving the mechanical energy

[TABLE]

and the electrostatic energy

[TABLE]

where $\psi_{u}$ is the electrostatic potential associated with a given deflection $u$ and provided by Proposition 3.1. The stumbling block in this direction is obviously the computation of the Fréchet derivative with respect to $u$ of the elastic energy $\mathcal{E}_{e}$, which involves the implicit dependence of $\psi_{u}$ on $u$ via the domain $\Omega(u)$. Nevertheless, this difficulty can be overcome by interpreting the derivative as the shape derivative of the Dirichlet integral of $\psi_{u}$, which can be computed by shape optimization arguments [69] and a transformation to the fixed rectangle $\Omega$ introduced in the previous section. More precisely, the rigorous result is the following.

Proposition 3.2** (Derivative of the Electrostatic Energy [91]).**

The functional $\mathcal{E}_{e}$ belongs to $C(S_{2}^{1})$, and if $\beta>0$ and $\alpha\in(3/2,2]$, then $\mathcal{E}_{e}\in C^{1}(S_{2}^{\alpha})$ with Fréchet derivative $\partial_{u}\mathcal{E}_{e}(u)=-g_{\varepsilon}(u)$ for $u\in S_{2}^{\alpha}$.

According to Proposition 3.2, the Fréchet derivative of $\mathcal{E}_{e}$ in $S_{2}^{2}$ is non-negative when $\beta>0$, which indicates a monotonicity property of $\mathcal{E}_{e}$. The latter is actually valid in a more general setting and is somewhat reminiscent of that of the classical Dirichlet energy [69].

Proposition 3.3** (Monotonicity of the Electrostatic Energy [91]).**

Consider two functions $u_{1}$ and $u_{2}$ in $S_{2}^{1}$ such that $u_{1}\leq u_{2}$. Then $\mathcal{E}_{e}(u_{2})\leq\mathcal{E}_{e}(u_{1})$.

Unfortunately, a monotonicity (or concavity) property of the derivative of $\mathcal{E}_{e}$, that is, of the right-hand side $-g_{\varepsilon}$ of the evolution equation (3.3), is not known. This is in clear contrast to the vanishing aspect ratio equation (1.9), where such a monotonicity is obvious and of utmost importance in the analysis. We shall see later in Section 4 and Section 5, however, how the information provided by Proposition 3.2 and Proposition 3.3 on the variational structure of (3.3) can be used to gain insight into qualitative aspects. Still, some care is needed if one wishes to exploit this variational structure. Indeed, it is quite clear that $\mathcal{E}_{e}$ is unbounded from below on $S_{2}^{1}$ and quantified by the estimates

[TABLE]

To prove the first inequality in (3.7) it suffices to note that (2.2) and the Cauchy-Schwarz inequality imply

[TABLE]

The second one is a consequence of the variational inequality satisfied by the weak solution $\psi_{u}$ to (2.1)-(2.2) (see [91, Lemma 2.8] for details). Clearly, $\mathcal{E}_{e}(u)$ is not finite for any deflection $u\in C^{1}(\bar{I})$ satisfying $u\geq-1$ and $u(x_{0})=-1$ for some $x_{0}\in I$. It is worth pointing out that the vanishing aspect ratio $\varepsilon=0$ leads to equality on both sides in (3.7) and thus, the bounds on $\mathcal{E}_{e}$ when $\varepsilon>0$ seem to be nearly optimal. That $\mathcal{E}_{e}$ features a singularity for certain touchdown deflections becomes also apparent from the striking identity [91]

[TABLE]

which is derived from (1.1)-(1.2) by means of Green’s formula and indicates that the blowup of $\mathcal{E}_{e}(u)$ resulting from the touchdown phenomenon is accompanied by a singular behavior of $x\mapsto\partial_{z}\psi_{u}(x,u(x))$ as already mentioned.

3.3. The Total Energy

Proposition 3.2 implies that the stationary solutions to (3.3) are the critical points of the energy $\mathcal{E}=\mathcal{E}_{m}-\lambda\mathcal{E}_{e}$ defined in (3.6) and that (3.3) may be seen as a gradient flow associated with $\mathcal{E}$. In particular, $\mathcal{E}$ is decreasing along sufficiently smooth solutions to the dynamic problem (3.3).

Corollary 3.4** (Total Energy [95]).**

• (i)

Let $u\in S_{2}^{4}$. Then $u$ is a stationary solution to (3.3) if and only if it is a critical point of the total energy $\mathcal{E}$.

• (ii)

Let $u\in C([0,T],S_{2}^{2+\nu})\cap C^{1}([0,T],H_{D}^{\nu}(I))$ solve (3.3) for some $\nu>0$ and $T>0$. Then

[TABLE]

for $t\in[0,T]$.

Retrieving tractable information from this variational structure is, however, not obvious since the total energy $\mathcal{E}$ is clearly the sum of terms with different signs and the possible pull-in instability becomes manifest in the non-coercivity of the energy $\mathcal{E}$. Nevertheless, we shall use it later on to construct stationary solutions by a constrained variational approach, thereby bypassing the non-coercivity of the energy. Also, combining (3.7) and (3.9) provides regularity estimates on $u$ as soon as one knows that it is separated from $-1$.

3.4. Maximum Principles

When available, the celebrated maximum principles are central tools in the analysis of differential equations and guarantee, roughly speaking, non-negativity (or even positivity) of solutions provided the boundary and/or initial data are non-negative and the equations satisfy suitable properties. Thus, it is not surprising that the deepest results regarding MEMS equations are for situations where a maximum principle holds true. While such a tool is well-established for second-order (elliptic) differential equations, less is known for fourth-order problems to which (3.3) belongs when $\beta>0$. We are mainly interested in the fourth-order operator $\beta\Delta^{2}-\tau\Delta$ in a bounded smooth domain $D$ of $\mathbb{R}^{d}$, $d\geq 1$, with $\beta>0$ and $\tau\geq 0$. Let us point out, though, that the validity of the maximum principle heavily depends on the boundary conditions and the domain itself. The solution to the fourth-order equation with pinned boundary conditions (1.7), i.e.,

[TABLE]

satisfies $\phi(x)>0$, $x\in D$, provided that $f$ is a smooth and nonnegative function in $\bar{D}$ with $f\not\equiv 0$. This is readily seen by writing (3.10) as a system of two second-order equations, each of them involving a second-order linear elliptic operator with homogeneous Dirichlet boundary conditions on $\partial D$, so that the maximum principle applies.222Note that this is not true in non-smooth domains, see [45]. A similar result, however, fails to be true in general when the pinned boundary conditions (1.7) are replaced by the clamped boundary conditions (1.4) in which case the same transformation no longer works as one equation of the second-order system comes without and the other with overdetermined boundary conditions. Nevertheless, it was observed by Boggio [7] that, when $D$ is the unit ball $\mathbb{B}_{1}$ in $\mathbb{R}^{d}$, $d\geq 1$, and $\tau=0$, the associated Green function can be computed explicitly and is positive in $\mathbb{B}_{1}$. An obvious consequence of this positivity property is that solutions to

[TABLE]

with a non-negative right-hand side $f\not\equiv 0$ are positive in $\mathbb{B}_{1}$. Actually, Boggio’s celebrated result provides the impetus for several further studies, including the extension of the maximum principle to other domains [28, 51, 52, 67, 132] and other operators [50, 97, 123, 133]. See also the monograph [45] for a detailed discussion of positivity properties of higher-order elliptic operators.

Our particular interest for the maximum principle with regard to MEMS equations is the inclusion of second-order terms in the differential operator accounting for stretching. As the following result shows, this is possible in one dimension or in radial symmetry.

Proposition 3.5** (Maximum Principle [7, 50, 123, 97]).**

Let $\beta>0$ and $\tau\geq 0$.

1. (i)

If $f:\bar{\mathbb{B}}_{1}\rightarrow[0,\infty)$ is smooth with $f\not\equiv 0$, then the solution to

[TABLE]

satisfies $\phi(x)>0$, $x\in\mathbb{B}_{1}$. 2. (ii)

If $f:\bar{\mathbb{B}}_{1}\rightarrow[0,\infty)$ is smooth with $f\not\equiv 0$ and, in addition, radially symmetric when $d>1$, then the solution to

[TABLE]

satisfies $\phi(x)>0$, $x\in\mathbb{B}_{1}$. 3. (iii)

The eigenvalue problem

[TABLE]

has an eigenvalue $\mu_{1}>0$ with a corresponding radially symmetric (normalized) eigenfunction $\zeta_{1}>0$ in $\mathbb{B}_{1}$, and $\mu_{1}$ is the only eigenvalue with a positive radially symmetric eigenfunction. Moreover, any radially symmetric eigenfunction corresponding to the eigenvalue $\mu_{1}$ is a scalar multiple of $\zeta_{1}$.

As already mentioned the first statement of Proposition 3.5 is proved in [7]. The maximum principle stated in Proposition 3.5 (ii) is established in [123] for $d=1$ by means of an explicit computation of the associated Green function and later extended in [50] (by a different proof relying on the maximum principle for second-order equations) also to include third order terms. The higher-dimensional case in radial symmetry and the application to the principal eigenvalue mentioned in Proposition 3.5 (iii) were recently shown for a wider class of of fourth-order operators that can be written as a composition of two second-order operators [97].

It is important to point out that the maximum principle fails for

[TABLE]

when $p>0$ is too large. The maximum principle for evolution problems of the form (3.3) is thus not valid, not even if $\gamma=0$. Therefore, when $\beta>0$, this central tool can mainly be helpful for the study of stationary solutions to (3.3).

4. Stationary Problem

We now focus our attention on classical stationary solutions to MEMS equations which in particular have a minimum value strictly greater than $-1$. As pointed out in the introduction the pull-in instability of MEMS devices is supposed to be revealed by a sharp threshold $\lambda_{\varepsilon}^{stat}>0$ for the voltage parameter $\lambda$ above which no such stationary solution exists while below there is at least one. In Section 4.1 we collect the present knowledge on this matter for the free boundary problem (2.1)-(2.5). As we shall see, however, the picture regarding the pull-in voltage is not yet completely established. For the vanishing aspect ratio model (1.9), the study is understandably more developed as briefly described in Section 4.2.

4.1. Stationary Solutions to the Free Boundary Problem

If not stated otherwise, $\varepsilon>0$ is fixed in this section and we often omit the dependence on this parameter.

Stationary solutions $(u,\psi)$ to the free boundary problem (2.1)-(2.5) satisfy equivalently

[TABLE]

where $g_{\varepsilon}(u)$ with $\psi=\psi_{u}$ is given in (3.1), its properties being stated in Proposition 3.1. Of course, $u>-1$ in $[-1,1]$.

4.1.1. Non-Existence

As expected on physical grounds, stationary solutions do not exist if the applied voltage is too large. This is rather easily seen in the second-order case, that is, when $\beta=0$ in (4.1). Indeed, in this case, any stationary solution $u$ is clearly convex. Consequently, the function $(x,z)\mapsto 1+z-u(x)$ is a supersolution to (2.1)-(2.2), and the comparison principle implies

[TABLE]

hence

[TABLE]

according to (2.2). Consequently, $g_{\varepsilon}(u)\geq 1$ in $I$ and we infer from (4.1) (with $\beta=0$) that $\tau\partial_{x}^{2}u\geq\lambda$ in $I$. Multiplying this inequality by the nonnegative eigenfunction $\phi_{1}(x)=\cos(\pi x/2)$ of the Laplacian $-\partial_{x}^{2}$ in $I$ with homogeneous Dirichlet boundary conditions and integrating over $I$ readily give

[TABLE]

since $u>-1$. This shows that $\lambda$ cannot exceed $\tau\pi^{2}/4$ (this bound can be improved, though).

Unfortunately, in the fourth-order case $\beta>0$, stationary deflections with clamped boundaries are not convex and (4.2) thus no longer holds. Nevertheless, the non-existence result is still true, but requires a more sophisticated proof. It is based on a nonlinear variant of the eigenfunction method involving the positive eigenfunction $\zeta_{1}$ from Proposition 3.5 associated with the operator $\beta\partial_{x}^{4}-\tau\partial_{x}^{2}$ subject to clamped boundary conditions and normalized by $\|\zeta_{1}\|_{1}=1$.

Theorem 4.1** (Non-Existence [93, 95]).**

Let $(\beta,\tau)\in[0,\infty)^{2}\setminus\{(0,0)\}$. There exist $\varepsilon_{*}>0$ (with $\varepsilon_{*}=\infty$ if $\beta=0$) and ${\lambda_{ns}}:(0,\varepsilon_{*})\to(0,\infty)$ such that there is no stationary solution $(u,\psi_{u})$ to (2.1)-(2.5) for $\varepsilon\in(0,\varepsilon_{*})$ and $\lambda>{\lambda_{ns}}(\varepsilon)$.

Of course, in Theorem 4.1 it is expected that $\varepsilon_{*}=\infty$ also if $\beta>0$.

Proof.

Let $u$ be a stationary solution to (4.1) with $\beta>0$ and recall that it satisfies $-1<u<0$ in $I$ by Proposition 3.5. Hence the solution $\mathcal{U}\in H^{2}_{D}(I)$ to

[TABLE]

satisfies $-1/2\leq\mathcal{U}\leq 0$. We multiply (4.1) by $\zeta_{1}(1+\alpha\mathcal{U})$ with $\alpha=\min{\{1,\varepsilon^{2}\}}\in(0,1]$ and integrate over $I$. After integrating several times by parts and using the properties of $u$, $\mathcal{U}$, and $\zeta_{1}$ including the lower bound $1+\alpha\mathcal{U}\geq 1/2$, we obtain

[TABLE]

for some $K_{1}>0$ depending only on $\beta$ and $\tau$. Next, setting

[TABLE]

we infer from Young’s inequality that, for $\delta\in(0,1)$,

[TABLE]

while computations in the spirit of (3.7) and (3.8) and the non-positivity of $u$ give

[TABLE]

and

[TABLE]

Combining these inequalities leads to

[TABLE]

and the choice $\delta^{2}:=\lambda\varepsilon^{2}/(2\alpha\beta)$ cancels the last term of the right-hand side of the above inequality. Consequently,

[TABLE]

from which we readily deduce that $\lambda$ cannot exceed a threshold value depending on $\varepsilon$ provided $\varepsilon<\varepsilon_{*}:=1/\sqrt{\|\partial_{x}^{2}\zeta_{1}\|_{1}}$. In particular, for $\varepsilon\in(0,\min\{1,\varepsilon_{*}/2\})$ we conclude from the above bound and the choice of $\alpha$ that

[TABLE]

which does not depend on $\varepsilon$. ∎

4.1.2. Existence

To establish the existence of stationary solutions for small values of $\lambda$ we introduce the operator $A:=\beta\partial_{x}^{4}-\tau\partial_{x}^{2}\in\mathcal{L}(D(A),L_{p}(I))$ with $D(A):=H_{D}^{4}(I)$ and $p=2$ if $\beta>0$ and $D(A):=W_{p,D}^{2}(I)$ with $p>2$ if $\beta=0$ and $\tau>0$. Note that $A$ is invertible. Given $q>2$, Proposition 3.1 implies that solutions to (4.1) are the zeros of the smooth function

[TABLE]

The Implicit Function Theorem yields a smooth curve $(\lambda,U_{\lambda})$ in $\mathbb{R}\times D(A)$ of zeros starting at $(\lambda,v)=(0,0)$ which can be maximally extended as long as the principal eigenvalue of the associated linearized operator $A+\lambda\partial_{u}g_{\varepsilon}(U_{\lambda})$ in $L_{p}(I)$ is positive.

Theorem 4.2** (Existence [95, 30]).**

Let $(\beta,\tau)\in\mathbb{R}_{+}^{2}\setminus\{(0,0)\}$. There is $\lambda_{s}(\varepsilon)>0$ such that for each $\lambda\in(0,\lambda_{s}(\varepsilon))$ there exists a stationary solution $(U_{\lambda},\Psi_{U_{\lambda}})$ to (2.1)-(2.5) satisfying $-1<U_{\lambda}<0$ in $I$ and the associated linearized operator $A+\lambda\partial_{u}g_{\varepsilon}(U_{\lambda})$ in $L_{p}(I)$ has a positive principal eigenvalue. In addition, $U_{\lambda}=U_{\lambda}(x)$ and $\Psi_{U_{\lambda}}=\Psi_{U_{\lambda}}(x,z)$ are even with respect to $x\in I$.

The sign property of $U_{\lambda}$ is a consequence of Proposition 3.5. Actually, one can show that these stationary solutions are locally asymptotically stable for the evolution [95, 30]. Obviously, the two critical values found in Theorem 4.1 and Theorem 4.2 satisfy $\lambda_{ns}(\varepsilon)\geq\lambda_{\varepsilon}^{stat}\geq\lambda_{s}(\varepsilon)$, but otherwise no further information is yet available. In particular, the occurrence of a sharp threshold in form of a pull-in instability (i.e. $\lambda_{ns}(\varepsilon)=\lambda_{s}(\varepsilon)$) is an open problem. Even the simpler problem of identifying the behavior of $(U_{\lambda},\Psi_{U_{\lambda}})$ as $\lambda\to\lambda_{s}(\varepsilon)$ seems currently out of reach.

4.1.3. Multiplicity

Theorem 4.2 leaves open the question whether more than one stationary solution exists for small values of $\lambda$ (as suggested by what is known for the vanishing aspect ratio equation, see Section 4.2 below). Based on the variational structure of (4.1) pointed out in Proposition 3.2 it is possible to show that multiple stationary solutions indeed exist for (some) small values of $\lambda$. Recall that the total energy $\mathcal{E}$ is non-coercive and a plain minimization thereof is thus not effective. However, when $\beta>0$ one can find additional critical points of $\mathcal{E}$ by considering a constrained minimization problem, where the idea is to minimize the mechanical energy $\mathcal{E}_{m}$ on a set of deflections with constant electrostatic energy $\mathcal{E}_{e}$.

Theorem 4.3** (Multiplicity [91]).**

Let $\beta>0$ and $\tau\geq 0$. Given $\rho>2$, there is $(\lambda_{\rho},u_{\rho},\psi_{\rho})$ with $\lambda_{\rho}>0$, $u_{\rho}\in H_{D}^{4}(I)$, and $\psi_{u_{\rho}}\in H^{2}(\Omega(u_{\rho}))$ such that $(u_{\rho},\psi_{u_{\rho}})$ is a stationary solution to (2.1)-(2.4) with $\lambda=\lambda_{\rho}$. Both $u_{\rho}=u_{\rho}(x)$ and $\psi_{u_{\rho}}=\psi_{u_{\rho}}(x,z)$ are even with respect to $x\in I$ and $-1<u_{\rho}<0$ in $I$. Moreover, $\lambda_{\rho}\rightarrow 0$ as $\rho\rightarrow\infty$ and $u_{\rho}\not=U_{\lambda_{\rho}}$ for all $\rho>2$ sufficiently large, the functions $U_{\lambda}$ being defined in Theorem 4.2.

Proof.

For a given $\rho>2$ define the (non-empty) set

[TABLE]

and

[TABLE]

Based on Proposition 3.3 it can be shown that the function $\mu$ is non-decreasing on $(2,\infty)$ with

[TABLE]

Hence, taking a minimizing sequence $(u_{k})_{k\geq 1}$ of $\mathcal{E}_{m}$ in $\mathcal{A}_{\rho}$ satisfying

[TABLE]

one can combine (4.3) with the energy inequality (3.7) to show that there is a constant $\delta(\rho)\in(0,1)$ such that

[TABLE]

Since the sequence \big{(}\|\partial_{x}^{2}u_{k}\|_{2}\big{)}_{k} is bounded due to (4.3), one may extract a subsequence of $(u_{k})_{k}$ which converges weakly to some $u_{\rho}$ in $H_{D}^{2}(I)$ still satisfying (4.4). Proposition 3.2 entails $\mathcal{E}_{e}(u_{\rho})=\rho$ so that $u_{\rho}$ is indeed a minimizer on $\mathcal{A}_{\rho}$ of $\mathcal{E}_{m}$ due to the latter’s weak lower semi-continuity. Consequently, there is a Lagrange multiplier $\lambda_{\rho}\in\mathbb{R}$ such that

[TABLE]

which, combined with classical elliptic regularity theory and Proposition 3.2, implies that $u_{\rho}\in H_{D}^{4}(I)$ is a classical solution to (4.1). Finally, owing to (3.8) and the properties of $u_{\rho}$ one may derive a bound for the Lagrange multiplier $\lambda_{\rho}$ of the form

[TABLE]

which shows that $\lambda_{\rho}\rightarrow 0$ as $\mathcal{E}_{e}(u_{\rho})=\rho\rightarrow\infty$ . In particular, $u_{\rho}\not=U_{\lambda_{\rho}}$ for all $\rho>2$ sufficiently large, since $\mathcal{E}_{e}(U_{\lambda})\rightarrow 2$ as $\lambda\rightarrow 0$. ∎

Theorem 4.3 now yields multiplicity of stationary solutions for certain small values of $\lambda$. Indeed, there is at least a sequence $\lambda_{j}\rightarrow 0$ of voltage values for which there are two different solutions $(u_{j},\psi_{u_{j}})$ (i.e. $\rho=j$ in Theorem 4.3) and $(U_{\lambda_{j}},\Psi_{U_{\lambda_{j}}})$ (i.e. $\lambda=\lambda_{j}$ in Theorem 4.2). Note that, by taking a different sequence $\rho_{j}\rightarrow\infty$ with $\rho_{j}\not=j$, one obtains different solutions $(u_{\rho_{j}},\psi_{u_{\rho_{j}}})$ but with possibly equal voltage values. We conjecture that the solutions constructed in Theorem 4.3 actually lie on a smooth curve. This is supported by the numerical simulations performed in [43] as well as by the results known for the simplified vanishing aspect ratio model

[TABLE]

supplemented with the boundary conditions (2.4), see Section 4.2 below. It is worth pointing out that the transformation for $\psi$ to a fixed domain introduced in Section 3.1 is also used to perform the numerical simulations in [43].

A multiplicity result as stated in Theorem 4.3 is not known up to now for the second-order case $\beta=0$ (but $\tau>0$), the main reason being the weaker regularity that comes along which appears to be insufficient to provide a control of the electrostatic energy $\mathcal{E}_{e}$ by the mechanical energy $\mathcal{E}_{m}$.

4.1.4. Vanishing Aspect Ratio Limit

Of course, the value $\lambda_{s}(\varepsilon)$ as well as the stationary solutions

[TABLE]

derived in Theorem 4.2 and Theorem 4.3, respectively, depend on the aspect ratio $\varepsilon$. It is then of primary interest to investigate the behavior of these solutions in the limit $\varepsilon\to 0$ and show that this limiting process produces stationary solutions to the stationary vanishing aspect ratio model (4.5). In particular, this puts the derivation of this model on firm ground which has been performed formally up to now. Roughly speaking, given a family $(u_{\varepsilon},\psi_{u_{\varepsilon}})_{\varepsilon}$ of stationary solutions to (2.1)-(2.5) with $\lambda=\lambda_{\varepsilon}$, the main steps of the proof are the following:

• (i)

Derive positive upper and lower bounds on $(\lambda_{\varepsilon})_{\varepsilon}$, which do not depend on $\varepsilon$ small enough.

• (ii)

Derive estimates on $(u_{\varepsilon})_{\varepsilon}$ in a suitable Sobolev space, e.g. in $W_{q}^{2}(I)$ for some $q\geq 2$, together with a positive lower bound on $(1+u_{\varepsilon})_{\varepsilon}$, again not depending on $\varepsilon$ small enough.

• (iii)

Estimate the difference between $\psi_{u_{\varepsilon}}$ and $(x,z)\mapsto(1+z)/(1+u_{\varepsilon}(x))$ in a Sobolev space by positive powers of $\varepsilon$, which, in turn, provides an estimate of $g_{\varepsilon}(u_{\varepsilon})-(1+u_{\varepsilon})^{-2}$ and shows that it converges to zero as $\varepsilon\to 0$.

Once these three properties are established, a compactness argument is used to obtain the following convergence result:

Theorem 4.4** (Vanishing Aspect Ratio Limit [93, 92]).**

There are a sequence $\varepsilon_{k}\to 0$, $\lambda_{0}>0$, and a (smooth) stationary solution $u_{0}$ to the stationary vanishing aspect ratio equation (4.5) with $\lambda=\lambda_{0}$ such that

[TABLE]

and

[TABLE]

Let us emphasize that Theorem 4.4 above is valid either when $\beta>0$ [92] or when $\beta=0$ and $\tau>0$ [93]. However, the starting point of the analysis and actually the derivation of properties (i)–(iii) listed above differ in both cases. In the second-order case $\beta=0$ stationary solutions for small voltage values can also be constructed by means of Schauder’s fixed point theorem [93] and this approach guarantees in particular a lower bound $\lambda_{min}>0$ for $\lambda_{s}(\varepsilon)$ independent of $\varepsilon\in(0,1)$. This thus gives a family of stationary solutions $(u_{\varepsilon},\psi_{u_{\varepsilon}})_{\varepsilon\in(0,1]}$ for each $\lambda\in(0,\lambda_{min}]$. Moreover, it yields a priori bounds on the stationary solutions of the form

[TABLE]

for each $\lambda\in(0,\lambda_{min}]$, where $\kappa_{0}\in(0,1)$ is independent of $\varepsilon\in(0,1)$. In that case, the family $(\lambda_{\varepsilon})_{\varepsilon}$ in Theorem 4.4 is given by $\lambda_{\varepsilon}=\lambda\in(0,\lambda_{min}]$ and the properties (i)–(ii) are satisfied. In the fourth-order case $\beta>0$, the outcome of Theorem 4.3 is the starting point. Given $\rho>2$ and $\varepsilon>0$, we set $(\lambda_{\varepsilon},u_{\varepsilon}):=(\lambda_{\rho}^{\varepsilon},u_{\rho}^{\varepsilon})$, the latter being constructed in Theorem 4.3 and depending on $\varepsilon$, and use the variational construction of $u_{\varepsilon}$ to derive the uniform estimates (i)–(ii). Property (iii) follows from (2.1) after multiplication by $\psi_{u_{\varepsilon}}$ and $\partial_{z}\psi_{u_{\varepsilon}}$. Note that the regularity of $\psi_{u_{\varepsilon}}$ is no longer the same in the $x$- and $z$-directions in the limit $\varepsilon\to 0$ since (2.1) becomes degenerate elliptic. Hence, a cornerstone of the proof is to obtain estimates for the trace of $\partial_{z}\psi_{u_{\varepsilon}}$ on the upper boundary of $\Omega(u_{\varepsilon})$ [93, 92].

4.2. Stationary Solutions to the Vanishing Aspect Ratio Model

The study of stationary solutions to the vanishing aspect ratio model (1.9) has a long history, starting with the seminal work [76] on the number of solutions to a broad class of semilinear Dirichlet problems in radial symmetry and then continued to cover second-order MEMS equations

[TABLE]

with open domain $D\subset\mathbb{R}^{d}$ (here $d$ is no longer restricted to $1,2$) and various variants hereof, including those listed in Section 2.4. As expected from physics, there is a critical value $\lambda_{0}^{stat}>0$ – depending on the domain’s geometry – such that no stationary solution exists if $\lambda>\lambda_{0}^{stat}$ and at least one stationary solution exists for $\lambda\in(0,\lambda_{0}^{stat})$. More precisely, owing to the monotonicity and convexity of the function $z\mapsto 1/(1+z)^{2}$, an iterative monotone scheme can be set up to construct a branch of stationary solutions $\{U_{\lambda}\,;\,\lambda\in[0,\lambda_{0}^{stat})\}$ to (4.8) with $U_{0}\equiv 0$ and $U_{\lambda^{\prime}}<U_{\lambda}$ in $I$ for all $0\leq\lambda<\lambda^{\prime}<\lambda_{0}^{stat}$. A consequence of this approach is that, given $\lambda\in(0,\lambda_{0}^{stat})$, each stationary solution $U_{\lambda}$ is also a maximal and stable solution in the following sense: Any other solution to (4.8) with the same value of $\lambda$ lies below $U_{\lambda}$, and the principal eigenvalue of the associated linearized operator $-\Delta-2\lambda/(1+U_{\lambda})^{3}$ in $L_{2}(D)$ is positive. The behavior of $U_{\lambda}$ as $\lambda\to\lambda_{0}^{stat}$ as well as a further and more precise characterization of the bifurcation diagram depend strongly on the dimension $d$. Concerning the former, owing to the monotonicity of the family $\{U_{\lambda}\,;\,\lambda\in[0,\lambda_{0}^{stat})\}$ with respect to $\lambda$ and the lower bound $U_{\lambda}>-1$, the function

[TABLE]

is always well-defined but might reach the value $-1$ at some points. However, this cannot occur for $d\in[1,7]$, and $U_{\lambda_{0}^{stat}}$ is then a classical solution to (4.8) with $\lambda=\lambda_{0}^{stat}$, see [47]. If $d\geq 8$, then touchdown points may exist in the sense that $\min U_{\lambda_{0}^{stat}}=-1$: this is in particular the case when $D$ is the unit ball $\mathbb{B}_{1}$ of $\mathbb{R}^{d}$, see [47]. As for the bifurcation diagram, there are exactly two stationary solutions for each $\lambda\in(0,\lambda_{0}^{stat})$ if $d=1$, see [76]. For dimensions $d\in[2,7]$ there is a curve $s\mapsto(\Lambda(s),U(s))$ of stationary solutions starting for $s=0$ at $(0,0)$ with infinitely many bifurcation or turning points so that there exists a unique (under additional geometric restrictions on $D$) stationary solution for sufficiently small voltage values $\lambda$, exactly two stationary solutions for $\lambda$ in a small left neighborhood of $\lambda_{0}^{stat}$, one extremal stationary solution at $\lambda=\lambda_{0}^{stat}$, and multiple stationary solutions for some values of $\lambda$ strictly between zero and $\lambda_{0}^{stat}$, see [33, 34, 35]. Additional properties of stationary solutions were investigated in dependence on the dimension of the ambient space and the geometry of the device. Proofs (and numerical evidence) for these results as well as regarding e.g. to estimates on the pull-in voltage or the pull-in distance (that is, the minimum of $1+U_{\lambda_{0}^{stat}}$) are to be found in [6, 129, 131, 127, 61, 38, 47, 33, 34, 4, 22, 65], see also [35] for an extended review. It is worth mentioning that the comparison principle is available in this case and a key tool for the analysis in most of the papers cited above.

Less attention has so far been dedicated to the fourth-order problem

[TABLE]

with $\beta>0$. The only result concerning (4.9)-(4.10) in a general domain $D$ – hence not relying on the maximum principle – seems to be [107], where stationary solutions have been constructed for small values of $\lambda$ by Banach’s fixed point theorem.

When restricting to the unit ball $\mathbb{B}_{1}$ of $\mathbb{R}^{d}$, $d\geq 1$, so that one may take advantage of the maximum principle in radial symmetry as stated in Proposition 3.5, similar arguments as in the second-order case provides a threshold value $\lambda_{0}^{stat}\in(0,\infty)$ such that, for $\lambda>\lambda_{0}^{stat}$, there is no stationary solution if $\tau=0$ and no radially symmetric stationary solution if $\tau>0$. Furthermore, there is a monotone branch of stable and maximal radially symmetric stationary solutions $\{U_{\lambda}\,;\,\lambda\in[0,\lambda_{0}^{stat})\}$ with $U_{0}=0$, see [13] for $\tau=0$ and [97, 94] for $\tau>0$. The function $U_{\lambda_{0}^{stat}}:=\inf_{\lambda\in(0,\lambda_{0}^{stat})}U_{\lambda}$ is again well-defined and its properties are investigated in [21] when $\tau=0$. In particular, it is shown that the minimum of $U_{\lambda_{0}^{stat}}$ in $\mathbb{B}_{1}$ is above $-1$ for $d\leq 8$ so that $U_{\lambda_{0}^{stat}}$ is a smooth solution to (4.9)-(4.10) in that case. In contrast, the minimum of $U_{\lambda_{0}^{stat}}$ in $\mathbb{B}_{1}$ is equal to $-1$ if $d\geq 9$. Additional information on the structure of the solution set is provided in [25] for $d\geq 3$.

In dimensions $d=1,2$ more precise information on the bifurcation diagram is available, even when including $\tau>0$:

Theorem 4.5** (Bifurcation Diagram [94]).**

Let $\beta>0$, $\tau\geq 0$, and consider $D=\mathbb{B}_{1}$ with $d\in\{1,2\}$. There exists $\lambda_{0}^{stat}>0$ such that there is no radially symmetric solution to (4.9)-(4.10) for $\lambda>\lambda_{0}^{stat}$, at least one radially symmetric solution for $\lambda=\lambda_{0}^{stat}$, and at least two radially symmetric solutions for $\lambda\in(0,\lambda_{0}^{stat})$, a stable and an unstable one. These solutions together with the corresponding values of $\lambda$ lie on an analytic curve and, as $\lambda$ approaches zero, the unstable solution converges to an explicitly computable function $\omega$ with minimum equal to $-1$.

Theorem 4.5 for $\tau\geq 0$ gives a positive answer to the question of existence stated in [35] of at least two solutions for each $\lambda\in(0,\lambda_{0}^{stat})$ in dimension $d=2$. When $\tau=0$ the existence of at least two solutions for each $\lambda\in(0,\lambda_{0}^{stat})$ is also shown in [62] in the two-dimensional case $d=2$. A second solution different from the stable one is constructed therein by the Mountain Pass Principle. Its behavior as $\lambda\rightarrow 0$ is studied and shown to coincide with the one described above. Let us also point out that the asymptotic shape $\omega$ was found in [114] but by formal asymptotic expansions.

Theorem 4.5 relies on the construction of an analytic curve $s\mapsto(\Lambda(s),U(s))$ in $\mathbb{R}\times C^{4}(\mathbb{B}_{1})$ defined for $s\in[0,\infty)$ such that $U(s)$ is for each $s\in[0,\infty)$ a radially symmetric solution to (4.9)-(4.10) with $\lambda=\Lambda(s)$. The curve emanates from $(\Lambda(0),U(0))=(0,0)$ with $(\Lambda(s),U(s))\rightarrow(0,\omega)$ as $s\rightarrow\infty$, where the end point $\omega$ is given as a solution of a boundary value problem in $\mathbb{B}_{1}\setminus\{0\}$ which can be computed explicitly. A plot of $\omega$ is shown in Figure 4, its qualitative behavior being the same for $d=1$ and $d=2$.

In addition, there is $s_{*}>0$ such that $\Lambda$ is an increasing function from $[0,s_{*}]$ onto $[0,\lambda_{0}^{stat}]$ and $\Lambda$ is decreasing in a right-neighborhood of $s_{*}$ with $\Lambda(s)\rightarrow 0$ as $s\rightarrow\infty$. Consequently, given $\lambda\in(0,\lambda_{0}^{stat})$ there are indeed two values $0<s_{1}<s_{*}<s_{2}$ with $\Lambda(s_{j})=\lambda$ for $j=1,2$ and $U(s_{2})\leq U(s_{1})$ in $\mathbb{B}_{1}$ with $U(s_{2})\not=U(s_{1})$.

The proof of Theorem 4.5 is performed in several steps. The implicit function theorem provides a smooth branch $\mathcal{A}_{0}$ of radially symmetric solutions $(\lambda,u)$ to (4.9)-(4.10) emanating from $(0,0)$ which can be extended to a global curve $\mathcal{A}$ by means of bifurcation theory for real analytic functions [12], a tool that was also used in [35, Sect. 6.2] for the second-order case $\beta=0$. It can then be shown that the branch $\{(\lambda,U_{\lambda})\,;\,\lambda\in[0,\lambda_{0}^{stat})\}$ of stable and maximal solutions to (4.9)-(4.10) constructed by the iterative monotone scheme described above actually coincides with $\mathcal{A}_{0}$. A priori estimates and the fact that $U_{\lambda}$ is decreasing in $\lambda$ entail that a radially symmetric stationary solution $U_{\lambda_{0}^{stat}}$ exists also at $\lambda=\lambda_{0}^{stat}$, which guarantees, on the one hand, that $\mathcal{A}_{0}\not=\mathcal{A}$ and, on the other hand, that the bifurcation result of [24] applies to extend the branch $\mathcal{A}_{0}$ beyond $(\lambda_{0}^{stat},U_{\lambda_{0}^{stat}})$. The final and most intricate step is to show that $\mathcal{A}$ connects $(0,0)$ to the end point $(0,\omega)$ and to identify the latter. As a consequence, the continuous curve $\mathcal{A}$ passes through $(\lambda_{0}^{stat},U_{\lambda_{0}^{stat}})$.

Remarks 4.6**.**

1. (1)

When $\tau>0$ the maximum principle is only valid for radially symmetric functions by Proposition 3.5 and we can only guarantee the non-existence of radially symmetric stationary solutions to (4.9)-(4.10) for $\lambda>\lambda_{0}^{stat}$. It is possible (though unlikely) that non-radially symmetric stationary solution exist above this threshold value. 2. (2)

The qualitative bifurcation diagram is expected to be completely different in higher space dimensions $d\geq 3$ as pointed out in **[25, 35]**.

For numerical results on (4.9)-(4.10) and investigations based on formal asymptotic methods we refer to [114, 113, 86]. In particular, by means of a shooting argument one can numerically verify that (4.9)-(4.10) admits no other solutions than found in Theorem 4.5 [113].

We end this section with a few words about (4.9) when it is supplemented with pinned boundary conditions (1.7), that is,

[TABLE]

with $\beta>0$ and $\tau\geq 0$. Interestingly, the maximum principle holds in arbitrary smooth domains [35, 45] and allows one at least in part to show similar results for the fourth-order problem as outlined above for the second-order case. In particular, there is again $\lambda_{0}^{stat}>0$ such that no solution to (4.11) exists for $\lambda>\lambda_{0}^{stat}$ while there is a branch of stable maximal solutions $\{U_{\lambda}\,;\,\lambda\in(0,\lambda_{0}^{stat})\}$ obtained from a monotone iterative scheme such that $U_{\lambda}$ is decreasing in $\lambda$, see [107]. If $d\leq 3$, then the minimum of $U_{\lambda_{0}^{stat}}:=\inf_{\lambda\in(0,\lambda_{0}^{stat})}U_{\lambda}$ is greater than $-1$ and $U_{\lambda_{0}^{stat}}$ is the unique classical solution to (4.9) subject to (1.7), see [66]. A similar result is available in [20, 23] when $\tau=0$ and $d\leq 6$. In the particular case of $D$ being a ball, the extremal solution is still regular in the sense that $\inf_{D}U_{\lambda_{0}^{stat}}>-1$ for dimension $d\leq 8$ [21, 118] while a touchdown in the sense of $\inf_{D}U_{\lambda_{0}^{stat}}=-1$ occurs for $d\geq 9$ (provided $\tau/\beta$ is small) [118]. If $D$ is a two-dimensional convex domain, then there is a second stationary solution for $\lambda\in(0,\lambda_{0}^{stat})$ constructed by the Mountain Pass Principle that approaches $-1$ as $\lambda\rightarrow 0$ [66, 62]. If $D$ is the three-dimensional unit ball, then (4.11) with $\tau=0$ has infinitely many fold points for the bifurcation branch corresponding to radially symmetric solutions [64]. For more precise and further results for (4.11) we refer to the abovementioned references [63, 66, 35, 107, 20, 64, 118, 21, 23, 62] as well as to [122, 14].

5. Evolution Problem

We next turn to the dynamical free boundary MEMS equations (2.1)-(2.5). Let us recall that it is expected on physical grounds that there is a (dynamical) pull-in voltage value determining the dynamics of the device in the following sense: Solutions corresponding to values of $\lambda$ below this threshold are global in time, while for values above the threshold solutions cease to exist globally as the elastic plate pulls in at some finite time $T_{m}<\infty$, i.e.,

[TABLE]

Rigorous proofs for such a qualitative behavior with a sharp threshold are obviously non-trivial, particularly for the free boundary problem (2.1)-(2.5). Nevertheless, some of the described behavior is established. As we shall see in the next sections, the fourth-order case $\beta>0$ and the second-order case $\beta=0$ are somewhat different with respect to what is proven so far for large voltage values. In the former case, non-global solutions necessarily exhibit the touchdown behavior described in (5.1) but it is not known (though expected) that such solutions indeed exist for large $\lambda$. In the latter case, solutions cannot exist globally in time for large $\lambda$ but may exhibit also an alternative singularity in form of a norm blow up. Roughly speaking, the reason for this difference is that the fourth-order case is technically more involved, but comes along with a higher regularity.

We first consider in Section 5.1 the local well-posedness and then in Sections 5.2-5.4 questions related to global existence. The dependence of the solutions on some parameters is then considered: in Section 5.5 the damping dominated limit $\gamma\rightarrow 0$ is studied and in Section 5.6 the vanishing aspect ratio limit $\varepsilon\rightarrow 0$.

5.1. Local Existence

We recall that one can formulate (2.1)-(2.5) as a single nonlocal problem of the form (3.3) only involving $u$, the corresponding electrostatic potential $\psi=\psi_{u}$ being determined with the help of Proposition 3.1. Specifically, recalling that $g_{\varepsilon}(u)$ is the square of the trace of the (rescaled) gradient of the solution $\psi_{u}$ to (2.1)-(2.2) on the upper boundary of $\Omega(u)$, see (3.1), the function $u$ solves

[TABLE]

with clamped boundary conditions

[TABLE]

and zero initial conditions

[TABLE]

the latter only assumed for a simpler exposition. To discuss the solvability of (5.2)-(5.4) we fix $2\alpha\in(0,1/2)$ and introduce the operator

[TABLE]

with $\ell=2$ if $\beta>0$ and $\ell=1$ if $\beta=0$ and $\tau>0$.

Let us first consider the case $\beta>0$ and $\gamma>0$. We reformulate (5.2)-(5.4) as a hyperbolic nonlocal semilinear Cauchy problem for ${\bf u}=(u,\gamma\partial_{t}u)^{T}$ of the form

[TABLE]

with ${\bf u}^{\prime}$ indicating the time derivative of ${\bf u}$. Hereby, the matrix operator

[TABLE]

with domain $D(\mathbb{A}_{\alpha}):=H_{D}^{4+2\alpha}(I)\times H_{D}^{2+2\alpha}(I)$ generates a strongly continuous group $e^{-t\mathbb{A}_{\alpha}}$, $t\in\mathbb{R}$, on the Hilbert space $\mathbb{H}_{\alpha}:=H_{D}^{2+2\alpha}(I)\times H_{D}^{2\alpha}(I)$ with an exponential decay, the latter being due to the damping term $\partial_{t}u$ in (5.2) (e.g. see the cosine function theory in [1, Sections 5.5 & 5.6]). Since $H_{D}^{2+2\alpha}(I)$ is continuously embedded in $W_{q,D}^{2}(I)$ for some $q>2$ and $\alpha\in(0,1/4)$, Proposition 3.1 entails that the function

[TABLE]

with

[TABLE]

is bounded and uniformly Lipschitz continuous when $\kappa\in(0,1)$. A classical application of the contraction mapping principle then yields a unique strong solution ${\bf u}\in C([0,T],\mathbb{H}_{\alpha})$ to (5.5) with ${\bf u}^{\prime}\in L_{1}(0,T;\mathbb{H}_{\alpha})$ on a (possibly small) time interval $[0,T]$, which can be continued as long as $\|{\bf u}(t)\|_{\mathbb{H}_{\alpha}}$ stays bounded and $u(t)>-1$. This in particular yields a criterion for global existence.

Let us point out that, unfortunately, this approach does not work in the second-order case $\beta=0$ with $\gamma>0$ in (5.2)-(5.4) since the right-hand side $g_{\varepsilon}$ is – roughly speaking – of order slightly above one (see Proposition 3.1) and thus exceeds half of the operator’s order.

The just described difficulty does not occur if $\gamma=0$ due to parabolic smoothing effects. Indeed, in this case we may consider $\beta>0$ and $\beta=0$ (with $\tau>0$ then) and handle equation (5.2)-(5.4) as a parabolic nonlocal Cauchy problem in $H_{D}^{2\alpha}(I)$ of the form

[TABLE]

possessing a classical solution in C\big{(}[0,T_{m}),H_{D}^{2\ell+2\alpha}(I)\big{)}\cap C^{1}\big{(}[0,T_{m}),H_{D}^{2\alpha}(I)\big{)}.

In summary, we have the following result regarding the local well-posedness of (5.2)-(5.4):

Theorem 5.1** (Local Existence [95, 30]).**

Assume that either $\beta>0$, $\gamma\geq 0$, $\tau\geq 0$ or $\beta=\gamma=0$, $\tau>0$ and consider $\lambda>0$.

• (i)

There is a unique (smooth) solution $(u,\psi_{u})$ to (5.2)-(5.4) on a maximal interval of existence $[0,T_{m})$ such that

[TABLE]

• (ii)

If the solution does not exist globally in time, that is, if $T_{m}<\infty$, then a touchdown singularity (5.1) occurs or some Sobolev norm of \big{(}u(t),\gamma\partial_{t}u(t)\big{)} blows up as $t\rightarrow T_{m}$.

As pointed out above, the regularity of the solution depends on the particular case one considers, that is, on whether or not the parameters $\beta$ and $\gamma$ are zero. For further use in Sections 5.5 and 5.6 let us emphasize that the solution $(u,\psi_{u})$ as well as the maximal existence time $T_{m}$ depend on both values $\varepsilon>0$ and $\gamma\geq 0$.

5.2. Refined Criterion for Global Existence

Clearly, the outcome of Theorem 5.1 (ii) is not yet satisfactory as it does not allow one to figure out whether the occurring finite time singularity is a touchdown or a norm blowup of $u(t)$ as $t\rightarrow T_{m}<\infty$, the latter being though a non-observed phenomenon in real world applications.

We shall thus address next the preclusion of a norm blowup for non-global solutions. This requires a priori estimates on solutions for which a promising source may be the gradient flow structure of (5.2)-(5.4) pointed out in Section 3.2. However, when $\beta=0$ the best estimates one can hope for coming from the functional $\mathcal{E}$ in (3.6) is an $H^{1}(I)$-estimate on $u(t)$ (by definition of the mechanical energy $\mathcal{E}_{m}$) which seems to be too weak, unfortunately, as can be seen from (3.7). In contrast, the additional regularity stemming from the fourth-order bending term $\beta\partial_{x}^{4}$ is more promising as it leads to sufficient $H^{2}(I)$-estimates on $u(t)$. The drawback when including fourth-order bending is, though, that no (dynamical) comparison principle is available as noticed in Section 3.2 and so even a preliminary upper bound on $u$ is not obvious (again in contrast to the second-order case $\beta=\gamma=0$ where the comparison principle applied to (5.2) immediately implies the bound $u(t)\leq 0$ since $g_{\varepsilon}(u)>0$). A further difficulty encountered in deriving a priori estimates is the non-coercivity of the corresponding functional $\mathcal{E}$ from (3.6). We now give more details on how to overcome these issues.

The strategy to prove that necessarily the touchdown singularity (5.1) occurs in the fourth-order case $\beta>0$ when $T_{m}<\infty$ is the following: Supposing that there are $T>0$ and $\kappa\in(0,1)$ such that

[TABLE]

Theorem 5.1 (ii) implies $T_{m}\geq T$ provided one can control a suitable Sobolev norm of $(u(t),\gamma\partial_{t}u(t))$, see the discussion in Section 5.1. To find such a bound one first observes that – even though there is no a priori upper bound on $u$ as pointed out above – a weighted $L_{1}$-bound of the form

[TABLE]

for sufficiently small $\gamma\geq 0$ can be derived, where $\zeta_{1}$ is the positive eigenfunction of the operator $\beta\partial_{x}^{4}-\tau\partial_{x}^{2}$ subject to clamped boundary conditions provided by Proposition 3.5. This estimate is obtained by testing equation (5.2) against $\zeta_{1}$ and combining the outcome of the computation with the lower bound on $u$. Using the estimate (3.7) on the electrostatic energy $\mathcal{E}_{e}(u(t))$ together with an interpolation inequality and a variant of Poincaré’s inequality, one deduces

[TABLE]

for $\delta\in(0,1)$ and $t\in[0,T]\cap[0,T_{m})$. Recognizing part of the mechanical energy $\mathcal{E}_{m}(u(t))$ on the right-hand side of this inequality, the weighted $L_{1}$-estimate on $u$ from above and Corollary 3.4 yield (for sufficiently small $\delta>0$)

[TABLE]

hence

[TABLE]

Applying Proposition 3.1 (ii) then gives a suitable bound on $g_{\varepsilon}(u(t))$ and thus on the right-hand side of the Cauchy problem (5.5), what readily yields the desired a priori bound on $(u,\gamma\partial_{t}u)$. Consequently, if $T_{m}<\infty$, then (5.1) occurs.

Summarizing, we have the following result regarding global existence for the fourth-order problem.

Theorem 5.2** (Criterion for Global Existence when $\beta>0$ [95]).**

Let $\beta>0$ and $\tau\geq 0$. There is $\gamma_{0}>0$ such that if $\gamma\in[0,\gamma_{0})$, then $T_{m}<\infty$ implies the touchdown behavior (5.1).

5.3. Global Existence

An interesting byproduct of the fixed point argument leading to local existence is that it simultaneously yields global existence for small voltage values $\lambda$ in any of the cases considered in Theorem 5.1.

Theorem 5.3** (Global Existence [30, 95]).**

Assume that either $\beta>0$, $\gamma\geq 0$, $\tau\geq 0$ or $\beta=\gamma=0$, $\tau>0$. There is $\lambda_{ge}(\varepsilon,\gamma)>0$ such that the solution to (5.2)-(5.4) with maximal interval of existence $[0,T_{m})$ provided by Theorem 5.1 is global for any $\lambda\in(0,\lambda_{ge}(\varepsilon,\gamma))$, that is, $T_{m}=\infty$. In addition, $u$ is bounded in $H^{2+2\alpha}(I)$ uniformly with respect to time, where $\alpha\in(0,1/4)$, and

[TABLE]

Finally, any cluster point in $H^{2+\alpha}(I)$ of $\{u(t)\ ;\ t\geq 0\}$ as $t\to\infty$ is a stationary solution to (5.2)-(5.4).

The damping term $\partial_{t}u$ in (5.2) is crucial for the validity of Theorem 5.3. Indeed, the proof of the latter mainly relies on the exponential decay of the corresponding underlying semigroup which entails that the contraction mapping principle used to establish Theorem 5.1 can be applied on any (finite) time interval when $\lambda$ is small (see [95, 30] for details). Note that no touchdown occurs in case of Theorem 5.3, not even in infinite time.

To show the last statement of Theorem 5.3 observe that the bounds on $u$ stated in Theorem 5.3 and (3.7) entail that $\mathcal{E}(u)$ is bounded uniformly with respect to time and that $\{u(t)\ ;\ t\geq 0\}$ is relatively compact in $H^{2+\alpha}(I)$. Combining the bound on $\mathcal{E}(u)$ with (3.9) shows that

[TABLE]

The compactness of the trajectory entails that there is a sequence $t_{n}\rightarrow\infty$ and $w\in H^{2+\alpha}(I)$ such that $u(t_{n})\longrightarrow w$ in $H^{2+\alpha}(I)$. Classical arguments of dynamical systems then guarantee that $(w,\psi_{w})$ is a stationary solution to (5.2)-(5.4). A further consequence of the last statement in Theorem 5.3 is that $\lambda_{ge}(\varepsilon,\gamma)\leq\lambda_{ns}(\varepsilon)$, the latter being defined in Theorem 4.1.

5.4. Finite Time Singularity

Showing the physically expected behavior that solutions to (5.2)-(5.4) cannot exist globally in time for large voltage values is more delicate and actually known only in the second-order case when $\beta=\gamma=0$. A classical technique [77] to prove that solutions to the vanishing aspect ratio model (1.9) only exist on a finite time interval for $\lambda$ large is to study the time evolution of

[TABLE]

wher $\phi_{1}$ is a suitably normalized positive eigenfunction of the Laplace operator with homogeneous Dirichlet boundary conditions, and show that $E_{0}$ reaches $-1$ in finite time [42, 61] which corresponds to the occurrence of a touchdown, see Section 5.7 for more precise statements. As described in Section 4.1, a similar approach is used in [93] for the stationary version of (5.2)-(5.4) and heavily relies on the convexity of $u$ with respect to $x$. Unfortunately, such a convexity property is not known for the evolution problem (5.2)-(5.4) (neither it is for (1.9)), the crucial bound (4.2) thus not being available. However, suitable integral estimates on the trace of the gradient of $\psi_{u}(\cdot,u)$ on the upper boundary of $\Omega(u)$ relating it to the behavior of $u$ allow one to study the time evolution of the nonlinear variant

[TABLE]

for which one can derive by Jensens’s inequality an energy inequality of the form

[TABLE]

which has no global solution when $\alpha\in(0,1)$ is suitably chosen and $\lambda$ is sufficiently large. Summarizing, we obtain the following result:

Theorem 5.4** (Finite Time Singularity [30]).**

Let $\beta=\gamma=0$. There is $\lambda_{sing}(\varepsilon)>0$ such that, if $\lambda>\lambda_{sing}(\varepsilon)$, then the maximal interval of existence $[0,T_{m})$ of the solution to the second-order parabolic problem (5.2)-(5.4) provided by Theorem 5.1 is finite, that is, $T_{m}<\infty$. Clearly $\lambda_{sing}(\varepsilon)\geq\lambda_{s}(\varepsilon)$, the latter being defined in Theorem 4.2.

The above mentioned integral estimates on which the proof of Theorem 5.4 relies are also at the heart of the study of non-existence of stationary solutions for the fourth-order free boundary problem as already outlined in Section 4.1.

5.5. The Damping Dominated Limit $\gamma\rightarrow 0$

When damping effects are assumed to be predominant, inertial effects are often neglected from the outset in the modeling. Most research papers hitherto are dedicated to the corresponding “parabolic” version of (5.2)-(5.4) or (1.9) with $\gamma=0$. We now discuss the behavior of the solutions to the “hyperbolic” problem (5.2)-(5.4) when $\gamma>0$ in the damping dominated limit $\gamma\to 0$. The starting point is Theorem 5.1 which provides local solutions to (5.2)-(5.4) in the fourth-order case $\beta>0$ likewise for $\gamma>0$ and $\gamma=0$.

Let us emphasize that not only the unique solution $(u,\psi_{u})$ to (5.2)-(5.4) provided by Theorem 5.1 depends on $\gamma\geq 0$, but also its maximal time of existence $T_{m}$, that is, $(u,\psi_{u})=(u_{\gamma},\psi_{\gamma})$ and $T_{m}=T_{m,\gamma}$. Obviously, considering the time singular limit $\gamma\to 0$ requires a common interval of existence which is independent of $\gamma$, that is, a lower bound $T>0$ on the maximal existence time $T_{m,\gamma}$. Based on this one can show that $(u_{\gamma},\psi_{\gamma})$ converges toward $(u_{0},\psi_{0})$ in a suitable sense as $\gamma\to 0$.

Theorem 5.5** (Damping Dominated Limit [98]).**

Suppose the assumptions of Theorem 5.1 with $\beta>0$ and let $(u_{\gamma},\psi_{\gamma})$ be the unique solution to (5.2)-(5.4) on the maximal interval of existence $[0,T_{m,\gamma})$ provided by Theorem 5.1 for $\gamma\geq 0$. There is $T>0$ such that

[TABLE]

where $\Phi_{u_{\gamma}}$ is the transformed electrostatic potential introduced in (3.5) (with $u$ replaced by $u_{\gamma}$). If $\lambda>0$ is sufficiently small, then $T_{m,\gamma}=\infty$ for any $\gamma$ small enough and the statement (5.6) is true for each $T>0$.

Deriving a lower bound for the maximal existence time independent of the parameter $\gamma$ is the main part of the proof of Theorem 5.5 and relies on an exponential decay of the energy associated with the damped wave equation (an idea inspired by [68]) which is independent of $\gamma$ near zero. The energy functional $\mathcal{E}$ introduced in Section 3.2 then again plays an important rôle as it provides compactness properties of $(u_{\gamma},\Phi_{u_{\gamma}})$ needed in order to obtain a cluster point, its identification as $(u_{0},\Phi_{u_{0}})$ then being due to parabolic singular perturbation theory outlined e.g. in [37, 82]. For further details we refer to [98].

5.6. The Vanishing Aspect Ratio Limit $\varepsilon\rightarrow 0$

A natural and interesting question concerns the relation between the vanishing aspect ratio equation (1.9) and the free boundary problem (5.2)-(5.4). In particular, whether the former can be obtained as a rigorous limit from the latter. In Theorem 4.4 we have seen that this is indeed the case for the stationary problem. We next show that this is true as well for the evolutionary problem and thus give a sound justification of the vanishing aspect ratio model.

Let us recall that the solution $(u,\psi_{u})$ to (5.2)-(5.4) on the maximal interval of existence $[0,T_{m})$ provided by Theorem 5.1 depends on the value of the aspect ratio $\varepsilon>0$ and we shall thus write $(u_{\varepsilon},\psi_{\varepsilon})$ and $T_{m}^{\varepsilon}$ to indicate this. As in the previous study of the damping dominated limit, the starting and simultaneously most difficult point is the derivation of a lower bound $T_{0}>0$ on $T_{m}^{\varepsilon}$ independent of $\varepsilon>0$. This relies on an estimate on the transformed electrostatic potential $\Phi_{u_{\varepsilon}}$ introduced in (3.5). Carefully exploiting the structure of equation (3.4) allows one to control the $H^{2}(\Omega)$-norm of $\Phi_{u_{\varepsilon}}(t,\cdot)$ in terms of $\varepsilon>0$ for any $t\in[0,T_{m}^{\varepsilon})$ for which $u_{\varepsilon}(s)\in S_{2}^{2+2\alpha}(\kappa)$ for $0\leq s\leq t$ (in the notation of Section 5.1). For such $t\in[0,T_{m}^{\varepsilon})$, one deduces the crucial trace estimate

[TABLE]

where $K>0$ depends on $\kappa\in(0,1)$ but not on $\varepsilon>0$. On the one hand, this entails an $H^{2\alpha}(I)$-estimate on the nonlinearity $g_{\varepsilon}(u_{\varepsilon}(t))$ independent of $\varepsilon$ and thus implies that there are $T_{0}>0$ and $\varepsilon_{0}>0$ such that $u_{\varepsilon}(t)$ stays in a set $S_{2}^{2+2\alpha}(\kappa/2)$ for $t\in[0,T_{0}]$, uniformly with respect to $\varepsilon\in(0,\varepsilon_{0})$. This in turn yields compactness of $(u_{\varepsilon})_{\varepsilon\in(0,\varepsilon_{0})}$ and hence a cluster point $u_{0}$ thereof. On the other hand, the trace estimate on $\partial_{\eta}\Phi_{u_{\varepsilon}}$ implies the convergence of the nonlinearity

[TABLE]

which shows that the cluster point $u_{0}$ is indeed a solution to the vanishing aspect ratio equation (1.9). Summarizing, one obtains:

Theorem 5.6** (Vanishing Aspect Ratio Limit [30, 92]).**

Let either $\beta>0$, $\gamma\geq 0$, $\tau\geq 0$ or $\beta=\gamma=0$, $\tau>0$. For $\varepsilon>0$ let $(u_{\varepsilon},\psi_{\varepsilon})$ be the unique solution to (5.2)-(5.4) on the maximal interval of existence $[0,T_{m}^{\varepsilon})$ as provided by Theorem 5.1. There is $T_{0}>0$ such that

[TABLE]

and

[TABLE]

for all $t\in[0,T_{0}]$, where $u_{0}$ is the solution to the vanishing aspect ratio equation (1.9). If $\lambda>0$ is sufficiently small, then $T_{m}^{\varepsilon}=\infty$ for any $\varepsilon$ small, and the statements (5.7)-(5.8) are true for each $T_{0}>0$.

The norm convergence stated in (5.7) can be improved, in particular in the case $\beta>0$ (and also when $\gamma>0$).

5.7. Dynamics in the Vanishing Aspect Ratio Model

As for the stationary vanishing aspect ratio model, we consider a more general geometry and let $D$ be a bounded domain of $\mathbb{R}^{d}$, $d\geq 1$, with sufficiently smooth boundary. The dynamics of the vanishing aspect ratio model

[TABLE]

is better understood, not only because it is simpler than the free boundary problem (5.2)-(5.4) but also as it is a model semilinear evolution equation with a singular reaction term. The study of the latter was actually motivated by the possible occurrence of the quenching phenomenon, which corresponds to a finite time blowup of the reaction term as the solution attains singular values of the nonlinearity, see [40, 99, 100] for instance. In fact, several techniques developed for that purpose have been subsequently applied to (5.9) and its variants. Obviously, the behavior of solutions to (5.9) depends upon whether or not the parameters $\beta$ and $\gamma$ are positive. Before briefly reviewing the available results it is worth pointing out that, in contrast to the free boundary problem (5.2)-(5.4), local existence in a suitable functional setting of a solution to (5.9) up to a maximal existence time $T_{m}>0$ readily follows from the properties of the linear part of the equation by a classical fixed point argument. In particular, given $T>0$, one has

[TABLE]

which provides an explicit criterion for global existence.

Let us now be more specific and begin with the parabolic second-order case

[TABLE]

corresponding to $\gamma=\beta=0$ in (5.9). In that case, the comparison principle represents a powerful tool. In particular, recalling that $\lambda_{0}^{stat}$ is the threshold value of the parameter $\lambda$ above which no stationary solution to (5.11) exists and below which there is a stable stationary solution $U_{\lambda}$, see Section 4.2, we have the following results: if $\lambda<\lambda_{0}^{stat}$, then the negativity of $U_{\lambda}$ and the comparison principle entail that $U_{\lambda}\leq u(t)\leq 0$ for all $t\in[0,T_{m})$. Since $\min\{U_{\lambda}\}>-1$, we readily deduce from (5.10) that $T_{m}=\infty$, see [42, 49]. This result actually implies that the pull-in voltages defined in Section 1.3 satisfy $\lambda_{0}^{stat}\leq\lambda_{0,0}^{dyn}$. Next, when $\lambda>\lambda_{0}^{stat}$, touchdown necessarily occurs and $T_{m}<\infty$ as shown in [42, 49, 61] with the help of a technique originally introduced in [87]. This entails that $\lambda_{0,0}^{dyn}\leq\lambda_{0}^{stat}$ and thus that they are equal as announced in (1.10). The case $\lambda=\lambda_{0}^{stat}$ is also studied in [35, 58] while additional information on the behavior near the touchdown time is provided in [42, 59, 61]. In particular, if $D=\mathbb{B}_{1}$ is the unit ball of $\mathbb{R}^{d}$ and $T_{m}<\infty$, then touchdown necessarily occurs at $x=0$ and only there.

As soon as either $\beta$ or $\gamma$ is not equal to zero, the comparison principle is no longer available for the evolution problem. This leads to a dramatic change in the analysis. Still, taking advantage of the existence of the branch of stable stationary solutions to (5.9) described in Section 4.2 when either $\beta=0$ or $\beta>0$ and $D=\mathbb{B}_{1}$, a technique developed in [87] allows one to prove that touchdown occurs as soon as $\lambda>\lambda_{0}^{stat}$, see [79, 80, 102] when $\beta=0$ and [94] when $\beta>0$ and $D=\mathbb{B}_{1}$. Consequently $\lambda_{0,\gamma}^{dyn}\leq\lambda_{0}^{stat}$, an inequality which is compatible with (1.11) and valid for all $\gamma\geq 0$. As for global existence it is shown in the aforementioned references that global solutions exist provided $\lambda$ is sufficiently small, without any restriction on the domain $D$. But no precise quantitative bound on the threshold value of $\lambda$ is known. We also refer to [109, 112] for numerical simulations and detailed descriptions of the touchdown behavior by asymptotic methods for (5.9) when $\beta>0$ and $\gamma=0$. In particular, there is numerical evidence that, if $d=1$, touchdown takes place simultaneously at two different points for $\lambda$ large enough, a situation which contrasts markedly with the second-order case where such a behavior is only observed with a non-constant permittivity profile [59]. Similar global existence results for small values of $\lambda$ are true for the vanishing aspect ratio model with $\beta>0$ and $\gamma>0$ when the clamped boundary conditions are replaced by pinned boundary conditions $u=\beta\Delta u=0$ on $\partial D$[60]. In fact, considering pinned boundary conditions does not help much for the evolutionary problem as far as the question of global existence is concerned since the comparison principle is not available in both cases, in contrast to the stationary setting described in Section 4.2. Finally, the behavior as $\gamma\to 0$ is analyzed in [101, 102] when $\beta=0$.

To end this overview of results pertaining to the time-dependent vanishing aspect ratio model for MEMS we shortly list some contributions investigating stationary and time-dependent extensions of (5.9), some of them having been introduced in Section 2.4:

• •

Non-constant permittivity profiles: It is already observed in [61, 127] that a non-constant permittivity profile (corresponding to $\chi=0$ and a non-constant $f$ in (2.7)) has a strong influence on the threshold value of $\lambda$ as well as on the touchdown phenomenon. Contributions to this issue may be found in [47, 48, 49, 57, 59, 58, 35, 81, 137] for the second-order case ($\beta=0$, $\tau>0$) and in [13] for the fourth-order equation ($\beta>0$).

• •

Capacitance: The resulting equations are the object of several research papers and we refer to [103, 104, 107, 54, 55, 56, 71, 72, 79, 78, 129, 131, 117] when $f\equiv 1$ and $\chi>0$ in (2.7). The case of a non-constant permittivity profile $f$ with $\chi>0$ is investigated in [14, 15].

• •

Fringing fields: The influence of fringing fields (2.8) has also attracted attention recently and has been analyzed in particular in [130, 27, 115, 116, 138, 136, 44, 125].

• •

Van der Waals or Casimir forces: The inclusion of (2.9) is studied in [89].

• •

When $\beta=0$ nonlinear effects due to large deformations have also been accounted for in the vanishing aspect ratio model (1.9) by replacing the Laplace operator $\Delta u$ by the mean curvature operator $\mathrm{div}\left(\nabla u/\sqrt{1+\varepsilon^{2}|\nabla u|^{2}}\right)$. So far only the stationary problem has been studied, the outcome being that the nonlinear diffusion alters significantly the structure of stationary solutions [9, 10, 8, 17, 124, 125].

• •

Control and inverse problems for MEMS using the vanishing aspect ratio model (1.9) are investigated in [16, 19].

• •

Finally, the extension to $(1+u)^{-p}$ for $p>1$ of the above mentioned results has been studied and shown to depend strongly upon both $p$ and $d$, see [70, 62, 20, 22, 17, 125, 140, 138, 117, 26] and the references therein.

6. The Three-Dimensional Problem

The results presented in Section 4 and Section 5 concern the free boundary problem (2.1)-(2.5) with a two-dimensional region $\Omega(u)$. As noted in Section 2.2 this is a special case of (1.1)-(1.5) when $D$ is assumed to be a rectangle and zero variation in the $y$-direction is presupposed. Of course, a more general (two-dimensional) geometry of $D$ is of great interest with respect to real world applications. But with a general two-dimensional domain $D$, the resulting $\Omega(u)$ is a three-dimensional Lipschitz domain and it is not clear that the trace of $\nabla\psi$ is well-defined on the boundary $\partial\Omega(u)$ thereof. Nevertheless, if $D$ is assumed to be a bounded and convex domain in $\mathbb{R}^{2}$ with a sufficiently smooth boundary, a similar, but technically advanced argumentation as outlined in Section 3.1 (see [90] for details) entails that solutions to the elliptic problem (1.1)-(1.2) are sufficiently smooth so that the trace of $\nabla\psi$ on $\partial\Omega(u)$ has a perfect meaning. However, the regularity of the square of the gradient trace of such solutions is weaker than before and the two-dimensional analogue of the function $g_{\varepsilon}$ from Proposition 3.1 is less regular than stated therein in the sense that it maps $S_{q}^{2}(\kappa)$ with $q>3$ only into $L_{2}(D)$ (instead of some Sobolev space with positive order). As a consequence, the results presently known for the three-dimensional problem are restricted to the fourth-order case $\beta>0$.

So, let $\beta>0$, $\gamma=0$, and consider a bounded and convex domain $D$ in $\mathbb{R}^{2}$ with a smooth boundary. Then one can prove a local existence result for the dynamic problem (1.1)-(1.5) with arbitrary $\lambda>0$ as in Theorem 5.1. Also Theorem 5.3 on global existence of solutions for small $\lambda$ values is still valid. As for stationary solutions, the same statement on existence of such solutions for small $\lambda$ as in Theorem 4.2 is true except for the negativity of $U_{\lambda}$. Nevertheless, this last property holds true if $D$ is a disc so that the maximum principle for radially symmetric functions is available according to Proposition 3.5. In this case, the stationary solutions provided by the analogue of Theorem 4.2 are radially symmetric and the non-existence result of Theorem 4.1 is valid as well. For details on all this we refer to [90].

Let us finally mention that the implicit function theorem is also used in [18] to study stationary solutions in the second-order case $\beta=0$ for small values of $\lambda$.

7. Summary and Open Problems

In the previous sections we reviewed the derivation of the free boundary problem (1.1)-(1.5) as well as the mathematical results obtained so far for its two-dimensional version (2.1)-(2.5). We have in particular shown that it is locally well-posed when either $\beta>0$ and $\gamma\geq 0$ or $\beta=\gamma=0$ and $\tau>0$, having also global solutions and stable stationary solutions for small values of $\lambda$ in both cases (Theorem 5.1, Theorem 5.3, and Theorem 4.2). These results remain true for the nonlinear elasticity model (2.6), see [30, 96], as well as when self-stretching forces are taken into account in addition to bending, that is, $\beta>0$ and the second-order contribution $-\tau\partial_{x}^{2}u$ in (2.3) is replaced by

[TABLE]

The local and global existence results also extend to initial conditions which are not identically zero under suitable smoothness assumptions, possibly supplemented with smallness conditions for global existence. We have also proved that there is no stationary solution for large values of $\lambda$ (Theorem 4.1) while, when $\gamma=\beta=0$, a finite time singularity occurs for the evolution problem (Theorem 5.4). This also holds true for the nonlinear elasticity model (2.6) with the same range of parameters $\gamma=\beta=0$ [29, 30] and extends as well to a broader class of initial data. We summarize the results obtained so far for (2.1)-(2.5) and its stationary counterpart in the table below, assuming that $\tau>0$.

[TABLE]

Though the above table contains quite a number of affirmative answers to the basic questions we have investigated recently, a few of them are, however, still not answered in a satisfactory or complete way. Below we list some open problems which we believe to be of importance for a better understanding of the free boundary model (2.1)-(2.5):

1. (1)

The existence of the threshold value $\lambda_{\varepsilon}^{stat}$ is yet unknown. 2. (2)

By Theorem 4.3 there are countably many values of $\lambda$ in a neighborhood of zero for which there are at least two stationary solutions to (2.1)-(2.5) when $\beta>0$. That this is true for all values of $\lambda$ in a neighborhood of zero would be a first step towards a finer description of the structure of stationary solutions. We actually conjecture that the bifurcation diagram in that case is similar to the one of the vanishing aspect ratio model stated in Theorem 4.5. No clue towards a proof of a similar result when $\beta=0$ is currently available. 3. (3)

According to Theorem 5.4 a finite time singularity occurs for all solutions to the evolutionary problem (2.1)-(2.5) when $\beta=\gamma=0$, $\tau>0$, and $\lambda$ is large enough. It is, however, unknown whether touchdown (in the sense of (5.1)) indeed takes place. 4. (4)

When $\beta>0$, it is yet unproved that a finite time singularity may occur for $\lambda$ large enough. 5. (5)

If touchdown occurs at a finite time $T_{m}$, nothing is known yet about the structure of the touchdown set $\mathcal{S}:=\{x\in(-1,1)\ :\ u(T_{m},x)=-1\}$ and the spatio-temporal scales of the evolution of $(u,\psi_{u})$ as $(t,x)$ approaches $(T_{m},x_{0})$ for $x_{0}\in\mathcal{S}$. 6. (6)

When $\beta>0$, self-stretching forces corresponding to (7.1) received little attention in the literature, even for the vanishing aspect ratio model, and their influence is far from being well understood.