Curvature-Dependent Electrostatic Field as a Principle for Modelling Membrane-Based MEMS Devices. A Review

The evolution of engineering applications is increasingly shifting towards the embedded nature, resulting in low-cost solutions, micro/nano dimensional and actuators being exploited as fundamental components to connect the physical nature of information with the abstract one, which is represented in the logical form in a machine. In this context, the scientific community has gained interest in modeling membrane Micro-Electro-Mechanical-Systems (MEMS), leading to a wide diffusion on an industrial level owing to their ease of modeling and realization. Physically, once the external voltage is applied, an electrostatic field, orthogonal to the tangent line of the membrane, is established inside the device, producing an electrostatic pressure that acts on the membrane, deforming it. Evidently, the greater the amplitude of the electrostatic field is, the greater the curvature of the membrane. Thus, it seems natural to consider the amplitude of the electrostatic field proportional to the curvature of the membrane. Starting with this principle, the authors are actively involved in developing a second-order semi-linear elliptic model in 1D and 2D geometries, obtaining important results regarding the existence, uniqueness and stability of solutions as well as evaluating the particular operating conditions of use of membrane MEMS devices. In this context, the idea of providing a survey matures to discussing the similarities and differences between the analytical and numerical results in detail, thereby supporting the choice of certain membrane MEMS devices according to the industrial application. Finally, some original results about the stability of the membrane in 2D geometry are presented and discussed.


Introduction
Recent industrial guidelines direct researchers and designers towards the development of low-cost devices to combine physical properties with low-level machine languages. Thus arises the need to design sensors and actuators to meet the multiple requirements of the most widespread industrial, civil and biomedical applications [1][2][3]. In this context, static and dynamic Micro-Electro-Mechanical-Systems (MEMS) technology has matured, especially in domains where miniaturized and integrated electromechanical systems are required [4][5][6]. Moreover, MEMS represent one of the most important achievements of engineering on an industrial scale [5,6]. Currently, the industrial applications of MEMS devices are extremely varied, from applications in the biomedical domain [2,3,7] and thermally driven systems [8,9] to elastic structures [1,10] gaining wide acclaim, is discussed in Section 7, providing food for thought and a comparison between the results obtained by both geometries in order, followed by a discussion of the conditions ensuring the simultaneous existence and uniqueness in both geometries (Section 8). Section 9 deals with the stability and the optimal control problems in 1D geometry, highlighting the range of possible values for V and the limitations concerning the V, which is necessary to win the mechanical inertia of the membrane and the maximum permissible value of V. Moreover, some remarks about the potential energy in 1D geometry are presented and discussed. Some original results of stability and optimal control in 2D geometry are also discussed in Section 10; the relative range of values of V is shown together with the limitations for the minimum V to overcome the membrane inertia and the maximum admissible V. Then, in the same section, relative to 2D geometry, some useful considerations on the values of V maximizing the energy variation are considered and compared with the results obtained in 1D geometry. Following this, numerical approaches for recovering the membrane profile in 1D geometry are discussed in Section 11, providing interesting results in terms of convergence areas and ghost solutions. In Section 12, the most important results obtained by applying numerical techniques on the 2D model are discussed, also highlighting the cases in which, in addition to the non-convergence of the procedure, the phenomena of instability are highlighted. Then, Section 13 offers a comparison between the limitations of the mechanical stress values obtained in both the 1D and 2D formulations. Finally, the conclusion and some future perspectives are presented in the last section of this survey. To facilitate the reading of the survey, the proofs relating to the results known in the literature are reported in the appendices, while proofs of the original results are reported in the text of the survey.

Some Theoretical Backgrounds
The starting point is a well-known dimensionless steady-state model, studied in [5,15]. It considers an MEMS composed of two parallel plates: one fixed and the other deformable, but clumped at the boundary of a region Ω = [−0.5, 0.5] (dimensionless conditions). When a voltage drop is applied, the deformable plate deflects from the rest state (characterized by u(x) = 0) towards the fixed plate (lower plate), located at height h = 1. The profile u(x), in the stationary case, was studied using the well-known fourth-order elliptic model [5,15,16] where f 1 is a bounded function carrying the dielectric properties of the material constituting the deformable plate; λ 1 is a quantity depending on V applied between the plates; ρ, γ and χ are related to the electric and mechanic properties of the material constituting the deformable plate. Moreover, if σ ≤ 2, one can take into account more general Coulomb potential. Moreover, the existence of at least one solution has been studied using the Steklov boundary condition, achieving Dirichlet and Navier boundary ones: u ν represents the outer normal derivative of u on ∂Ω and, ifd = 0, one obtains the Navier boundary conditions, while, ifd = +∞, one obtains the Dirichlet boundary conditions [5,15].

The Cassani-d'O-Ghoussoub Model and Some Theoretical Backgrounds
In R 3 with a system of Cartesian axes O x y z (see Figure 1a), let us consider an electrostatic-elastic system with length 2L, composed by a pair of parallel plates (one fixed and the other one elastic but clumped at the edges). The plates are located at a mutual distance h but lying orthogonal to z . When V is applied, the elastic plate moves towards the fixed plate (for it V = 0). Thus, the electrostatic potential φ satisfies ∆φ = 0 between the plates, such that φ = V on the elastic plate and φ = 0 on the fixed plate. Thus, indicating by ∆ ⊥ the Laplacian operator, with respect to x and y only, the deflection w of the elastic plate satisfies the equation [5]: where T is the mechanical tension of the plate D is the flexural stiffness and 0 is the permittivity of the free space. In (3), − 0 2 |∇φ| 2 represents a source term due to E coupling the solution of the elastic problem to the solution of the electrostatic one [5,16]. Exploiting the following scaling factors (3) becomes the following nonlinear coupled partial differential equation system [5,16]: 2 ⊥ w(x) = −λ 2 2 |∇ ⊥ Φ| 2 + ∂Φ ∂z 2 Φ = 1 on elastic plate, Φ = 0 on fixed plate (5) where the importance of tension and rigidity is δ = D (2L) 2 T and the aspect ratio of the system becomes = h 2L . Finally, setting λ 1 = λ 2 , we write: in which = 0 (2L) 2 2h 3 T considers the electro-mechanical properties of the material constituting the deformable plate. , in dimensionless conditions and by (4) becomes 0 2T > 10 12 . If a membrane replaces the deformable plate, the thickness with D and δ is negligible. Then, as → 0, the first equation of (5) becomes ∂ 2 Φ ∂z 2 = 0 so that Φ = z w from which the second equation of (5) becomes [5,16]: 2 in Ω w(−0.5) = w(0.5) = 1.

Electrostatic and Mechanical Pressures
Once V is applied, an electrostatic pressure, [23,24,27] p el = 1 2 takes place. Thus, if k is a dimensionless coefficient of proportionality, it makes sense to write p = kp el where p is the mechanical pressure which, if its amplitude is sufficient to win the mechanical inertia of the membrane, pushes the membrane towards the fixed plate.
General Formulation of the 1D Model (17) can be written in a more general manner. In fact, being Ω = [−0.5, 0.5] and u : Ω → R, we suppose that u(x) ∈ C 2 (Ω) is the solution of the following general problem (in the Dirichlet form): where f ∈ C 0 (Ω × R × R). Then, considering (17), it follows that f x, u(x), du(x) Apparently, (17) does not present the singularity that characterizes (9). (17) we would achieve d 2 u(x) dx 2 = 0. Thus, considering Remark 4, |E| = 0 produces a linear u(x), which represents a physically unacceptable condition.

Remark 5.
It is worth noting that when the membrane is too thin, the formation of wrinkles around the membrane is highly probable. However, in this work, in the curvature formulation (see (16)) the basic hypothesis is that the membrane profile u(x) ∈ C 2 (Ω). This implies that abrupt local variations of the profile are excluded a priori. In other words, wrinkles around the membrane are not allowed. Obviously, removing the hypothesis that u(x) ∈ C 2 (Ω) will result in formulating the curvature so that any wrinkles around the membrane can be taken into account.

2D Circular Membrane MEMS Device: The Differential Model
Let us consider two parallel disks with radius R, mutual distance h (in dimensionless condition, h = 1) and a circular membrane of the same radius but clumped on the edge of the lower disk, which acts as a support for the membrane. Axial symmetry in the geometry of the membrane inside the device is observed, and considering the z axis a rotation axis, u of the membrane can be thought of as a surface obtained by rotating the curve C, located on the vertical plane rz in the first quadrant and around the axis z when 0 ≤ r ≤ R. Thus, u only depends on the radial coordinate r so that this problem can be considered as a 1D problem wherein the independent variable x is replaced by the radial coordinate r.
E between the disks generates a p el (see (10) in which the coordinate x is replaced by the radial coordinate r), deflecting the membrane. Finally, when the membrane deforms towards the upper disk, the electrostatic capacitance the distance between the membrane and the upper disk varies locally. Therefore, considering the Remark 1 (here is still valid) and just considering the radial part of the Laplace operator [28], model (17) becomes However, (12) holds (x is replaced by r) so that (21) can be written as: Finally, exploiting the expression of mean curvature [25,29] (22) becomes: For details on how to obtain (24), see Appendix A.
Remark 6. Also in this case, as highlighted in Remark 5, in the formulation of the mean curvature (23), the hypothesis that u(x) ∈ C 2 (Ω) remains, so that the possibility of the formation of wrinkles around the membrane is excluded a priori.

General Formulation of the 2D Model
As in the 1D framework, (24) can be written in a more general way, considering Ω = [−0.5, 0.5] and a singularity located at −0.5. Thus, u(r) : (−0.5, 0.5) → R, with u(r) ∈ C 2 (Ω). Then, (24) is a particular case of the general following model: with F ∈ C 0 ((−0.5, 0.5] × R × R) and B, m ∈ R. Moreover, setting B = 0, and m = 0, we achieve (24). The need to rewrite both the problems 1D and 2D ( (17) and (24), respectively) in general terms (as in (20) and (26)) lies in the fact that this generalization allows for using general results on the existence and uniqueness of the solution, as consolidated in the literature, relating to the class of boundary value problems [16,18,20]. Both (17) and (24) do not allow obtaining the solutions explicitly. Therefore, we must be satisfied with obtaining any conditions that ensure the existence and uniqueness of the solution. However, it must be considered that (17) has no obvious singularities, while (24) explicitly manifests a singularity when r = 0 [16,18,20].

On the Existence of at Least One Solution for the 1D Model
As highlighted in [16], it is not possible to obtain the solution in explicit form for (17), so one looks for any conditions ensuring existence and uniqueness. Usually, an important tool is the Banach-Caccioppoli fixed point theorem [28], which guarantees the existence and uniqueness of a fixed point for certain maps of metric spaces on themselves, providing a constructive method to find them. This result has the advantage of proving simultaneously the existence and uniqueness of the solution. However, the conditions to be verified are stringent, and, therefore, for certain Boundary Value Problems (BVPs), such verifications are often impractical. In [16], this theorem was not applicable as an alternative way was chosen, obtaining an important result of the existence of the solution based on the Schauder-Tychonoff fixed point theorem and, sequentially, establishing conditions of uniqueness [16,18,28]. In particular, to achieve a result of existence for (17) by this procedure, it is necessary to start by the definition of two suitable functional spaces [16].
It is well-known that (17) (or its general form (20)), by differentiation, can be transformed into du(s) ds ds (27) where 0 < u < h − d * and G(x, s) is a suitable Green function, the main properties of which are discussed in Appendix B. Then, in [16] the existence of the solution for T(u) = w, with u ∈ P 1 , was proved exploiting the Schauder-Tychonoff fixed point theorem applied to w = T(u) from P to P. The following two results were achieved, the proofs of which are detailed in Appendices C and D, respectively.
then T(u(x)) defined by (27) is an operator from P to P.
Theorem 2. (17) admits at least one solution in P.

On the Existence of at Least One Solution for the 2D Model
In [20], the following result of existence of the solution for (24) was achieved; the proof is presented in Appendix E.
Theorem 3. Let us consider (24) and two twice continuously differentiable functions, u 1 (r) and u 2 (r), defined on [0, 1], with u 1 (r) < u 2 (r) such that for r ∈ (0, 1). Moreover, 1 r du(r) is a continuous function (except for r = 0), satisfying the Lipschitz condition in U × (−∞, +∞), with U = {(r, u) : 0 < r < R and u 1 (r) ≤ u(r) ≤ u 2 (r)}. If with k the constant of proportionality between the displacement of the membrane at the center of the plate u 0 and the mechanical pressure p, there exists at least one solution for (24).
As observed in [20], the greater k is, the lower the value of θλ 2 will be and, thus, dr 2 will be small, so that the concavity of the membrane rises (the greater k 2 is, the greater the influence of p el will be). Then, p will rise increasing the deformation of the membrane. In addition, Figure 3b shows, in the plane d * − θλ 2 , the area of existence of at least one solution: the line of equation θλ 2 = R 2 d * 2 2V 2 0 k in Figure 3b (blue line), separated the area of existence of at least one solution from the area where at least a solution was not ensured.

Conditions Ensuring the Existence of at Least one Solutions for 1D and 2D Geometries: A Comparison
As specified in [16], (28) was the algebraic condition, ensuring the existence of at least one solution for (17). It depended on λ (linked to V) that won the mechanical inertia of the membrane. Moreover, θ took into account the electromechanical properties of the material constituting the membrane. On the other hand, in [20], the specified algebraic condition, (31), ensures the existence of at least one solution for (24). We note that (31) depends on both the mechanical characteristics of the membrane (presence of θ) and V (not only the voltage to win the mechanical inertia of the membrane). In this context, it seems interesting to understand which algebraic condition is weaker with respect to the other one. From (28), we write: and being λ 2 > λ 2 , from (32), it makes sense to write: so that, from both (31) and (33), we achieve The following results holds.
is verified.
Proof. Inequality (35) is immediately verified by substituting the numerical values for each parameter.

Remark 7.
By Proposition 1, the algebraic condition ensuring the existence of at least one solution in 2D geometry is stronger than the algebraic condition in 1D geometry because, in 2D geometry, not only is V considered necessary to overcome the mechanical inertia of the membrane but also the other possible values of V. This is confirmed by the fact that in 1D geometry, V is not explicitly present. Moreover, in 1D geometry, H in (35) exceeds the inf{θλ 2 } with a consequent increase of θλ 2 making the numerical results unrealistic.

On the Uniqueness of the Solution for both 1D and 2D Models
7.1. On the Uniqueness of the Solution for the Model (17) In [16], the uniqueness of the solution for (17) was proved to be always guaranteed as stated by the following result, the proof of which is detailed in Appendix F.  (17) is unique. Moreover, the following properties hold: (2) u(x) is symmetric with respect to the origin; In [16], the uniqueness of the solution for (17), being always guaranteed (see Theorem 4), did not depend on the electromechanical properties of the membrane, unlike the existence of the solution that was conditioned by these properties. If the electromechanical properties of the membrane governed the existence of the solution (i.e., the stiffer the material, the more difficult it is for the membrane to move towards the non-deformed plate and the more the material accumulates electrical charges, the greater the possibility that |E| is more intense), the uniqueness of the membrane deflection was guaranteed regardless of the capacity to accumulate electric charges and the membrane's stiffness. Although this has been proved mathematically in [16], physically it appears to be rather lacking. Versaci et al. in [19] have significantly remedied this gap through an algebraic condition governing the uniqueness of the solution for (17) stronger than the algebraic condition governing the existence. In other words, the existence and uniqueness of the solution for (17), proved in [19], has a more realistic physical-mathematical meaning than the proof given in [16]. In particular, the result obtained in [19] is as follows (the sketch of the proof is shown in Appendix G).

Theorem 5. If
Then, (17) admits the uniqueness of the solution.
In addition, the uniqueness of the solution for (17) depended on the electromechanical parameters of the membrane, but its inertia did not appear. This confirms that when V is applied, the membrane moves, if V overcomes the inertia of the membrane. Thus, the condition of existence of at least one solution depends on λ 2 . Nevertheless, the condition that guarantees the uniqueness of the solution (28) is independent of λ 2 .

On the Uniqueness of the Solution for the Model (24)
Unlike (17), the uniqueness of the solution for (24) was not guaranteed. This important result, studied in detail in [20], is condensed in the following Theorem (the proof is detailed in Appendix H). Theorem 6. Let us consider (24) and suppose that the conditions of the Theorem 3 are satisfied. Moreover, u 1 (r) and u 2 (r) satisfy the given boundary conditions. Then, the uniqueness of the solution u(r), such that u 1 (r) ≤ u(r) ≤ u 2 (r), is not ensured.

1D Geometry
In [16,19], the following results have been proved and Appendix I details the proof. Remark 8. (36), proved in [19], depends on the electromechanical properties of the membrane. In [16] the uniqueness, always guaranteed, did not depend on those properties, since the uniqueness was proved independently of those properties reducing the risk of obtaining ghost solutions. Finally, we note that h − d * does not appear in (36); thus, existence and uniqueness are not dependent on the critical distance. Table 1 summarizes these results.

2D Geometry
Unlike 1D geometry, reference [20] provided an algebraic condition that guarantees the existence of the solution depending both on d * and k. However, uniqueness was not guaranteed. In other words, even if a number of different deflections are allowed, they never reach the upper plate avoiding producing the electrostatic discharge between the two plates. Table 1 summarizes these results.

1D and 2D Geometries: A Comparison
As for the comparison of the algebraic conditions assuring the existence of the solution for both geometries, in this section we deepen which of the two conditions assuring simultaneously the existence and uniqueness of the solution for both geometries is the weakest. From (36) we can write so that, taking into account (31), the following result holds. (31) and (37), it follows that

Proposition 2. From
Proof. As for the Proposition 1, it is sufficient to carry out a dimensional analysis to prove (38).
From Proposition 2 we deduce that the algebraic condition assuring both the existence and the uniqueness of the solution in 2D geometry is more stringent than the one in 1D geometry. Thus, the same observations discussed in Remark 7 continue to apply.

Stability and Optimal Control Problems in 1D Geometry
In Reference [26] authors studied whether the movement of the membrane in 1D geometry, when V is applied, admits stable equilibrium configurations. Furthermore, since the membrane has an inertia while moving and considering that it should not touch the upper plate, in [26], the range of possible values of V in 1D geometry was achieved. Finally, using concepts based on the variation of potential energy stored in the device, optimal control conditions were obtained. In this section, we will present the main results obtained in [26] and some original results regarding 2D geometry.

Stability and Optimal Control in 1D Geometry
In [26], (17) was transformed into its corresponding system of two first order differential equations Particularly, it was set u 1 (x) = u(x) and u 2 (x) = du(x) dx , it was posed [5,26]: in [−0.5, 0.5] with θλ 2 = 0, obtaining the equilibrium point The following result about the stability in 1D geometry was achieved in [26]. The proof exploited the first Lyapunov criterion based on the linearization of (40) in the neighborhood of the equilibrium state. For details, see Appendix J. In [16], it was highlighted that the unique unstable equilibrium point obtained is, in fact, the point considered to be the most critical because it concerns the value of u(x) (located at x = 0) closest to the upper plate of the device. Since V is the main cause of membrane movement towards the upper plate, it is important to know the sup{V}, which ensures that the membrane does not go beyond the unstable equilibrium position. Moreover, having the membrane inertia while moving, in [26], the minimum value of V allowing the deflection of the membrane was obtained knowing the range of possible value of V to understand, by (10), if the p el obtained leads to instability phenomena. Finally, in [26], it was highlighted that the knowledge of the range of possible values of V means to know the range of possible values of λ 2 (see (6)), serving, on the one hand, as a tuning parameter for the device [5] and, on the other, to guarantee the convergence of any numerical approach to recover the profile of the membrane [17,19]. In Sections 9.2-9.4, some important results presented in [26] have been reviewed. Particularly, Section 9.2 discusses the range of possible values for (V min ) inertia to overcome the inertia of the membrane (refer to Appendix K), while Section 9.3 presents some results regarding the (V max ) permissible so that the membrane does not reach the upper plate. Finally, Section 9.4 illustrates the relationship between (V min ) inertia and (V max ) permissible as studied in [26].

V min to Overcome the Inertia of the Membrane in 1D Geometry
(42) has an interesting physical significance. In fact, as T increases (i.e., stretching the membrane more at the edges), a higher value of V is required to overcome the inertia. Furthermore, the greater the distance between the plates, the greater V will be to overcome inertia.

(V max ) permissible in Order That the Membrane Does Not Reach the Upper Plate in 1D Geometry
Proposition 5. Concerning (40), (V max ) permissible is bounded as follows: From (43), it is easy to deduce that the distance between the plates is decisive in the limitation of (V max ) permissible . In fact, increasing h (distance between the plates) requires a greater V to push the membrane towards the upper plate. Furthermore, an increase in k reduces the (V max ) permissible because p generated is higher.

Remark 9.
In [26], to prove (43), indicated by u 0 , the deflection in the center of the membrane, the following algebraic conditions were used.
Please refer to [26] for details of the proofs.
We note that, in [26], an interesting range of possible values for V was achieved. In particular, the two following propositions detail the contents (See Appendix L). Proposition 6. The following inequality holds: Proposition 7. Thus, the range of admissible values for V to win the mechanical inertia of the membrane and to remain far from the upper plate, is as follows: 9.4. Relationship between (V min ) inertia and (V max ) permissible in 1D Geometry Finally, in [26], the relationship between (V min ) inertia and (V max ) permissible in 1D geometry has been obtained as testified by the following Proposition, the proof of which is detailed in Appendix M. Proposition 8. The following inequality holds: Remark 10. (47) suggests us a more manageable algebraic condition. In fact, being h = 10 −9 , 0 = 8.85 × Again, as highlighted in [26], k 2 ≈ 1 since p el p and, moreover, T ≈ 1000 Pa, then (48) can be written as: Thus, (A84) in the proof of Proposition 8 makes physical sense. Moreover, in (A84) , is a constant depending on the geometrical parameters of the device. Therefore, one can write H so that (A84) becomes: Thus, fixing T (i.e., having chosen the material constituting the membrane), depends on G that changes with the geometry of the device.
Conversely, if one chooses the geometry of the device (i.e., G is fixed), the material constituting the membrane determines

Some Remarks about the Potential Energy in the Device in 1D Geometry
As shown in [26], if the membrane is at rest, its distance from the upper plate is equal to h, . Finally, the total variation of the potential energy, ∆W, is: Therefore, from (51), we can write: In [26] was proved that and which combined with each other it follows (see [26]): which in dimensionless conditions became: In [26], the value of the V maximizing ∆W was achieved, as detailed in the following proposition, and the proof is detailed in Appendix N. Figure 4): Remark 11. E, due to the application of V, determines the deflection of the membrane, also determining the amount of ∆W in the device. In such a context, in [26], a limitation for ∆W was obtained starting from |E| as detailed in the following Section.

A Limitation for ∆W Obtained Starting from |E| in 1D Geometry
From both (9) and (17), and considering that [26], proved that In [26], 0 Therefore, in [26], it was proved that condition (52) is stronger than condition (58). In the following section ( Figure 5), some original results on the stability and the optimal control in 2D geometry are discussed in detail.
This is not surprising because 2D geometry is generated by the rotation of a curve lying on a plane. Therefore, for symmetry reasons, the equilibrium positions in 1D and 2D geometries coincide. However, unlike the 1D geometry, here the unique equilibrium point is stable, as detailed in the following proposition.

On the Stability of the Critical Point
Proposition 12. (60) represents a stable equilibrium configuration.
Proof. (59), is writable as: Thus, (59) can be matricially written as:u whereu(r) = [u 1 (r)u 2 (r)] T . To linearize the system, we use the change of variable (A68). Therefore, from (61), considering (A68) and also that u 0 1 and u 0 2 do not depend on r, the following can be stated: From which, developing in Taylor series both f (u 1 (r), u 2 (r)) and g(u 1 (r), u 2 (r)) (neglecting the terms of orders higher than the linear one) and setting τ = ξ 2 + η 2 , it follows that solving which gives ξ(r) = C 1 r C 2 and η(r) = C 3 r , so that ξ(r)(η(r)) C 1 = C 4 , which represents a hyperbola, where C 1 , C 2 , C 3 , and C 4 are all constant. Matricially, (65) is writable as: (66) admits stable equilibrium position only if A does not have eigenvalues with positive real part and if any eigenvalues with real part zero have a unit index. Furthermore, if z 0 = [z 0,1 , z 0,2 ] T , then z(r) = e Ar z(0) = e Ar z 0 . Moreover, |A| = 0, then at least one eigenvalue is zero. Thus, the origin is not an isolated equilibrium point, so that there is at least one line (plane) of equilibrium points: the only critical point obtained is represented by (60). Therefore, as the variation of d * changes the equilibrium point, theoretically, infinite equilibrium points, ∀d * , could take place. In our case, the eigenvalues of matrix A are Moreover, since the number of eigenvalues of A counted with their algebraic multiplicity is equal to the order of A and the geometric multiplicity of each eigenvalue is equal to with the algebraic multiplicity, A is diagonalizable. Thus, e Ar can be written as in which t k and s k are the left and right eigenvectors, corresponding to λ k , respectively. Then, in our case which represents, on the plane z 1 z 2 , an equilateral hyperbola (see the red lines in Figure 6). Furthermore, we achieve z 2 (r) = z 0,1 +rz 0,2 r − z 1 (r) r , which represent straight lines passing through a fixed point (point A in Figure 6), as r → R. Then, as r increases (for r > 0), the hyperbola is traversed such that the point D (see Figure 6 when r = R) tends to the point B (when r = R). Ifż = Az is stable, then the critical point ofu(r) = f(u(r)) is also stable [5,30], and the critical point (60) is an equilibrium stable point for (62).
identifies a point very close to the upper disk (with a high risk of the membrane touching it), it is still a stable point. Electrostatically, it can be justified as follows.
Thereby fixing V, p el ≈ 0 V 2 2d * , so that it does not swing.

Admissible Values for V in 2D Geometry
V min and the Problem to Win the Mechanical Inertia of the Membrane Here, a condition to which (V min ) inertia must satisfy is presented in the following Proposition. Proposition 13. If (31) holds, thus (V min ) inertia to win the membrane mechanical inertia satisfies: Both in 1D and 2D geometry (V min ) inertia are strongly dependent on h, T and θ. The dependence on h is imperative because the greater the distance between the plates, the greater the V to overcome the inertia of the membrane. Furthermore, the greater the T, the greater the mechanical tension of the membrane (i.e., the membrane is tighter), so that greater V is needed to overcome the inertia of the membrane. Finally, if θ is higher, it means that the influence of |E| on the device is greater so that a smaller value of V is sufficient to overcome the membrane inertia (see Remark 1).

(V max ) permissible in 2D Configuration
We present here some useful propositions [20].

Proposition 14.
As in 1D geometry, (44) also holds. It is sufficient to replace x by r.

Remark 12.
Since the proof of Proposition 14 is the same as that given for 1D geometry, (V max ) permissible in 2D geometry is as for 1D geometry (see (43)).
Proposition 17. In 2D geometry, the range of admissible values for V is: Proof. It follows from propositions (13) and (5).
In 1D and 2D geometry, the range of admissible values for V, although different in the formulations (see (46) and (71), respectively), admit the same functional dependence with respect to T, θ and d * . This ensures that the electromechanical properties of the membrane, the critical safety distance and mechanical stress to which the membrane is initially subjected affect the range of possible values for V in the same way for both geometries.
10.5. Relationship between (V min ) inertia and (V max ) permissible in 2D Geometry Proposition 18. In this case, the following inequality holds: Proof. From (43), and considering that d * ≈ 0.1d, it follows .
Therefore, once T is fixed (i.e., the material constituting the membrane has been chosen), On the other hand, once d is fixed (i.e., once B has been chosen), the material constituting the membrane determines Figure 7 left). We note that, in both geometries, the link between (V max ) permissible and (V min ) inertia , in both cases, depends on T (even if with a different functional link). This confirms the fact that whatever the geometry, the larger the T is, the greater the (V min ) inertia .

Some Optimal Control Conditions in 2D Geometry
If the membrane of the device is at rest, the distance between the membrane and the upper disk is d. Therefore, the electrostatic capacitance of the device is C = 0 πR 2 h , so that the potential energy of the device can be evaluated as . If the membrane deforms, is a bounded and continuous function depending on φ. Therefore, the total variation of the potential energy, ∆W = W f inal − W initial becomes , being a bounded and continuous function, allows to write Conversely, from (67), Finally, combining (77) and (78), we achieve the range of the admissible values for ∆W, that is: which can be written as (where F 1 and F 2 constant) analogous to the (56) relating to 1D geometry. Figure 7right depicts the zone of possible values for ∆W, corresponding to the zone below the red straight line and above the blue curve. Obviously, this zone shrinks as k increases.

On the Values of V that Maximize ∆W
Being h−u(r) dφ and considering (81), one obtains C < 4 0 Dd * R 2hd * −k 0 V 2 . Thus, where h(V) is positive, so that the unique stationary point for h(V) is V * = 2hd * (πR−4Dhd * ) 0 πkR . Furthermore, V * is a point of maximum for h(V) (discarding the negative root which represents a non-physical condition to achieve). Furthermore, it is also easy to verify that 10.8. From |E|A to a Useful Limitation for ∆W Starting from (9) and (24), we can write In addition, it follows However, W f inal = 1 2 0 |E| 2 , it makes sense to write and considering that 1 − u − d * < h − d * and du(r) dr < H, with H = 99 positive constant (see [16]), (84) becomes (76) and (85), we can write
Then, by Proposition 19, it follows which represents the limitation for ∆W (depending on θ), obtained from |E|.
The study of analytical models in 1D and 2D geometry [16,26] has determined algebraic conditions, ensuring the existence and uniqueness for the solution for both models. Then, by means of suitable numerical procedures, solutions are obtained that, if they satisfy the aforementioned algebraic conditions of existence and uniqueness, do not represent ghost solutions.

Zeros of F(η): The Dekker-Brent Approach
This approach, used in [22], utilizes the bisection procedure to solve a non-linear equation. For each iteration, b k approximates temporary zero; a k represents the "contra-point" such that F(a k ) and F(b k ) have opposite signs so that [a 0 , b 0 ] contains the solution; b k−1 represents b at the previous iteration. Therefore, two temporary values are evaluated: the first is obtained by the secant method and the second by the bisection approach; otherwise. If s, as result of the secant procedure, falls in (b k , m), then s = b k+1 , otherwise m = b k+1 . Thus, the new value of the contra-point is chosen so that F(a k+1 ) and F(b k+1 ) have a different sign, so that a k+1 = a k , otherwise, a k+1 = b k . Finally, if |F(a k+1 )| < |F(b k+1 )|, a k+1 is a best approximation of the solution, so that a k+1 and b k+1 are exchanged. Sometimes, b k converges slowly, so Brent proposes a modification of this approach using a test that must be satisfied before the result of the secant method is accepted for the next iteration [23]. If δ is a tolerance and if the previous step has been used in the bisection procedure, |δ| < |b k − b k−1 | and |δ| < |s − b k | < 1 2 |b k − b k−1 | have to be applied, otherwise the bisection procedure is used again. If the previous step uses interpolation, |δ| < |b k−1 − b k−2 | and |δ| < |s − b k | < 1 2 |b k−1 − b k−2 | are applied to decide if to perform the interpolation or the bisection. The Brent approach, ensuring that at the kth iteration the bisection method is used at most for 2 log 2 times, utilizes inverse quadratic interpolation instead of a linear one (as in the secant procedure).

Obtaining the Solution
η k , at each iteration, was evaluated in [17,19,22] by solving the related Initial Value Problem (IVP). Thus, a suitable stop criterion is used to verify if η k → η as k → ∞. The solutions are achieved exploiting both ode23 and ode45 MatLab ® R2017a routines (accuracy and adaptivity parameters defined by default). However, the main difficulty in achieving the solutions is related to the integration of unstable initial value problems: thus, the solutions of the BVP could be insensitive from the variations of the boundary values. However, the solutions of the IVP achieved by the shooting procedure are computed by the variations of the initial values [22].
where p(r) is an interpolation polynomial of degree lower than s interpoling [x n,i , F(x n,i ), y(x n,i )], i = 1, 2, ..., s, x n,i = x n + τ i h, i = 1, ..., s, 0 ≤ τ 1 < ... < τ s ≤ 1. If the Lagrange method is used, thus: where L j (r) are the fundamental Lagrange polynomials. Therefore, plugging (93) into (92), one achieves so that, (94) is forced for all the x n,j , so that u n,j at collocation node points are obtained, for i = 1, ..., s, by If τ s = 1, then y n+1 = y n,s , otherwise y n+1 = y n + ∑ s j=1 F(x n,j , y n,j ) x n+1 x n L j (r)dr. Collocation methods are reliable tools, although may not be suitable if high accuracy is required [23].
Particularly, (97) means that is exact for all polynomials in which the degree is lower than p (if (97) is satisfied, then the procedure has quadrature of order q). Analogously for condition (98): if it is satisfied, then t+c i h x F(s)ds ≈ h ∑ s j=1 a ij F(x + c j h) are exact for all polynomials in which the degrees are lower than q. It is worth noting that all methods satisfying condition (98) having c i , i = 1, 2, ..., s distinct are collocation procedures [30].

The Three-Stage Lobatto IIIa Formula
This procedure requires that c i be chosen as roots of [30] s is the number of the stage), achieving that c 1 = 0 and c s = 1 ∀s. Therefore, the quadrature formula is exact for any polynomial in which the degree is less than 2s − 2 [32]. Two definitions are necessary for the following.

Definition 2.
Let us define the mesh-grid: defining, on it, the step-size h m = r m+1 − r m .

Remark 13. p(r) satisfy the boundary conditions in
This approach can be considered as a collocation procedure and it is equivalent to the three-stage Lobatto IIIa implicit RK procedure [32], the Butcher tableau of which is [30] Thus, the three-stage Lobatto IIIa formula becomes [30]: and Remark 14. Again, this procedure is achievable from (91) by the Simpson quadrature formula to approximate the integral between x n and x. Obviously, when the procedure is applied to a quadrature problem, it reduces (103) to the well-known Simpson formula [31]: which are non-linear equations solvable by a MatLab ® routine. In addition, MatLab ® , ∀r ∈ (a, b), evaluates the cubic polynomial by means of its special routine bvpval [30].

Remark 16.
As known, a BVP could have more than one solution [31]. Thus, it is important to give an initial guess for both the initial mesh and the solution. Obviously, the MatLab ® solver makes the mesh adaptively achieve a solution using a reduced number of mesh points [30].
However, a good initial hypothesis could be very difficult. Thus, the MatLab ® solver checks a residue defined as [30] res(r) = p (r) − F[r, p(r)]. Obviously, if res(r) is small, then p(r) is a good solution. Moreover, the case of the well-conditioned problem, p(r) is next to y(r). In [23], MatLab ® R2017a bvp4c solver has been exploited since it implements the collocation technique by a piecewise cubic p(r), with coefficients determined requiring that p(r) be continuous on (a, b). Furthermore, both mesh and estimation error are based on the computation of the residual of p(r), the control of which is exploited to manage inadequate guesses for both mesh and solution [30]. Moreover, this routine provides a very reduced computational complexity to compute the Jacobian J = ∂F i ∂y . Finally, bvp4c is a vectorized solver, so that it is able to strongly reduce the run-time vectorizing F(r, y(r)) [31].

Four-Stage Lobatto IIIa Formula
It is derived as an implicit RK method. Its Butcher tableau is [30]: Like the three-stage formula, this approach is also a polynomial collocation procedure, but it provides solutions in which the accuracy is of the fifth-order, belonging to C 1 ([a, b]). However, unlike bvp4c that exploits the analytical condensation procedure, MatLab ® solves the four-stage Lobatto IIIA formula by a finite difference procedure (by means of bvp5c solver) and solves the algebraic equations directly. Furthermore, unlike bvp4c that handles the unknown parameters directly, bvp5c augments the system with trivial differential equations for unknown parameters [30,31]. Figure 8 left depicts the numerical results for u(x) evaluated in [22] exploiting different values of θλ 2 using the shooting procedure implemented by the MatLab ode23 routine. It can be noted that the minimum value of θλ 2 ensuring convergence is 0.63. Similar results have been achieved exploiting the other numerical procedures. For the shooting procedure and ODE solvers (indicated by Shoot&23 and Shoot&45), it was set u 1 (0) = 1 and u 2 (0) = 1.2 as initial guess when θλ 2 = 0, 63, 1, 4 and u 1 (0) = 0.1 and u 2 (0) = 0.2 as initial guess when θλ 2 = 2, 3. Regarding the relaxation procedure & Keller-Box scheme (Rel&Box), both initial guesses were set as u 1 (0) = u 2 (0) = 1 [22]. Finally, u 1 (0) = u 2 (0) = 0 were set for the collocation procedure and Lobatto formulae (indicated by Col&III and Col&IV). Table 2 presents a comparison of the results achieved when θλ 2 = 0.63, 1, 2, 3, 4. Finally, when θλ 2 = 4, the same value max(u(x)) = 0.029918, with J = 13, J = 52, J = 4000, J = 4 number grid points, respectively, was achieved. The results proved that each procedure showed a good performance, even if different orders of accuracy and different grid points were used. Noting that both the relaxation method and Keller-Box scheme reveal robustness and accuracy, the latter provides results as accurate as those of the shooting and collocation method, since it involves more grid points for each iteration. However, although the shooting method is not so robust as the relaxation and collocation procedures, it is faster and well-implemented in MatLab. Moreover, the collocation procedure gives a solution by a few numbers of grid points because the profile is smooth. Finally, the Relaxation procedure and the Keller-Box scheme are more robust with respect to the other approaches.

Convergence of the Numerical Approaches
In [22], indicated by [(θλ 2 ) conv ] S ode23 , the range of θλ 2 that guarantees convergence by the shooting method (using ode23 MatLab ® ) routine, it was obtained that [(θλ 2 ) conv ] S ode23 = [0.63, +∞), so that, if [(θλ 2 ) no conv ] S ode23 = [0, 0.63), the convergence is not guaranteed. Moreover, using the Keller-Box scheme, it was found that the range of θλ 2 ensuring convergence was [(θλ 2 ) conv ] S Keller−Box = [0.592, +∞). As mentioned above, of [(θλ 2 ) no conv ] Keller−Box = [0, 0.592) the convergence of the Keller-Box scheme is not guaranteed. Moreover, the range of θλ 2 ensuring convergence when the shooting procedure was implemented by the ode45 MatLab ® routine was [(θλ 2 ) conv ] S ode45 = [0.63, +∞), so that the range that did not guarantee convergence was [(θλ 2 ) no conv ] S ode45 = [0, 0.63). Finally, exploiting both III and IV Lobatto IIIA formulas, in [22], it was found that [(θλ 2 ) conv ] Three Stage Lobatto = [1.181, +∞) and [(θλ 2 ) conv ] Four Stage Lobatto = [1.181, +∞). Then, for IV Stage Lobatto IIIa formulas, it was obtained that the range of θλ 2 that did not ensure convergence was [(θλ 2 ) no conv ] Four Stage Lobatto = [0, 1.181). These conditions are summarized in Table 3. If all the procedures work parallelly, the minimum value of θλ 2 ensuring the convergence of at least one numerical procedure was obtained as follows [22]: In other words, for values greater than 0.592, the convergence of at least one numerical solution is ensured. Then, for θλ 2 ≥ 0.63 convergence is ensured for all the numerical procedures considered. However, even if a numerical solution is obtainable, one must be sure that this solution does not represent a ghost solution.

Convergence and Ghost Solutions
As studied in [22], from (36), one obtains θλ 2 ≥ 18 so that if θλ 2 ∈ [18, +∞) we have that both existence and uniqueness are ensured. Moreover, from (6) [16], we can write λ 2 = 0.25 0 V 2 Moreover, combining (36) and (108), one obtains 1 + sup du(x) Being 0.63 18 1 + sup dy(x) dx 6 we obtain 0.63 (110) highlights that the thicker the membrane, the higher V to be applied to the device for overcoming the inertia of the membrane. Moreover, since 18 18 1 + sup dy(x) dx 6 , one can write: Thus, both (110) and (111) identify, on the plane formed by T and V, zones of convergence where ghost solutions could also take place. As depicted in Figure 8 right [22], (110), below the blue curve, the non-convergence area is identified while, between the blue and red curves, at least one numerical procedure converges (ghost solutions could take place). Finally, above the red curve, the area in which convergence is guaranteed without ghost solutions is identified. To achieve solutions that are not ghost solutions is very important because it gives us the possibility to achieve ranges of possible values for V, |E| and T defining operating conditions to which the device has been subjected [17,19,22].

Range of Parameters for the Correct Use of the Device
In this section, we discuss the correct use of the device as studied in [22]. In other words, once the material constituting the membrane has been chosen (that is, T is fixed), what is the range of possible values for V and |E|? Vice versa, fixing both V and |E| (that is, fixing the intended use of the device), which material is more suitable for the membrane? Thus, starting from both (36) and (108), one obtains giving the interval of admissible values for sup , once θ, Tand Vare known. Moreover Multiplying (112) by λ 2 and considering that (1 − u(x)) 2 < 1, |E| 2 < sup{|E| 2 , β 1 = 0 /2T and (1−u(x)) 2 from which and By (115), fixing θ and T, one obtains V 2 sup{|E 2 |} (operative electrostatic parameters of the device). Vice versa, one obtains so that, starting from |E| and V, T and θ are achieved.

On the Applicability of the Numerical Procedure
As is known, BVPs are much more difficult to solve than IVPs [33]. Unlike the IVPs (having a unique solution), a BVP could not have a solution, could have a finite number or could have infinitely many. To solve BVPs in 2D geometry, shooting procedures and a collocation one can be exploited. The first one combines a numerical procedure based on the solution of a corresponding IVP for ordinary differential equations and one for the solution of non-linear algebraic equations. However, the shooting procedure could require the integration of an unstable IVP, so that the solution of a BVP could be insensible to the changes in the boundary value and the solution of the IVP could be sensible to the changes in the initial values. On the other hand, the collocation procedure, even if they are efficient and reliable tools, often could not be suitable for high accuracy. In [22,23], the authors used the collocation procedure based on piecewise polynomial functions for solving (24) because just one coefficient was singular and, moreover, the solutions were smooth [33]. Thus, they proposed the collocation procedure implemented in MatLab ® , by its routine bvp4c, because its code implements the three-stage Lobatto IIIa formula, which is a collocation formula exploiting a collocation polynomial, which provides a C 1 continuous solution that is fourth-order accurate.
12.2. Numerical Procedure and Convergence: Interesting Ranges of θλ 2 and V 2 k

θλ 2 and Its Characteristics Ranges
As specified in [23], instability phenomena of the membrane can arise when V grows too much. Thus, it is important to know the range of V generating instability. With V being linked to θλ 2 (see (31)), it follows that the knowledge of the behavior of the membrane when θλ 2 increases gives us, when both d * and k are fixed, the range of V producing instability of the membrane in absence/presence of ghost solutions. In [23] it has been observed that θλ 2 , depending on both the electromechanical properties of the material constituting the membrane and V, once the ranges of stability/instability of the membrane are known, it is possible to know the operation parameters in the convergence area respecting (31) and the engineering areas of applicability of the device. In such a context, it was reasonable to consider that p el and p are equivalent (i.e., k = 1 and negligible losses). This is correct because when V is applied |E| and p el are generated inside the device. In [23] all the simulations have been carried out by the bvp4c MatLab ® solver exploiting both the default relative and absolute error tolerances obtaining 100 as the optimal number of grid points on [0, R] = [0, 1] (by a greater number of grid points the performance did not improve). Moreover, different solutions starting with different initial guesses were obtained. Three cases occurred.
Case 1 ∀θλ 2 ∈ (10 −6 , +∞) the numerical procedure converged without instabilities. In other words, when θλ 2 increased from 10 −6 , d 2 (u) dr 2 increased from negative values towards zero. Thus, the concavity of the deformed membrane decreased avoiding instabilities next to the edge. This was confirmed by (31): the higher θλ 2 is, the lower V 2 will be (increasing θλ 2 , the membrane does not deform too much when a low V is applied). Figure 9, as achieved in [23], displays an example of recovering with θλ 2 = 0.5, and initial guesses u 1 ≤ 2.446 and u 2 = 0. With V being reduced, the membrane moves just a little so that instabilities do not appear. The procedure behaves differently when the initial guess of u 1 increases (with u 2 = 0). In fact, Figures 10-12 [19.979, +∞), respectively (initial guess for u 2 is zero).    One can note that when the initial guess for u 1 and u 2 = 0 increases, the recovering of the membrane is symmetrical but erratic ( Figure 10) until the profile assumes a bell-shape ( Figure 11). The erratic behavior is also present when the initial guess is increasing ( Figure 12). However, although Figures 10-12 show simulations that numerically are valid, since u 1 is greater than d, they are not realistic.

No Convergence Convergence and Instability Convergence and Stability
Moreover, λ 2 depends on both V and the electromechanical properties of the material constituting the membrane (see (6) and (7)). Then, combining (121) with (6) one achieves [23]: Again However, in dimensionless conditions, (h − u(r)) 2 < 1, d = R = 1 and |E| 2 < sup{|E| 2 }; thus, from (123), one can write θλ 2 = 4T 2 sup{|E| 2 } . Moreover, as known, θλ 2 is writable as [17,20] If the numerical procedure does not converge, it follows that 10 −7 > θλ 2 > 2 0 V 4 T 2 sup{|E| 2 } so that: Then, from (124), one obtains: Once the intended use of the device has been chosen (that is, once the pair {V, sup{|E| 2 } has been fixed), from the inequality (125), inf{T} is computable. In other words, once {V, sup{|E| 2 } has been chosen, the material of the membrane is selected. Conversely, if the material of the membrane has been chosen (that is, T has been fixed), it is possible to achieve the pair {V, sup{|E| 2 } satisfying the inequality (126) (that is the intended use of the device is achieved).

Electromechanical Properties of the Membrane and Exploitation of the Device: A Useful Comparison between the 1D and 2D Formulations in Convergence Conditions
The following result holds.
From the Proposition 20 it follows that, with the same V and sup{|E| 2 } (i.e., fixed the intended use of the product), the machine voltage T is higher in 2D geometry. This is due to the fact that, in 2D geometry, it is necessary to take into account all the contributions due to the mechanical stresses relative to each vertical plane passing through the vertical axis of symmetry. Recall that each of these planes is affected by 1D geometry.

Conclusion and Perspectives
In recent literature, the surveys published on MEMS membrane electrostatic devices are abundant and many of them are of high quality. The discussions contained therein range from the design procedures of the individual devices to the techniques for analyzing the behavior of the devices under the most varied operating conditions. However, recently, a new line of research has emerged on MEMS membrane devices based on the observation that on each point of the membrane E is always orthogonal to the tangent to the membrane at the point in question, so that |E| is to be considered proportional to the curvature K of the membrane. This approach allows modeling the problem by means of elliptic semi-linear second-order differential models relative to 1D and 2D geometries (the latter with circular symmetry). Particularly for the 1D geometry, the explicit singularities present in models known in the literature are not present, while for the 2D geometry, the singularity is present on the axis of symmetry passing through the center of the circular plates of the device. In this survey, in addition to having presented and discussed both the differential models mentioned above, the main results of existence and uniqueness of the solution for both models were presented, discussed and compared in the first part, verifying which conditions are stronger than the remaining. The second part of the survey was entirely devoted to stability and optimal control problems for both geometries, underlining the similarities and differences that emerge from the study of both geometries. In particular, important conditions that manage the range of possible values for V have been studied and compared providing original results as far as the 2D geometry is concerned. The last part of this survey, since the analytical models do not provide explicit solutions, was entirely dedicated to the comparison of the results obtained from the numerical procedures applied to reconstruct the membrane profile. Particularly, respecting the analytical conditions of existence and uniqueness of the solution and for both geometries studied, the ranges of admissible values for T to which the membrane must be subjected before deformation under the effect of V were compared. The comparisons highlighted how both geometries provide limitations and a range of possible values for the fundamental quantities (such as, for example, V and T) dependent on the same parameters but with different functional bonds due to the different membrane shapes in the two geometries. This, at least qualitatively, confirms that both models studied have been correctly formulated, highlighting a certain uniformity of behavior of the control parameters. Finally, we observe that, even if the approach of considering |E| is proportional to K, as validated by the results obtained in the absence of ghost solutions, it appears clear that the formulations used for K require refinement; therefore, in the future, it would be desirable to exploit more sophisticated formulations for membrane curvature in order to take into account the symmetries present in the geometries studied. Moreover, by removing the hypothesis that u(x) ∈ C 2 (Ω), new scenarios open up regarding the reformulation of the curvature, thus taking into account any formation of wrinkles around the membrane. It is worth underlining the fact that the numerical models used enjoy a limited computational load. Then, the writing of the hardware of what has been implemented in the software appears desirable in order to make the elaboration usable for any real-time applications.
Author Contributions: All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

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

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Mean Curvature and 2D Modellization
Here, a surface S generated by rotating around z a curve C lying on a vertical plane normal to the xy plane, which forms an angle t with the zx plane, is considered (Figure 3a) [29]. To simplify the procedure, let us suppose that C is parametrized with r (generic parameter). Thus, P(r) = ( f (r), 0, g(r)), r ∈ I ⊂ R, in which f (r) and g(r) are regular functions, which satisfy the following inequality Remark A1. P(r), being a so-called natural parametrization, ensures that the curve C is regular everywhere. Thus, by rotation, S is also regular.
Thus, one can easily achieve from which, the coefficients of the first fundamental form become: Observing that F = 0 everywhere, thus the coordinate lines are orthogonal to each other (everywhere), so that, one can write: Therefore, from (A2), Moreover, the coefficients of the second fundamental form become: To obtain both the principal curvatures, k 1 (t, r) and k 2 (t, r), one needs to solve (e − kE) and k 2 (t, r) = d 2 f (r) r 2 . Then, the mean curvature, dr d 2 g(r) r 2 and since C belongs to the plane y = 0, one can set: in which both . Therefore, by H(r), (22) becomes: Finally, one can divide both sides of the equation of (A7) by , so that eqrefmedia can be written as (24).

Appendix B
As known, differentiating (20), one obtains: where 0 < u < h − d * and G(x, s) is a suitable Green's function [34]. Moreover, The existence of the solution of the equation T(u) = w, where u ∈ P 1 , can be proved by the Schauder-Tychonoff fixed point theorem applied to w = T(u) from P to P [34]. In fact, from (A10), one defines the positive operator T as follows: The suitable Green's function is [34] G(x, s) = (s+0. Due to the construction of G(x, s), it follows that T(u(−0.5)) = T(u(0.5)) = 0 verifying that T(u) ≥ 0. Moreover, it is easy to verify that µ(x, u(x)) > 1 in [−0.5, 0.5]. This is due to the fact that E deforms the membrane satisfying the condition that locally it must win the inertion of the deformation, so that µ(x, u(x)) assumes values greater than the unity. Therefore, to overcome the inertia of the membrane, a minimum voltage λ > 0 should be applied, so that λ < λ < sup{λ} and sup{λ} being a bounded quantity, it follows that 1/λ 2 < +∞. Therefore, also considering the (A13), one can write: Moreover, we can deduce that: To prove that T(u) ∈ P, one considers that 4(h − d * ) 1  Therefore, it follows that θλ 2 > Hθλ 2 2(h−d * ) from which: which represents the condition guaranteeing that T(u) : P → P.

Appendix D. Proof of Theorem 2
By Theorem 1 and taking into account that the following compact immersions C 2 0 [−0.5, 0.5] → C 1 0 [−0.5, 0.5] and P 1 → P holds, applying the Schauder-Tychonoff fixed-point theorem, u(x) = T(w(x)) admits at least a fixed point u(x) = T(u(x)) in P 1 . In other words, there exists at least a solution for (17).

Appendix E. Proof of Theorem 3
Let us introduce the following Lemma.

Appendix H. Proof of Theorem 6
We premise the following Lemma.

Appendix I. Proof of Theorem 7
As seen, (28) guarantees that (17) admits at least one solution. Moreover, from Theorem 4, it follows that the solution is also unique. Therefore, to ensure both existence and uniqueness, it is imperative to solve the system: . Therefore, H should be lower than a notably small amount that is incompatible with the definition of H [19]. Thus, (A63) is equivalent to (36).
Remark A3. (A70) makes sense because u 1 and u 2 , as proved in [16], are analytical functions allowing the linearization by Lemma A5. In a system as (A73), if there is at least one eigenvalue with a positive real part or a not simple eigenvalue with zero real part, then the system is unstable [26].
Proof. It follows from Lemma A6.
Remark A4. The membrane, when V is applied, deforms reaching a distance, in u(0), equal to h − d * . However, this equilibrium is an unstable configuration.
Furthermore, (A83) can be written as: Appendix N. Proof of Proposition 9 . f (V), for usual values of d, d * , k and 0 , and considering (46), is a continuous function. In fact, if in f (V), 2dd * − k 0 V 2 = 0, we would avchieve V = 2dd * k 0 ≈ 0.1503V. However, this value is lower that inf(V min ) inertia , thus V ≈ 0.1503V is not within the range of possible values for V. To calculate the stationary points for f (V), we set . The negative root must be discarded, since this value of V would deform the membrane symmetrically with respect to the lower plate (condition physically impossible to obtain).