On the Quenching Behaviour of a Semilinear Wave Equation Modelling MEMS Technology

. In this work we study the semilinear wave equation of the form with homogeneous Dirichlet boundary conditions and suitable initial conditions, which, under appropriate circumstances, serves as a model of an idealized electrostatically actuated MEMS device. First we establish local existence of the solutions of the problem for any (cid:21) > 0 : Then we focus on the singular behaviour of the solution, which occurs through ﬁnite-time quenching, i.e. when || u ( · ; t ) || 1 → 1 as t → t (cid:3) − < ∞ , investigating both conditions for quenching and the quenching proﬁle of u: To this end, the non-existence of a regular similarity solution near a quenching point is ﬁrst shown and then a formal asymptotic expansion is used to determine the local form of the solution. Finally, using a ﬁnite diﬀerence scheme, we solve the problem numerically, illustrating the preceding results. where λ is a positive parameter. Problem (1.1) can model the deformation of an elastic membrane inside an idealized electrostatically actuated MEMS.

"MEMS" stands for micro electro-mechanical systems, and refers to precision devices which combine mechanical processes with electrical circuits. MEMS devices range in size from millimetres down to microns, and involve precision mechanical components that can be constructed using semiconductor manufacturing technologies. The devices are widely applied as sensors and have fluid mechanical, optical, radio frequency (RF), data storage, and biotechnology applications. In particular, examples of microdevices of this kind include microphones, temperature sensors, RF switches, resonators, accelerometers, data-storage devices etc., [7,37,41].
The key part of such a MEMS device usually consists of an elastic plate suspended above a rigid ground plate. In the simplest geometry, the elastic plate (or membrane) is rectangular and held fixed at two ends while the other two edges remain free to move. An alternative configuration could entail the plate or membrane (no longer necessarily rectangular) being held fixed around its entire edge. When a potential difference V d is applied between the membrane and the plate, the membrane deflects towards the ground plate. Under the realistic assumption that the width of the gap, between the membrane and the bottom plate, is small compared to the device length, then the deformation of the elastic membrane is described by a dimensionless equation of the form (1.2b) where u = u(x, t) stands for the (dimensionless) deflection of the membrane, ϵ 2 = inertial terms damping terms and λ = Here u 0 (x) and u 1 (x) represent the initial deflection and velocity, respectively, of the elastic membrane. The function f d (x, t) describes the varying dielectric properties of the membrane; for simplicity we assume here that f d (x, t) ≡ 1. Furthermore, T m stands for the tension in the membrane, L c is the width of the parallel plates, each of them denoted by Ω, l c is the unperturbed width of the gap between the membrane and the ground electrode, and ε 0 is the permittivity of free space. The boundary condition represents the membrane being kept in its unperturbed position along its edge. When the damping terms dominate, i.e. when ϵ 2 ≪ 1, then (1.2) reduces to the parabolic problem which has been extensively studied in [7,10,13,21].
(1.3b) In general, the two parallel plates are of arbitrary shape. However, if the parallel plates are thin and narrow homogeneous strips of fixed width L c , see [1,34], then, with suitable scaling, (1.3) can be reduced to the one-dimensional model (1.1). The one-dimensional problem can also be used as a simple model to get better insight into the operation of devices with more general geometries, and especially for the two-dimensional radially symmetric case, i.e. when Ω is a disk, which will be investigated in a forthcoming paper. For certain MEMS-type devices, e.g. resonators and some devices with applications in data storage and optical engineering, [14,39,40,41], the rectangular geometry is practical and it is in fact used. The investigation of the one-dimensional model (1.1) is thus of importance in its own right.
When the MEMS device is connected in series with a voltage and a fixed capacitor one can derive a non-local model of the form where the parameter γ represents the ratio of a fixed capacitance to a reference capacitance, see [10]. Model (1.4), depending on the contribution of the inertial and damping terms, gives rise to the non-local parabolic problem which has been studied in [16,17,19,35], or to the hyperbolic non-local model whose behaviour for the one-dimensional case was investigated in [22]. Recently some authors initiated the investigation of fourth-order models, using the bilaplacian operator which models the moving part of a MEMS device as an elastic plate (with non-zero thickness), rather than as a simple, thin, membrane, [18,25]. There are also papers investigating the quenching behaviour for fourth order parabolic equations, [24,32]. Some recent works investigate the wave equation with damping, [9,30,31].
For a more detailed account of the modelling of MEMS devices, see the books [7,37,41]. From the above, it is clear that the applied voltage V d controls the operation of the MEMS device. Indeed, when V d takes values above a critical threshold V cr , called the pull-in voltage, this can lead to the phenomenon of touch-down (or pull-in instability as it is also known in MEMS literature) when the elastic membrane touches the rigid ground plate, possibly causing destruction of the device in some applications. (The designers of such MEMS devices consequently need to tune the voltage load so that stays away from the pull-in voltage.) Equivalently, this means that there should be some critical value λ cr , depending upon the initial data, of the parameter λ above which singular behaviour should be expected for the solution of problem (1.1). Looking at the nonlinear term of problem (1.1), one can notice that singular behaviour is possible only when u takes the value 1, a phenomenon known in the literature as quenching, see also Section 4. From the point of view of applications it is important to determine whether quenching occurs and, if it does, to clarify when, how and where it might happen. We address two of these questions in this manuscript.
Many authors have investigated the occurrence of quenching for the hyperbolic problem (1.1), [5,6,20,29,38]. However, to the best of our knowledge, the behaviour close to quenching, i.e. the quenching profile, has not been studied previously. In the current work, we first prove some quenching results for problem (1.1) which improve some of the results in [5] for the one-dimensional case and in [29,38] for higher dimensions. Although we mainly focus on the one-dimensional case, our quenching results can be easily extended to higher dimensions as we note in the text.
The outline of the current work is as follows. In Section 2 the local and global existence of solutions to problem (1.1) are studied, while the steady problem is briefly looked at in Section 3. Then, in Section 4, we establish some conditions under which the solution u of (1.1) quenches in finite time, see Theorems 4.2, 4.4 and 4.6. Section 5 is devoted to the investigation of the question of existence of a regular self-similar quenching solution and we finally give a negative answer, see Theorem 5.1. This result is rather surprising since, in standard semilinear wave equations with nonlinearities leading to blow-up, the local behaviour close to blow-up is usually of self-similar type, see [2,3,11,33]. The result is also in contrast to the existence of a self-similar quenching profile for the corresponding parabolic problem, [10]. In Section 6 we go on to use an asymptotic expansion to obtain the local form of the quenching profile as const.×(x − quenching point) 4 3 . Finally, in Section 7, a moving mesh adaptive method is used to obtain a numerical solution of the problem, corroborating the results regarding the quenching profile. We close the paper with a short discussion of our results.
Moreover, the total energy of any weak solution of (1.1) is preserved, i.e.
We notice that Theorem 2.2 implies that the solution of problem (1.1) ceases to exist by quenching. For local existence results in the higher dimensions N = 2, 3, see [29,38]. For a different existence proof see [30].

The Steady-State Problem
The steady-state problem of (1.1) is For the steady problem it is known that there exists a critical value λ * such that problem (3.1) has exactly two solutions (the minimal solution w and the maximal one w) for any λ < λ * , moreover, there is a unique solution 0 < w * < 1 for λ = λ * and no solution for λ > λ * (see [12,27]). We can actually calculate the critical value λ * . If we set W = 1 − w then (3.1) becomes Multiplying both sides of (3.1) by W ′ and integrating from m = min{W (x), ) .
This gives, equivalently, which implies The latter yields, on setting x = 1 so that W = 1, By the above relation we conclude that the maximum of m = m(λ) is attained for λ = λ * ≈ 1.4, see Figure 1.
The computation of the value λ * , with a different way, is also given in [10].

Finite-Time Quenching
By Theorem 2.2 we derive that the solution of (1.1) ceases to exist only when u reaches the value 1 at some point (x, t) ∈ [0, 1] × (0, ∞], i.e. for finite or infinity quenching time. This phenomenon is usually called quenching or touch-down (or pull-in instability) since it corresponds to the situation where the elastic membrane touches down on the rigid plate in the MEMS device. Rather more mathematical discussions of this phenomenon can be found in Sections 5 and 6.
We now present two results regarding the finite-time quenching of solution u(x, t) of (1.1). The first one proves that finite-time quenching occurs when the parameter λ is too big for steady-state solutions to exist. This resembles a result valid for the corresponding parabolic problem, [13,21]. The occurrence of quenching for λ > λ * resembles also the results obtained in [18,25] for the fourth order wave equation.
Proof. For the proof we use the spectral method developed in [23]. We assume that for u(x, t) < 1 almost everywhere in (0, 1) for any 0 < t < T max . (4.1) We first provide some results for the associated linearized eigenvalue problem Set µ 1 = µ 1 (λ; w), the principal eigenvalue of problem (4.2), then µ 1 (λ; w) > 0 and µ 1 (λ; w) < 0 for any 0 < λ < λ * , [12], which indicates the stability of w and the instability of w, see also Figure 1. Moreover, µ 1 → µ * 1 = µ 1 (λ * ; w * ) = 0 as λ → λ * −, as stated in Theorem 1.3 of [12]. Let ϕ * be the eigenfunction corresponding to the eigenvalue µ * 1 = 0, taken to be strictly positive, [12], and normalized so that i.e. ϕ * satisfies (4.5) Now define the functional Multiplying both sides of equation (4.5) with the eigenfunction ϕ * , integrating over the interval [0, 1], using Green's identity and equation (4.4), we then obtain and we then derive, on combining Sobolev's and Poincaré's inequalities (2.2), that where C 0 is a positive constant depending only upon λ and the initial data. Note that due to (4.7) the first term of the right-hand side of (4.6) is estimated from below by On the other hand, the integrand of the second term of the right-hand side of (4.6) is non-negative since Thus, we obtain the differential inequality where and It is readily seen that the positive root of the equation G(t) = 1 is thus lim t→t 1 A(t) = 1− for some t 1 ≤ t + by (4.9). However, the latter, since Theorem 4.2 improves the result of Theorem 3.2 of [5] where quenching was proved only for λ > λ * + , for some λ * + > λ * , and left a gap for the range (λ * , λ * + ]. Moreover, the result of Theorem 4.2 can be easily extended to the practically important two-dimensional case and indeed to three dimensions. The proof follows exactly the same steps. The existence of λ * < ∞ for the higher-dimensional steady-state problem for Ω being a bounded domain of R N , N = 2, 3, is guaranteed by the results in [7,12], where the C 2 −regularity of the extremal solution w * (x) = w(x; λ * ) is also proved. In [7,12], it is additionally proved that the principal eigenvalue of the linearized problem is µ * = 0. Moreover, a lower estimate of the form (4.7) still holds due to Sobolev's inequality holding for N = 2, 3. Therefore the quenching result in [38] can also be improved. In addition, the estimates of the quenching time presented in the next remark are also applicable for N = 2, 3.

Remark 4.3.
The upper bound of the quenching time obtained in Theorem 4.2 can be used to estimate the quenching time, from above, in the asymptotic limit of λ → λ * +.
For u 1 ≥ 0 and not identically zero, so that For λ close to the critical value, any upward perturbation leads to quenching in an order-one time, or less. With u 1 identically zero, so that The middle estimate of the quenching time agrees with one which holds for the corresponding parabolic problem, see [13,23].
We cannot easily get a good bound in the same way in the opposite limit of λ → ∞. This is due to C 0 being potentially unbounded.
We should note, however, that in this one-dimensional case we can proceed slightly differently. Writing We deduce that if u falls to −1 at some time t 1 , quenching must then occur before t 1 + 1.
Combining this result with the estimate got from (4.10) (based on assuming that u remains greater than −1) gives the quenching time estimate t * 1 for λ → ∞.
Since u(x, t) represents the deflection of the elastic membrane inside MEMS device, one expects touch-down to occur when the initial deformation u 0 (x) of the elastic membrane is big enough and/or there is movement towards the rigid plate, meaning that u 1 (x) is positive. This expectation is verified by the following. Proof. Again we proceed as in [23]. Let us assume that the maximum existence time of problem (1.1) is 0 < T max ≤ ∞. For any 0 < λ ≤ λ * set u(x, t; λ) = w(x; λ) + z(x, t; λ), (note that w = w * for λ = λ * ). Then z satisfies (4.11) Let (µ 1 , ϕ 1 ) be the principal eigenpair of the problem: where ϕ 1 is considered to be positive and normalized according to (4.3) and µ 1 is known to be non-negative, since w is unstable, see [7,12].
has non-negative principal eigenvalue as well as a lower estimate of the form (4.7) still valid.
We close this section with a quenching result applying for higher dimensions N > 3, where the estimate (4.7) obtained via Sobolev's inequality is no longer valid. For a similar result see also [10,18].

19) then the solution of problem
Proof. We define the functional where ψ 1 > 0 is the eigenfunction of (4.19) normalized so that ∫ Ω ψ 1 dx = 1. Differentiating F (t) twice and using integration by parts together with equation (4.20a) and Jensen's inequality, yield the differential inequality with associated initial conditions It can be easily seen that for λ > λ * which, by virtue of (4.22) and (4.23), guarantees that (4.24) But relation (4.24) implies that lim t→t 2 F (t) = 1 for some t 2 ≤ t + where Finally, (4.21) implies that ||u(·, t)|| ∞ → 1 as t → t * ≤ t + .

Remark 4.8.
For convenience the proofs of all the quenching results given in the current section concern smooth solutions. However, the same results can be proved for the weak solutions defined by Definition 2.1 under the assumption that u 0 (x), u 1 (x) ∈ L 2 ([0, 1]).

Non-Existence of Regular Similarity Solutions
In many cases the existence of a similarity solution can provide us with a description of the solution profile during quenching. However, as we will see in this section, we do not have such a similarity solution for our problem. Our analysis will also be a guide towards obtaining an asymptotic expansion describing the quenching profile in the following section. For simplicity, we take the quenching time to be t = 0 and position to be x = 1 2 both in this section and in Section 6, provided we consider initial data symmetric with respect to x = 1 2 . Also for simplicity we may consider λ = 1. We take an alternative form of the local hyperbolic problem by setting U = 1 − u. Thus we have with 0 < U < 1 and quenching occurring when U = 0. We set U = (−t) α v(η), for η = (x − 1 2 )/(−t) and we have ∂η/∂t = (x − 1 2 )/t 2 . Then the terms in equation (5.1a) become: To eliminate time we must take α − 2 = −2α or that α = 2 3 and we obtain the relevant equation for v, . Proof. We set v = a + V and then for V = V (η) we get Noting that V = η is a solution of (5.3) if the g term is neglected, we set V = ηq and obtain Using the integating factor η 2 (1 − η 2 ) This gives Then we get

ds.
Thus due to the fact that V = ηq we have In order to obtain regularity at η = 0, with v ′ (0) = V ′ (0) = 0, (i.e. demand the symmetry condition) we have and thus which implies that the solution, V , develops a singularity at η = 1.
Hence any regular symmetric solution must have v(0) > a, i.e. V (0) > 0. From (5.2), it is clear that v is then decreasing for η small and positive. Since v must remain positive, either it must reach a positive local minimum, say v * , at some point η * in (0, 1), or v remains decreasing throughout [0, 1], taking some positive value v 0 at η = 1. We examine the former case first.
As before, Given that V has a minimum at η = η * , V ′ (η * ) = 0 so and hence In particular, Turning again to (5.2), it is seen that for v and V to have local minima at η = η * , v < a and hence V is negative at that point. It follows that 3 g(s)ds ≥ A c + G(η * ) > 0 for η > η * and from (5.4) we again see that V ′ → ∞ as η → 1− so that V develops a singularity at η = 1.
For the other case, if v is to be regular, it will have a first derivative, say v 1 , at η = 1. Then (5.2) gives . We see that for v 0 > a, v 1 > 0 so that v is locally increasing, contradicting the assumption of v decreasing in [0, 1]. Taking v 0 = a, v 1 = 0, and we regain the trivial solution v ≡ a (contradicting v(0) > a). We are left with 0 < v 0 < a and v 1 < 0, so that now v decreases in a neighbourhood of η = 1.
If v is to be smooth, we can differentiate (5.2) to get Supposing that there is a first point η * > 1 where v ′′ = 0, so that v ′′′ ≥ 0 at η = η * , v > 0 (for the solution to still exist) and v ′ < 0 at that point. Then (5.9) gives v ′′′ (η * ) < 0, another contradiction. This means that v ′′ < 0 for η ≥ 1, so that v must fall to 0, and the solution ceases to exist, at a finite value of η.
Remark 5.2. The local behaviour of solutions which are singular at η = 1 can be determined formally. We write η = 1 + σ so that for σ small we are close to η = 1 and the equation has the form We assume that v has the form of a power-series expansion v ∼ v 0 + v 1 σ α + . . . , for some constants v 0 , v 1 and α. Then the equation becomes The leading-order terms, which must balance, are either the first and second, for α ≤ 1, or the third and fourth, for α ≥ 1.
Taking the first two terms to be small gives v 0 = a. This is the special case of v ≡ a. With all the terms of the same size, α = 1 and we obtain, to first-order These are the locally regular, but non-trivial, solutions noted in the above theorem.
We see from (6.4) that the radius of convergence of the power series for p(η) is (in general) 1, consistent with p having a singularity at η = 1. However, the series terminates, with p being a polynomial, p ± n (η), and hence smooth for all η, for α = α ± n with α − n = 2n − 1 for n = 1, 2, 3, . . . and α + n = 2n + 2 3 for n = 0, 1, 2, . . . . For the dominant, slowest decaying, behaviour, with spatial variation, i.e. dependence upon η, we need α = α − 1 = 1 and then a 1 = a 0 /3 and v − 1 = ce −τ (η 2 + 3) on writing a 0 = 3c. On the other hand, a slower shrinking solution still, but without η dependence, is given by α = α + 0 = 2 3 , so that we have v + 1 = de − 2 3 τ , on writing a 0 = d. Therefore we have that the solution to the v problem near the quenching time has the form v ∼ a + ce −τ ( and therefore for U we get approximately for large η as well as large τ , or, in terms of t and x, 2 3 ] .
i.e. v ηη is negligible compared with η 2 v ηη . The neglected term corresponds to the diffusion term, U xx of the original equation U tt = U xx − 1/U 2 . Therefore to determine the outer solution we must solve the equation Multiplying both sides of (6.6) with dU/dt and integrating results in ( dU dt where b is a constant of integration. Therefore we have We set U = 2 b tan 2 θ with dU = 4 b tan(θ) sec 2 (θ)dθ and we obtain In addition we have sec(θ)dθ = ln (tan(θ) + cos(θ)) . Thus .
Finally we have that The quantity inside the logarithm can be written in the following way Therefore we have that
Expanding the quantities in the brackets we obtain ] .
After doing the appropriate eliminations we get, to leading order, that

This implies that
and we obtain an expression for the outer approximation. An alternative way for solving the equation for the outer solution is the following: We have the equation ) . Thus again .
The same leading-order approximation results. Matching. We have the approximation for the inner region being in the form, for v ∼ a + ce −τ ( 2 3 ] . In addition the approximation for the outer region has the form (2 + bU ) while for x → 1 2 and U → 0 we obtain 3c or that the profile of the solution at the quenching point is This gives us an (x− 1 2 ) 4 3 dependence of the solution profile near the quenching point x = 1 2 . Note also that if we rescale to put the factor λ back into the equation, so that (5.1a) is replaced by U tt = U xx − λ/U 2 , then, according to the above analysis, we have that We note that this asymptotic behaviour for the semilinear hyperbolic problem differs substantially from that for the corresponding parabolic problem, see [15] for results specifically for MEMS devices and [8] for general results on the monotonic quenching of solutions of problems which can be written in the form u t − u xx = λ(1 − u) −β . For the parabolic problem, centre manifold techniques have been used in showing (i) that the spatially uniform quenching solution is unstable and (ii) that the quenching profile differs from that suggested by the apparently obvious similarity solution (const.×|x − x * | 2/3 for β = 2 and x * the blow-up point) by a factor of | ln |x − x * || −1/3 . For the hyperbolic PDE, although the simple-minded guess of a self-similar solution again suggests an |x − x * | 2/3 profile, our formal asymptotics now give local behaviour like |x − x * | 4/3 , a quite different power of distance, apparently without logarithmic dependence.
These formal asymptotic results are not yet proved. It is possible that a centre manifoldtype approach might be useful in trying to do so.

Numerical Solution
We now carry out a brief numerical study of problem (1.1), with a variety of initial conditions. A moving mesh adaptive method, based on the techniques suggested in [4], is used. This captures the behaviour of the solution near a singularity.
More specifically we take initially a partition of M + 1 points in [0, 1], 0 = ξ 0 , ξ 0 + δξ = ξ 1 , · · · , ξ M = 1. For the solution u = u(x, t), we introduce a computational coordinate ξ in the interval [0, 1] and we consider the mesh points X i to be the images of the points ξ i (uniform mesh) under the map x(ξ, t) so that X i (t) = x(iδξ, t). Given this transformation, we have, for the approximation of the solution u i (t) ≃ u(X i (t), t), that For the solution of the above system a standard ODE solver can be used such as the matlab function "ode15i", see [4]. In the figures of this section, the results of various numerical simulations are presented. In the numerical method we took M = 161.
The parameter used in the differential equation was λ = 1.5, this value being chosen so as to be slightly greater than the approximate value 1.4 found in Section 3 for λ * , above which quenching should occur by Theorem 4.2.
In Figure 2, u(x, t) is plotted against x ∈ [0, 1] and t ∈ [0, T ] for T = 1.1547. In the next plot, Figure 3, u(x, t) is plotted against x for various times. The uppermost line corresponds to the solution of the problem near quenching. The initial conditions were u 0 = 0 and u 1 = 0.
Similar numerical simulations were carried out for smaller values of λ, still with zero initial data. Taking λ approaching 1.4 from above, the same quenching behaviour was observed. This indicates that, for u 0 = u 1 = 0, λ cr ≈ 1.4. The same was seen to happen for non-zero initial data u 0 lower than w still with u 1 = 0.
From the numerical solution of problem (1.1) we have that near quenching time t * , ln ) ∝ ln (t * − t) with constant of proportionality 2 3 . This is demonstrated in Figure 4, where the fluctuations for t beyond t * − e −30 are apparently due to numerical errors.
A similar plot of ln u(x, t * ) against ln(x − 1 2 ), in Figure 5, shows that u(x, t * ) behaves like C(x − 1 2 ) 4 3 near quenching. The agreement is also illustrated in Figure 6, where the solid line shows the numerical solution of problem (1.1) at the quenching time t * , while the 3 . The constant c = 2.1 is chosen in such a way so that there is agreement of the plots at the boundaries, x = 0, 1.

Similarity Solution.
It is also of some interest to investigate numerically the behaviour of possible similarity solutions, even though we have seen that they do not give local behaviour near quenching. We recall that we have taken 1 − u = U = (−t) α v(η), η = (x −    ) (solid curve) against ln (t * − t) for λ = 1.5. The straight line (dashed) has slope 2 3 and indicates good agreement between 1 − u( 1 2 , t) and const.×(t * − t)   for j = 1, for j = 2 and solution attains a singularity at η = 1, for v(0) = 1 we get the constant solution v ≡ a, while for v(0) = a+1 the solution falls to zero at some η < ∞ -as indicated by the analysis of Section 5.

Discussion
In the current work we have investigated the quenching behaviour of a one-dimensional semilinear wave equation modelling the operation of an electrostatic MEMS device. After establishing local existence, we have proved that the solution u of the equation quenches in finite time, i.e. ||u(·, t)|| ∞ → 1 as t → t * − < ∞ whenever the parameter of the problem λ > λ * , where λ * is the supremum of the spectrum of the associated stationary problem. Although this type of result is fairly standard for related parabolic problems, it is (as far as we are aware) new for nonlinear hyperbolic equations. We also showed that quenching occurs in the parameter range 0 < λ < λ * if the initial conditions are large enough. Similar results were also found for the practically important two-dimensional case, and indeed for the three-dimensional case. Furthermore, in the second part of this work the quenching profile of the solution was studied. In particular, the existence of self-similar solutions was investigated and our main result in this direction was the surprising one that no nonconstant regular self-similar solutions occur. With the aid of this result we studied the profile of the solution near a quenching point and, by use of formal asymptotics, we got that the solution resembles a curve of the form (x − quenching point) 4/3 . Finally, numerical solutions of the problem confirmed the results on the quenching profile.