Interfacial instability of thin ﬁlms in soft microﬂuidic conﬁgurations actuated by electro-osmotic ﬂow

We analyze the interfacial instability of a thin ﬁlm conﬁned between a rigid surface and a pre-stretched elastic sheet, triggered by non-uniform electro-osmotic ﬂow. We derive a nonlinear viscous − elastic equation governing the deformation of the elastic sheet, describing the balance between viscous resistance, the dielectric and electro-osmotic eﬀects, and the restoring eﬀect of elasticity. Our theoretical analysis, validated by numerical simulations, shows several distinct modes of instability depending on the electro-osmotic pattern, controlled by a non-dimensional parameter representing the ratio of electro-osmotic to elastic forces. We consider several limiting cases and present approximate asymptotic expressions predicting the electric ﬁeld required for triggering of the instability. Through dynamic numerical simulations of the governing equation, we study the hysteresis of the system and show that the instability can result in an asymmetric deformation pattern, even for symmetric actuation. Finally, we validate our theoretical model with ﬁnite-element simulations of the two-way coupled Navier equations for the elastic solid, the unsteady Stokes equations for the ﬂuid, and the Laplace equation for the electric potential, showing very good agreement. The mechanism illustrated in this work, together with the provided analysis, may be useful in toward the implementation of instability-based soft actuators for lab-on-a-chip and soft-robotic applications.


I. INTRODUCTION
Interfacial instabilities of thin liquid films with a free-surface subjected to temperature gradients or chemical gradients have been extensively studied for over six decades [1,2]. These Marangoni instabilities are induced by forces acting on the liquid−fluid interface and do not exist in fluidic systems with solely no-slip liquid−solid interfaces. In recent years there has been a growing interest in similar type of problems involving elastic sheets suspended on thin liquid films, where elastic effects determine the evolution of the interface geometry. In particular, much attention was given to investigate the fluid and solid mechanical instabilities in blistering, which forms when an injected viscous fluid peels an elastic sheet from a solid surface [3]. Few viscous−elastic instabilities of the fluid−elastic interface have been reported to date, including snap-through instabilities induced by the interaction between a buckled elastic arch and viscous flow [4], as well as wrinkling of lubricated elastic sheets under compression [5][6][7]. However, given the similarity between free-surface thin-film equations and those describing the dynamics of a lubricated elastic sheet, one would expect similar (closely analogous to Marangoni) interfacial instabilities to occur in elastic-plated thin films and no such studies have been presented to date.
Configurations involving viscous flows bounded by elastic structures are relevant to a wide spectrum of applications such as fabrication of flexible microelectro-mechanical systems [8,9], suppression of viscous fingering instabilities [10][11][12], impact mitigation [13], fabrication of microfluidic devices [14][15][16][17][18][19] and soft robotics [20][21][22][23]. In particular, Inamdar and Christov [24] studied the transient fluid−structure interaction in a two-dimensional elastic micro-channel and developed a one-dimensional lubrication model, which accounts for bending and nonlinear induced tension, as well as the inertia of solid and liquid. Meanwhile, Martínez-Calvo et al. [25] extended the steady analysis of Christov et al. [19] for a slender geometry to the transient case by accounting for fluid and solid inertia in the lubrication and Kirchhoff−Love equations, respectively.
Viscous flows bounded by elastic substrates are also often encountered in lab-on-a-chip and microfluidic devices, in which electro-osmotic flow (EOF) is a commonly used driving mechanism. EOF is the bulk fluid motion arising from the interaction of an externally applied electric field with the net charge at a solid−liquid interface. In previous works, we have suggested the use of non-uniform electro-osmotic flow as an actuation mechanism to create desired dynamic deformations in a lubricated elastic sheet by inducing internal pressure gradients in the fluid [26,27]. Considering small deformations and strong pre-stretching of the lubricated elastic sheet, we examined the linear and weakly nonlinear viscous−elastic interaction driven by non-uniform EOF, which exhibited stable behavior [27]. Since the pressure, formed due to EOF, scales inversely with the thickness of the liquid film, sufficiently negative pressures can TF IG. 1. Schematic illustration of the examined configuration, showing the coordinate system and relevant physical and geometric parameters. A thin viscous liquid film of initial thicknessh0 is confined between a rigid surface and a pre-stretched elastic sheet of lengthlm and thicknesshm, supported at its boundaries. Non-uniform electro-osmotic flow, induced by an EOF slip velocityũEOF(x,t), creates negative pressures within the viscous fluid, resulting in viscous−elastic interaction, which leads to deformationd(x,t) of the elastic sheet. Above a certain threshold of the electric field, the system exhibits instability which collapses the elastic sheet onto the floor.
trigger instability of the liquid−elastic interface, which acts to diminish the thickness of the film. We have recently demonstrated this concept for the simplified case of a plate−spring model, accounting only for temporal dynamics between two rigid plates [28].
In this theoretical work, we analyze the complete nonlinear viscous−elastic interaction in the case of large deformations and examine the spatiotemporal evolution of the interfacial instability of a pre-stretched elastic sheet subjected to non-uniform EOF. In Sec. II, we present the problem formulation and the equations governing the viscous−elastic dynamics for constant voltage and constant current actuation modes. We provide their scaling and summarize the key non-dimensional parameters and assumptions used in the derivation of the model. Focusing on the case of a constant applied current, in Sec. III we present a linear stability analysis of the system and further provide an analytic expression of the threshold electric field for the onset of instability in a gravity-dominant regime. Using dynamic numerical simulations, in Sec. IV we consider both a constant voltage and a constant current actuation modes and examine the interfacial instability under various physical conditions, showing the existence of hysteresis for the onset of instability. We further explore the effect of bending on the onset of instability and provide closed-form expressions for the threshold electric field in tension-versus bending-dominant regimes. In Sec. V, we demonstrate that the system may exhibit distinct modes of instability depending on the magnitude and the spatial form of the electro-osmotic pattern. Specifically, we show that the instability can result in an asymmetric deformation pattern, even for a symmetric actuation. In Sec. VI, we perform finite-element numerical simulations to validate the results of our theoretical model and show a very good agreement between the two. We conclude with a discussion of the results in Sec. VII.

II. PROBLEM FORMULATION AND GOVERNING EQUATIONS
We study the viscous−elastic dynamics and interfacial instability of a thin liquid film subjected to non-uniform EOF and confined between a flat rigid surface and a pre-stretched elastic sheet of lengthl m , thicknessh m , densitỹ ρ m , Young's modulusẼ Y , and Poisson's ratio ν. Figure 1 presents a schematic illustration of the configuration and the Cartesian coordinate system (x,z), whosex axis lies at the lower flat surface andz is perpendicular thereto. We denote dimensional variables by tildes, normalized variables without tildes and characteristic values by an asterisk superscript.
The fluid density and viscosity are respectivelyρ andμ, the fluid velocity isũ = (ũ,w) and fluid pressure isp. The total gap between the plates ish(x,t) =h 0 +d(x,t), wheret is time andh 0 is the initial gap. The deformation field d(x,t) is induced due to an internal pressure formed by a non-uniform and time-varying electro-osmotic slip velocitỹ u EOF (x,t) on the rigid surface. The characteristic velocities in thex andẑ directions are respectivelyũ * andw * , and the characteristic pressure, deformation, and time are respectively denoted asp * ,d * , andt * . We assume that surface roughness prevents complete contact between the elastic sheet and the bottom surface, and denote this by a minimal gaph r . Our assumption that the elastic sheet remains wetted even when in contact with the surface is similar to the pre-wetting modeling introduced in work by [9].

A. Governing equations in dimensional form
Considering a shallow fluid layer and negligible fluidic inertia, represented by small reduced Reynolds number, the fluid motion is governed by the lubrication equations [29] ∂ũ ∂x whereg is the acceleration of gravity acting in −ẑ direction. These equations are subjected to the electro-osmotic slip and the no-penetration boundary conditions at the bottom surface, as well as the no-slip and the kinematic boundary conditions at the fluid−elastic interface, whereũ EOF (x,t) is the electro-osmotic slip velocity, which in the thin-double-layer limit is given by the Helmholtz− Smoluchowski equation [30],ũ whereinε is the fluid permittivity,ζ is the zeta potential, andẼ is the applied electric field. We examine configurations where the electric fieldẼ(x,t) is created either by applying a constant electric currentĨ or a constant voltage drop V. For comparison between the two actuation modes, we consider that initially they both induce the same electric fieldẼ 0 . The explicit expressions of the electric fieldẼ(x,t) for both actuation modes arẽ with details of the derivation provided in appendix A. We note that in the lubrication limit of shallow configurations, the electric field is independent of thez coordinate to the leading order in [31]. Following standard lubrication theory, from (2) and (3), the evolution of the fluid−elastic interface,h(x,t), with respect to the rigid substrate is related to the fluidic pressure,p(x,z,t), by the Reynolds equation (see Ref. [29]), where the last term represents spatial variations in electro-osmotic fluxhũ EOF , which drives the fluid−structure interaction. We neglect the weight and the inertia of the elastic sheet, γ =ρ mhml 2 m /Tt * 2 1, and focus on the case of a strongly pre-stretched elastic sheet with tensionT , assumed to be much larger than any internal tensionT in ∼ (d * /l m ) 2Ẽ Yhm [32] formed in the sheet during the deflection, α = (d * /l m ) 2Ẽ Yhm /T 1. Upon application of the electric field, a pressure is induced on the elastic sheet by the non-uniform electro-osmotic flow as well as by direct traction by the Maxwell stresses. We assume that the permittivity of the elastic sheet and air are negligible compared to that of the fluid and therefore neglect their contributions to the Maxwell stress, accounting only for a dielectric contribution of the fluid. Based on these assumptions, and under the assumption of small slopes, |∂d/∂x| ∼d * /l m 1, the pressure in the fluidp(x,z,t) can be expressed by a combination of elastic bending and tension stresses, and Maxwell stresses and the hydrostatic pressure [32][33][34], whereB =Ẽ Yh 3 m /12(1 − ν 2 ) is the bending stiffness, whereinẼ Y andh m are assumed to be constants, and the last term is an upward-directed dielectric contribution arising from the Maxwell stresses.
Combining (6) and (7) yields the nonlinear governing equation for the deformation Invoking current conservation and electroneutrality in the bulk fluid yields ∂(h(x,t)Ẽ(x,t))/∂x = 0 (see appendix A) and thus the last term on the left-hand side of (8) can be expressed as where the last equality stems from (5). Using (4)-(5) and (9), (8) takes the form Equation (10) is a nonlinear viscous−elastic governing equation, which accounts for bending, tension, gravitational and dielectric effects and describes the deformation field due to non-uniform EOF, acting either in a constant current or a constant voltage actuation mode. In addition, (10) clearly indicates that the driving mechanism for viscous−elastic interaction is non-uniform EOF due to heterogeneous zeta-potential distribution. This is in contrast to the study of Mukherjee et al. [18], which considered the relaxation of an initially deformed micro-channel under a uniform zeta potential.
The governing equation (10) holds only forh >h r . In regions whereh =h r , additional contact forces come into play, preventing penetration of the elastic sheet onto the surface, and (10) is no longer valid. Since these forces are not a priori known, we simply apply a kinematic condition ∂d/∂t = 0 in these regions. We emphasize that for the purpose of calculating the pressure, (7) can be used only in regions where the elastic sheet is free of contact with the surface, whereas (6) can be used uniformly everywhere, once the deformation field is obtained.

B. Scaling analysis and non-dimensionalization
Scaling by the characteristic dimensions, we introduce the non-dimensional variables where the characteristic velocity in thex direction,ũ * , is given by the Helmholtz−Smoluchowski slip condition as u * = −εζ * |Ẽ 0 |/μ, whereinζ * is the characteristic negative value of the zeta potential, so thatũ * is positive. From order-of-magnitude analysis of the continuity andx-component momentum equations (2), we obtainw * = ũ * and p * = 12μũ * / 2l m . We note that as the fluid motion is driven by EOF through the electro-osmotic slip velocity, the characteristic pressure is independent of the viscosity,p * = −12εζ * |Ẽ 0 |/ 2l m [35]. Since the deformations we are interested in here are on the order of the initial fluid thickness, it is convenient to scale the deformation byh 0 (d * =h 0 ), so that the fluid layer thickness h can be expressed as h = 1 + d.
In this study, our main focus is on a tension-dominant regime, and therefore the appropriate scaling for the viscous−elastic time scale is based on tension and is obtained by balancing the first and the second term on the left-hand side of (8), yieldingt * = 12μl 4 m Th 3 It is worth noting that analogous expressions for a viscous−elastic time scale were previously obtained by Elbaz and Gat [36] for the case of viscous fluid flow in an elastic cylinder and by Martínez-Calvo et al. [25] for the case of a start-up flow through a deformable micro-channel.
C. Viscous−elastic governing equations for constant current and constant voltage actuation modes Substituting (11) and (12) into (10), we obtain the non-dimensional viscous−elastic governing equation for the deformation where we have introduced the function F(t) and the non-dimensional parameter, which we refer as an elasto-electro-osmotic number E EOF , indicating the relative contribution of electro-osmotic and elastic restoring forces to the deformation of elastic sheet. We note that an elasto-electro-osmotic number can be either positive or negative due to the sign ofẼ 0 and is similar to the capillary number encountered in free-surface thin-film flows [1]. Three additional positive non-dimensional parameters appear in (13). The first two parameters, B and G, determine the relative importance of bending and gravity versus tension forces, respectively, and are defined as The last non-dimensional parameter ϕ appearing in (13) is defined as [28] ϕ =Th and is independent of the applied electric field in contrast to E EOF . We note that the product ϕE EOF , given by −Ẽ 0h 2 0 /(12ζ * l m ), is a non-dimensional parameter that represents the ratio of the dielectric to electro-osmotic effects. Table 2 lists the physical parameters for a typical microfluidic configuration withh 0 = 100 µm andl m = 5 mm, showing that ϕ is of O(10 −3 − 10 −2 ). Moreover, since ϕ and E EOF scale as (h 0 /l m ) 4 and (h 0 /l m ) −3 , ϕ can attain much smaller values for more shallow typical configurations, while keeping ϕE EOF 1, corresponding to the negligible contribution of the dielectric forces. Therefore, in this work, we restrict our analysis to the case of ϕ 1 (more strictly speaking ϕ = 0) and neglect the contribution of the dielectric forces. We note that for very large electric fields, the dielectric contribution becomes apparent and its effect can not be neglected in the analysis. Further investigation would be required to access the effect of this contribution on the interfacial instability.
We consider the following boundary conditions at the edges of the elastic sheet and the initial condition d(x, t = 0) = 0. The first two conditions correspond to no deflection and no moment at the boundaries, whereas the last condition is obtained from (7) by assuming the fluidic pressure has a zero gauge value at the boundaries, p(x = 0, 1, z = 1, t) = 0. The corresponding flow field can be described using the stream function ψ, given by and related to the velocity field throughũ = (∂ψ/∂z, −∂ψ/∂x). For h > h r , the pressure gradient ∂p/∂x can be calculated either using (6) or (7), whereas when the film thickness reduces to h = h r , the elastic balance (7) is no longer valid and ∂p/∂x is obtained from (6).
For completeness, we here provide a list of the non-dimensional numbers used in the problem, where the parameters in (20a) correspond to the shallowness of the fluid layer as well as the relative importance of the fluid inertia, the solid inertia, the internal tension and the dielectric contribution, wherein our assumptions also require 1, Re 1, γ 1, α 1 and ϕ 1. The non-dimensional numbers appearing in (20b) determine the relative importance of bending and gravity as well as the relative magnitude of electro-osmotic forcing.
In this work, we examine the evolution of the deformation field and the interfacial instability, using a particular case of a spatially non-uniform cosine zeta-potential distribution as a test case, and explore the effect of the wavenumber k on the resulting deformation and the onset of instability.

III. LINEAR STABILITY ANALYSIS FOR THE CASE OF CONSTANT CURRENT
We here focus on the case of a constant current actuation mode (equation (13) with F(t) = 1) and examine the linear stability of the corresponding steady-state solutions. We consider small perturbations of the sheet from its equilibrium deflection d ss (x) by letting where f (x)e σt is the disturbance of the deformation from its steady state, σ is the growth rate of the perturbation and s 1. Substituting (22) into (13), at the leading order we obtain the equation for steady-state deformation, At the first order in s , we find that the eigenfunction f (x) satisfies the eigenvalue problem, from which the growth rate σ of the perturbation, being the eigenvalue, can be evaluated. In the following sections we determine the steady-state deformation d ss (x) and the corresponding eigenfunctions f with eigenvalues σ by solving numerically the steady-state boundary value problem (23) and the corresponding eigenvalue problem (24) subjected to the boundary conditions (18). Additional details of the numerical method are provided in appendix B.

A. Gravity-dominant regime
In this section we consider a gravity-dominant regime, where the hydrostatic pressure dominates over the tension and bending stresses, and obtain an analytical expression for the threshold elasto-electro-osmotic number, E EOF,CR , corresponding to the growth rate of σ = 0 which determines neutral stability.
Assuming that the elastic bending is small compared to the tension, G 1 B, we first neglect both bending and tension terms in (23) and obtain subjected to the boundary conditions where the subscript G denotes the solution obtained solely from the gravitational contribution. For the case of a cosine zeta-potential distribution (21), a closed-form analytical solution of (25)−(26) for the steady-state deformation is given by We note that all other roots are not physically relevant (one solution gives always negative film thickness, while the other two are complex). It follows from (27) that a solution exists provided that and thus the threshold value of the elasto-electro-osmotic number is Each such critical elasto-electro-osmotic number, E EOF,CR , corresponds also to a maximum deformation, d max,CR , beyond which the system becomes unstable. Here, d max,CR = −1 for G 1 B. As expected, since we have neglected the highest derivatives in (23), the solution (27) cannot satisfy all boundary conditions along x = 0 or x = 1. To explore the effect of neglecting the tension term on the deformation and the onset of instability, we keep both tension and hydrostatic terms and solve numerically (23) with B = 0 and k = 1, subjected to the first four boundary conditions in (18). Figure 2(a) presents the steady-state deformation along thex axis for various values of E EOF , with G = 100. Solid lines represent the closed-form asymptotic solution (27), whereas dashed lines represent the numerical solution which takes into account both tension and hydrostatic terms. Similar to the results shown in the study of Tan et al. [37] on steady thermocapillary flows driven by non-uniform heating, in the case of pure hydrostatic stress there is a cusp at d(1/2) = d max,CR = −1 in the shape of the elastic sheet for E EOF = E EOF,CR,G that disappears when the external tension is included. Furthermore, accounting for tension results in a smoothing effect on the shape and in a reduction of deformation magnitude. Figure 2(b) presents the growth rate σ as a function of E EOF /|E EOF,CR,G | for both cases, and clearly indicates that the steady-state deformations shown in Fig. 2(a) are stable, since the negative growth rate tends to zero only as E EOF approaches E EOF,CR,G , given by (29). As can be inferred from the results of Fig. 2, while neglecting tension in the gravity-dominant regime results in an overestimated maximum deformation, it accurately predicts the threshold value of the elasto-electro-osmotic number at which the system is at neutral stability.
In appendix C, we provide the variation of the growth rate σ with E EOF for higher values of wavenumber k in a tension-dominant regime. We show that for k ≥ 2 the growth rate approaches zero both for positive and negative values of E EOF , implying that the instability will occur regardless of the direction of the applied electric field. We also present and discuss the variation of the growth rate σ with k 2 for several values of E EOF , showing that σ approaches a constant value of −π 4 as k increases regardless of the value of E EOF .   1)), we obtain that in this case, the threshold elasto-electro-osmotic number scales as E EOF,CR ∼ a 1 + a 2 G, where a 1 and a 2 are constants. Using the asymptotic limits G 1 and G 1, we determine the coefficients a 1 and a 2 for the case of k = 1. For G 1, numerical results show that the threshold elasto-electro-osmotic number E EOF,CR is almost independent of G and thus we solve (23) with B = G = 0 and obtain E EOF,CR = a 1 = −6.51. On the other hand, for G 1 the threshold elasto-electro-osmotic number scales linearly with G and is accurately predicted by the asymptotic solution (29), i.e. a 2 = −πG/4, illustrated in Fig. 3(a) as a gray solid line. Therefore, the threshold elasto-electro-osmotic number E EOF,CR scales as represented in Fig. 3(a) by a solid black line, showing good agreement with numerical results (black dots). Figure 3(b) presents the corresponding critical maximum deformation d max,CR as a function of G. Comparing the results of the Fig. 3(b) to the pure hydrostatic case shows that tension yields a reduction of the threshold value of maximum deformation for the onset of instability. The critical maximum deformation d max,CR attains approximately a constant value up to G = O(1) and then for G > 1 monotonically increases with G throughout the investigated range.

IV. DYNAMIC SIMULATIONS
To investigate the viscous−elastic dynamics and the spatiotemporal development of the instability, we solve numerically the nonlinear evolution equation (13) using finite differences. Further details of the numerical procedure are presented in appendix B. In this section, we focus on the case where the established electro-osmotic flow drives the fluid from the center of the system outward, corresponding to a zeta-potential distribution described by (21) with k = 1. In all numerical simulations, hereafter we set h r = 10 −2 .

A. Deformations due to constant current and constant voltage actuation modes
In Fig. 4 we show the deformation resulting either from a constant current ( Fig. 4(a,c,e)) or a constant voltage ( Fig. 4(b,d,f)  a plate−spring model [28]. However, the present model allows for the first time to observe the spatial behavior of the elastic sheet as it approaches, and finally contacts the bottom surface. Figures 4(e) and 4(f) present the shape of the steady-state deformation along thex axis for several values of E EOF . The numerical analysis reveals several differences between the two actuation modes. Firstly, the critical maximum deformation d max,CR , below which the system is unstable, as well as the corresponding magnitude of threshold value E EOF,CR are smaller for a constant current (d max,CR = −0.44; E EOF,CR = −6.51) than for a constant voltage (d max,CR = −0.57; E EOF,CR = −10.1). Secondly, we observe that the part of the elastic sheet where h = h r , which we denote by ∆x and refer as the contact width, is significantly larger for the case of a constant current as compared to the case of a constant voltage. These differences stem from the source term in (13), (14), which remains constant for the case of applied current, but deceases as the elastic sheet descends for the case of applied voltage. In contrast to the case of a constant voltage where the pulling effects of EOF and restoring effects of elasticity weaken as h decreases, being nonlinearly coupled through F, in the case of a constant current the electro-osmotic forcing remains constant, while the restoring effects of elasticity decrease with h. Therefore, for the case of a constant current, the electro-osmotic driving force overcomes the restoring effects of elasticity at smaller absolute values of E EOF,CR and d max,CR , thus triggering an instability. This means that triggering the instability in the case of a constant voltage requires a comparatively higher initial electric field, which will then decease due to increase in electric resistance of the configuration as the elastic sheet is pulled downward.

B. Investigation of contact width ∆x
To highlight the effect of the elasto-electro-osmotic number on the transient behavior and the magnitude of the contact width, ∆x, we consider for simplicity the case of a strongly pre-stretched elastic membrane (B = G = 0), actuated by a constant voltage, and solve numerically the governing equation (13) for a wide range of E EOF values.
As an illustrative example, Fig. 5 the dynamic nature of ∆x. After the onset of instability, as the film thickness reduces to h = h r , the contact width ∆x increases until the system reaches a steady state, where the electro-osmotic forcing is balanced by the restoring effects of elasticity as well as the additional contact forces that come into play. Figure 5(b) presents the contact width ∆x as a function of |E EOF |, above the threshold value E EOF,CR = −10.1. As expected, the contact width ∆x monotonically increases with |E EOF |, yet this dependence weakens as |E EOF | increases.

C. Hysteresis for the onset of instability
In the previous sections we showed the existence of the onset of an interfacial fluid−elastic instability above a certain threshold value of |E EOF,CR |. In addition to instability, the nonlinearity of the system that arises from the inverse dependence of the EOF force on the film thickness, results in hysteresis. In this section, we show that the transition between stable and unstable states may therefore occur at different values of E EOF,CR depending on the current state of the elastic sheet. For example, Fig. 6 illustrates that after the onset of instability, the elastic sheet may remain adhered to the floor after a decrease of electric field up to a different critical value of |E EOF,CR |. Figure 6(a,b) presents the maximum deformation d max at steady state as a function of E EOF for constant current and constant   (31) and (32). Gray and black symbols correspond to constant current and constant voltage actuation modes, respectively. All calculations were performed using k = 1 and G = 0.
voltage actuation modes, showing a hysteresis loop. Black dots represent the maximum deformation of an initially flat elastic sheet which starts to descend in response to an applied electric field. As the electric field is increased above the threshold value, the instability occurs and the elastic sheet approaches the bottom floor. However, if at this point we decrease the magnitude of electric field, the elastic sheet does not ascend but remains in contact with the floor until significant reduction of the applied forcing is reached, below which the sheet rises and achieves a stable non-contacting steady-state position, represented by gray dots.
To quantify the effect of gravity on the hysteresis, Fig. 6(c,d) presents the magnitude of the threshold elasto-electroosmotic number E EOF,CR as a function of G for constant current and constant voltage actuation modes, resulting from these two initial states of the system. For G 1, the magnitude of E EOF,CR corresponding to the adhered initial state is smaller than the magnitude of E EOF,CR corresponding to the undeformed initial state, and both elasto-electro-osmotic numbers are independent of G. However for G 1 (G > 100), the two threshold elasto-electro-osmotic numbers become identical and, as expected, scale linearly with G.

D. Effect of bending on the onset of instability
In Secs. IV A -IV C, we neglected the influence of bending and gravitational effects and considered a membrane (tension-dominant) regime with B = G = 0. Aiming to elucidate the effect of bending on the onset of instability, in this section we consider a finite value of B and obtain the threshold elasto-electro-osmotic number, E EOF,CR , and the corresponding deformation, d max,CR , by solving numerically the sixth-order governing equation (13). We determine the threshold value E EOF,CR using a bisection method and restricting our resolution up to two decimal places. In addition, for simplicity, we eliminate the effect of gravity by setting G = 0. Figure 7(a) presents the magnitude of the threshold elasto-electro-osmotic number E EOF,CR as a function of B, for the cases of a constant current and a constant voltage. We note that the presented results and behavior for E EOF,CR , which correspond to tension-bending regime, are qualitatively similar to the results shown in Fig. 3(a) for tension-gravity regime. Similarly to the tension-gravity regime, performing a scaling analysis of (13) (with x ∼ 1 and d ss = O(1)), we obtain that the threshold elasto-electro-osmotic number E EOF,CR scales as E EOF,CR ∼ a 1 + a 2 B, where a 1 and a 2 are constants. To find the value of a 1 , corresponding to the case of B 1, we first neglected the bending and gravity contributions and then solved the resulting fourth-order nonlinear equation, yielding a 1 = −6.51 represented in Fig. 7(a)   numerical results (gray and black dots). We chose to present the results on a log−log plot in order to verify that there are no singularities as B approaches zero. Figure 7(b) presents the corresponding critical maximum deformation d max,CR as a function of B, for a constant current (gray dots) and constant voltage (black dots) actuation modes. As opposed to the tension-gravity regime considered in Sec. III, where for G > 1 d max,CR monotonically increases with G (see Fig. 3(b)), for the tension-bending regime considered here, the resulting threshold value of maximum deformation d max,CR attains approximately a constant value, indicating remarkably weak dependence on B throughout the investigated range.

V. SYMMETRIC AND ASYMMETRIC DEFORMATION PATTERNS RESULTING FROM A SYMMETRIC ACTUATION
Aiming to examine more complex deformation patterns, we consider the case of a strongly pre-stretched elastic sheet (B = G = 0) actuated by a constant voltage, and focus on the zeta-potential distribution (21) with k = 3. While for k = 1 the instability may occur only for negative values of E EOF , for k = 3 (more generally for k ≥ 2) the system may exhibit the instability both for positive and negative values of E EOF , since negative pressures may still arise. Figure 8 illustrates several distinct modes of instability depending on the sign and the magnitude of E EOF (or an applied electric field), and indicates that the instability can result in an asymmetric deformation, even for a symmetric actuation. Figure 8(a,b) presents the time evolution of the deformation field and the corresponding pressure distribution for E EOF = 32. For positive values of E EOF , an initially sinusoidal-shaped deformation transits to first-mode deformation behavior, illustrated in Fig. 5(a), due to a high negative gauge pressure that develops at the center of the system.
On the other hand, for negative values of E EOF , corresponding to Figs. 8(c,d) and 8(e,f), we observe a completely different deformation behavior after the onset of instability at E EOF,CR = −153.3 is reached. Clearly, while the baseline system is symmetric, our numerical results reveal that there is a range of E EOF , −254 ≤ E EOF ≤ −153.3, for which small non-uniformities in the system, here simulated by numeric round-off errors, result in more rapid growth of the instability on one side, which changes the system in such a way that prevents the collapse on the other side. We note that small changes in the grid result in inversion of the instability to the other side, indicating that our numerical solver is not biased in one direction. Figure 8(c,d) presents the time evolution of the deformation field and the corresponding pressure distribution for E EOF = −200, showing that initially sinusoidal and symmetric deformation and pressure fields take an asymmetric form after the onset of instability. However, as illustrated in Fig.  8(e,f), as the magnitude of E EOF increases above the limiting value E EOF = −254, the system transitions back from the asymmetric behavior to the symmetric one, characterized by a symmetric deformation pattern with two contact regions where h = h r .
Figures 9(a) and 9(b) present the resulting streamlines at steady state underneath the deformed elastic sheet for E EOF = −200 and E EOF = −275, obtained from (19) using (6). As opposed to the case of a symmetric deformation, ( Fig. 9(b)) where the corresponding volume flux vanishes at steady state, the case of an asymmetric deformation ( Fig. 9(a)) is characterized by a negative net flux from right to left.

VI. FINITE-ELEMENT NUMERICAL VALIDATION
To validate the results of our theoretical model, we performed finite-element numerical simulations with the commercial software COMSOL Multiphysics (version 5.0, COMSOL AB, Stockholm, Sweden). Complete details regarding the governing equations, boundary conditions, domain discretization, and physical parameters employed in the finite-element numerical simulations are provided in appendix D. Figure 10 presents a comparison of finite-element simulation results and theoretical model predictions for the case of a constant applied voltage. Figure 10(a) presents the maximum deformation at steady state as a function of E EOF , showing very good agreement between the theoretical predictions (black dots) and the finite-element simulations (red crosses). The threshold value of E EOF for the instability is also well predicted, with a value of E EOF,CR,Th = −11.2 versus E EOF,CR,FE = −11.3 for the theoretical model and finite-element simulations, respectively. Figure 10 Fig. 10(b,c,d) that at early times, corresponding to relatively small deformations, there is an excellent agreement between the theoretical model predictions and finite-element simulation results for all values of E EOF . Furthermore, far below or above the threshold value of E EOF , the theoretical model describes accurately the transient dynamics and provides a very good prediction for steady-state deformation. In the vicinity of E EOF,CR , due to the small difference in the threshold value of E EOF , we observe that our theoretical model slightly underpredicts the time required to an elastic sheet to collapse onto the floor.
We also performed finite-element numerical simulations for the case of a constant applied current, showing similar agreement with the theoretical model, as presented in Fig. 12 of appendix D.

VII. CONCLUDING REMARKS
In this work, we examined the interfacial instability of a thin film confined between a rigid surface and a pre-stretched elastic sheet, triggered by the pressure formed due to EOF. Applying the lubrication approximation to the flow field and modeling the elasticity by the Euler−Bernoulli beam approximation, we derived a nonlinear viscous−elastic governing equation describing the deformation of an elastic sheet, for constant current and constant voltage actuation modes. Our theoretical analysis revealed that the instability is controlled by a non-dimensional elasto-electro-osmotic number, representing the ratio of electro-osmotic to elastic forces. Through dynamic numerical simulations of the governing equation, we illustrated several distinct modes of instability depending on the electro-osmotic pattern. Furthermore, we demonstrated that this instability can result in an asymmetric deformation pattern, even for symmetric actuation. Finally, we performed finite-element simulations to validate the theoretical model predictions, showing very good agreement.
The study of a lubricated elastic sheet can be viewed in analogy to the classical studies of free-surface thin films, yet some key differences must be highlighted. For example, while instability in thin-film problems ultimately results in rupture of the film, in the case of an elastic sheet, contact is reached with the surface yet the interface remains continuous. An underlying assumption of our modeling is that the membrane remains wetted even when in contact with the surface, similar to the pre-wetting film thickness introduced in the work of Lister et al. [9], yet the question of proper modeling of this contact region (e.g. the effect of van der Waals forces) remains open. Furthermore, the surface-membrane-liquid contact lines obtained in the case of an elastic instability (i.e. a solid-solid-liquid contact, in contrast to the solid-liquid-fluid contact in free surfaces) warrants additional investigation. Lastly, the use of an elastic sheet instead of a free surface opens up new degrees of freedom, such as spatial variation of the membrane thickness or elasticity, which when coupled with the instability mechanism may result in new and interesting dynamics.
While throughout this work we neglected internal tension and considered a strongly pre-stretched elastic sheet, the model can also be extended to a non-pre-stretched elastic sheet. In such a case, based on preliminary simulations, we expect that due to nonlinear coupling between the internal tension and the deflection, the instability will occur at much lower values of E EOF (based on the internal tension), yet at much larger deformation as compared to the pre-stretched case.
Manipulation of fluids using EOF is currently widely encountered in microfluidic devices, that are often fabricated from soft materials such as poly(dimethylsiloxane) (PDMS). The mechanism illustrated in this work may pave the way for implementation of instability-based soft actuators for lab-on-a-chip and soft-robotic applications. EOF is also extensively used as a driving mechanism in nanochannels, where even relatively rigid walls (e.g. glass covers) may results in deformations that are significant relative to the height of the channel. The presented results lay the theoretical foundation for control of the EOF-driven instability in such devices, providing the key features required to either induce or prevent the instability. nonlinear ordinary differential equations as well as eigenvalue problems with various boundary conditions. We obtained the numerical solutions for d ss (x) by beginning with very small values of E EOF and using the asymptotic solution in the limit of small elasto-electro-osmotic number as an initial guess. Solutions for other values of E EOF were then computed through numerical continuation. A second numerical method was employed to explore the dynamic behavior of the instability by solving the nonlinear evolution equation (13). To solve numerically the governing equation (13), we first discretized spatial derivatives in (13) using a second-order central difference approximation with uniform grid spacing, leading to a series of ordinary differential equations for the evolution of d i (t) = d(x i , t). We then integrated forward in time the resulting set of ordinary differential equations using MATLAB's routine ode15s.
We cross-validated our boundary value (first method) and time-dependent (second method) numerical solvers in Figs. 4(a) and 4(b), showing very good agreement.
Appendix C: Variation of growth rate with EEOF and k Figure 11(a) presents the growth rate σ as a function of E EOF for three values of the wavenumber, k, showing that for k = 1 the growth rate monotonically increases as |E EOF | increases and approaches zero only for a negative value of E EOF with E EOF,CR = −6.51. For k = 2, the growth rate varies symmetrically with regards to E EOF = 0, and achieves zero both for positive and negative values of E EOF , E EOF,CR = ±31.3. For k = 3, the growth rate also approaches zero both for positive and negative values of E EOF , though having non-symmetric dependence on E EOF , similarly to the results shown in Fig. 8. Figure 11(b) presents the growth rate σ as a function of k 2 for three values of E EOF , showing that, as expected, for all values of E EOF , σ approaches a constant value of −π 4 as k increases. This behavior can be explained as follows: in the limit of k 1 (kπ 1) and E EOF = O(1), we expect the deformation to be small and thus can expand the deformation (22) terms of a small parameter E EOF /kπ 1, as Substituting (C1) into (23)−(24) yields two uncoupled linear equations for the steady-state deformation d ss (x), and for the eigenvalue problem, subjected to the boundary conditions (18). The values of σ are eigenvalues of (C3) and are given as σ = −B(nπ) 6 − (nπ) 4 − G(nπ) 2 , where n = 1, 2, 3..., indicating that for the case of E EOF /kπ 1 the perturbations always decay and the deformation is stable. For B = G = 0, σ simplifies to σ = −(nπ) 4 , with the maximum growth rate −π 4 that is independent of E EOF and k when k 1, consistent with the results shown in the Fig. 11(b). We performed two-dimensional finite-element numerical simulations with the commercial software COMSOL Multiphysics (version 5.0, COMSOL AB, Stockholm, Sweden) by coupling the Fluid−Structure Interaction module to Electrostatics or Electric Currents modules.
In the Fluid−Structure Interaction module, we fully coupled the unsteady Stokes equations with gravitational body force for the flow to the unsteady Navier equations for the elastic deformation. Additionally, assuming constant permittivity and conductivity, we solved the Laplace equation for the electric potentialṼ in the time-varying domain  using the Electrostatics module (for constant voltage) or Electric Currents module (for constant current). For the case of a constant voltage, we applied a Dirichlet boundary condition,Ṽ = V =Ẽ 0lm , atx = 0 and for the case of a constant current, we prescribed a constant current densityj and applied a Neumann boundary condition, j/σ =Ẽ·x = −∂Ṽ /∂x =Ẽ 0 , atx = 0. The boundary conditions used in the finite-element simulations are summarized in table 1. We discretized the domain using a rectangular mesh with 150 uniformly distributed elements in the longitudinal dimension and 33 uniformly distributed elements in the transverse dimension, three of which reside inside the elastic sheet. We employed the third-order (cubic) discretization for the flow field, the solid's displacement field, and the electric potential, as well as the second-order (quadratic) discretization for the pressure field, resulting in 240 704 degrees of freedom. Additionally, we performed tests to assess the grid sensitivity at this resolution and established grid independence. Finally, in all dynamic finite-element simulations, the solver was forced to take at least a single  time step every 0.1 s, and we stopped the unstable simulations when the film thickness reduced to h = 1 µm. Figure 12 presents the comparison of finite-element simulation results and theoretical model predictions in the case of a constant applied current, showing very good agreement. Tables 2 and 3 summarize the typical values of physical parameters and corresponding non-dimensional numbers used in our numerical finite-element simulations, indicating that the assumptions of the theoretical model, i.e.