Development of soft electromechanical materials is critical for several tantalizing applications such as human-like robots, stretchable electronics, actuators, energy harvesting, among others. Soft dielectrics can be easily deformed by an electric field through the so-called electrostatic Maxwell stress. The highly nonlinear coupling between the mechanical and electrical effects in soft dielectrics gives rise to a rich variety of instability and bifurcation behavior. Depending upon the context, instabilities can either be detrimental, or more intriguingly, exploited for enhanced multifunctional behavior. In this work, we revisit the instability and bifurcation behavior of a finite block made of a soft dielectric material that is simultaneously subjected to both mechanical and electrical stimuli. An excellent literature already exists that has addressed the same topic. However, barring a few exceptions, most works have focused on the consideration of homogeneous deformation and accordingly, relatively fewer insights are at hand regarding the compressive stress state. In our work, we allow for fairly general and inhomogeneous deformation modes and, in the case of a neo-Hookean material, present closed-form solutions to the instability and bifurcation behavior of soft dielectrics. Our results, in the asymptotic limit of large aspect ratio, agree well with Euler's prediction for the buckling of a slender block and, furthermore, in the limit of zero aspect ratio are the same as Biot's critical strain of surface instability of a compressed homogeneous half-space of a neo-Hookean material. A key physical insight that emerges from our analysis is that soft dielectrics can be used as actuators within an expanded range of electric field than hitherto believed.

## Introduction

Soft materials, such as polymers and many soft biological materials, play an important role in our daily life. They can be easily deformed to large strain values due to intrinsically low elastic stiffness. Meanwhile, surface instabilities like wrinkles [1,2] and creases [35] are often observed under mechanical compression or constrained swelling. Soft dielectrics, an important subclass of soft materials, can achieve significantly large deformation when they are subject to electrical stimuli. Soft dielectrics find applications in human-like robots [6,7], stretchable electronics [8], actuators [911], energy harvesters [1214], among others. Large deformations of soft dielectrics are often accompanied by electromechanical instabilities including pull-in instability [1517], wrinkling and the creasing [18], the electro-creasing to cratering instability [19], electro-cavitation [20], among others.

Historically, instabilities are often thought to cause “failure” and usually avoided. The pull-in instability, for example, is suppressed [2127] in order to enhance the actuation strain and the electrical energy density of soft dielectrics. More recently, research has increasingly also been directed at how electromechanical instabilities of soft dielectrics can be harnessed for various applications such as giant actuation strain, dynamic surface patterning, and energy harvesting [28,29].

A commonly used actuator, for example, is a film of dielectric elastomer coated with compliant electrodes on its surfaces. Upon application of a voltage difference between the two electrodes, the Maxwell stress from the electric field compresses the film in the thickness direction, causes expansion in the plane, and creates a large actuation strain. The thinning of the film increases the intensity of the electric field in the material. When the film thickness decreases to a certain threshold value, the film is unable to sustain the Maxwell stress and the pull-in instability occurs. Exploitation of soft dielectric films in applications requires a thorough understanding of large deformation mechanics and the electromechanical instabilities induced by voltages and mechanical forces. To this end, numerous theoretical analyses [10,16,22,3033] have been carried on this subject matter.

In a prior work [22], Zhao and Suo analyzed the electromechanical stability of a film of dielectric elastomer subject to tensile forces in its plane and a voltage difference across its thickness. From the principle of minimum energy, they studied the stability of the homogeneously2 deformed film by examining the positive definiteness of the Hessian matrix. They showed that prestress can significantly enhance the stability of the homogeneously deformed film and markedly increase the actuation stretch. We remark that Zhao and Suo [22] assumed a homogeneously deformed film throughout their equilibrium state and stability analysis. Subsequently, this assumption of homogeneous deformation has been widely used in other works [10,24,25,30,31,34,35].

The aforementioned assumption of a homogeneous deformation imposes the restriction that the upper and bottom surfaces of the dielectric thin film remain perfectly plane. Hence, nonhomogeneous deformation and the effects of the geometry of the dielectric film, like the thickness or the aspect ratio, on the electromechanical instability are excluded. In a recent work, Dorfmann and Ogden [33] investigated the instability (buckling) of an infinite plate of electroelastic material by analyzing its incremental elastic deformation. In another work, Dorfmann and Ogden [36] studied the surface instability of a half-space subject to both mechanical compression and an electric field normal to its surface.

In this work, we present a complete linearized bifurcation analysis for electromechanical instability in a finite block of a soft dielectric material subject to physically reasonable boundary conditions,3 and compare it with the response of a thin film and half-space. An elastic finite block is often used to study the mechanical behavior of elastic materials at finite strain, such as the instability [37] and post-buckling [38] of a mechanical compressed elastic block, and the buckling of a compressible magnetoelastic block [39]. Unlike a half-space [1,36] or an infinite long plate [33], a finite block has measurable length quantities, such as the aspect ratio, and allows for physically well-defined boundary conditions on all its surfaces. Hence, the effect of the boundary conditions due to finite dimensions on the electromechanical instability can be addressed in the present work by analyzing a finite block. Compared to the in-plane biaxial dead loads on the dielectric film, as used in past works [22], we employ displacement-controlled boundary conditions on the two sides of the finite block which allows a facile consideration of both tension and compression.4 Based on the implicit function theorem [40,41], we present an analysis of the onset of bifurcation from the trivial solution of a finite block of dielectric elastomer subject to mechanical loads (compression or extension) and a voltage across its thickness. Although our analysis of electromechanical instability is applicable to a general elastic dielectric elastomer, we present closed-form expressions for the special case of ideal neo-Hookean dielectrics.

The paper is organized as follows. In Sec. 2, we present the general formulation for the electrostatic problem of a finite block of a dielectric elastomer subject to electromechanical loads. The linear bifurcation analysis is presented in Sec. 3, where the incremental boundary-value problem is obtained by linearizing the equations of equilibrium with respect to deformation and the polarization. In Sec. 4, we obtain the solutions of the homogeneous deformation and the incremental boundary-value problem, and discuss the onset of bifurcation from the trivial solution. Finally, in Sec. 5, we compare our analytical results with Euler's predictions for the buckling of both mechanically and electromechanically compressed slender block and discuss the pertinent physical insights.

## Formulation

### Domain and Boundary Conditions.

Consider a finite block of an elastic dielectric (see Fig. 1). Assuming plane-strain condition in the X3 direction, the dielectric block in the reference configuration can be represented by

Fig. 1
Fig. 1
Close modal
$ΩR={X∈ℝ2:0≤X1≤l1,−l22≤X2≤l22}$
(1)
where X1 and X2 are the Cartesian coordinates, l1 is the length, and l2 is the height of the dielectric block. The boundary $∂ΩR$ of ΩR consists of four parts
$Sl={X∈ΩR:X1=0}, Sr={X∈ΩR:X1=l1}Su={X∈ΩR:X2=l22}, Sb={X∈ΩR:X2=−l22}$
(2)
The deformation of the block is expressed by a smooth function $x:ΩR→ℝ2$, and the constraint of incompressibility requires a unit Jacobian, such that
$J=detF=1$
(3)

where $F=∇x=(∂x/∂X1)e1+(∂x/∂X2)e2$ is the deformation gradient in two dimensions, and $ei$, i = 1, 2, are the unit vectors in the Xi directions.

A few comments regarding the boundary conditions are in order. For a dielectric elastomer film, a voltage is usually applied across the top and bottom surfaces and in-plane tensile dead loads are introduced. Such kinds of boundary conditions [22] are used for dielectric elastomers that work in uniaxial actuation mode, such as in spring-roll actuators [10]. In this paper, we apply a potential difference between the upper and lower surfaces, but mimic the physical situation where the dielectric block is controlled by a loading device that stretches or compresses the block in the X1 direction. For example, the loading device can be made of two well-lubricated, rigid plates (see Fig. 1) that are in contact with the side surfaces $Sl∪Sr$ with rollers. During the compression/tension process, there is no rotation of the two plates in order to ensure that the compressive/tensile stresses on the rigid plates are always along the X1 direction. The mechanical and electrostatic boundary conditions on the side surfaces $Sl∪Sr$ are defined by
$x·e1=λX1, D̃·e1=0 on Sl∪Sr$
(4)
where $λ>0$ is the prescribed stretch along the X1 direction and $D̃$ is the nominal electric displacement. Meanwhile, the mechanical and electrostatic boundary conditions on the upper and lower surfaces $Su∪Sb$ are
$t̃e=0, ξ=ξb on Su∪Sb$
(5)

where $t̃e$ is the surface traction and ξ is the voltage. Here, the prescribed voltages are $ξb=V$ on the upper surface $Su$ and $ξb=0$ on the lower surface $Sb$ (see Fig. 1).

### Equations of Electrostatics of a Deformable Media.

In the deformed domain Ω, the true electric field is denoted by e, the true electric displacement by d, and the polarization by p. In the absence of free charges, currents, and magnetic fields, the Maxwell equations reduce to
$curl e=0, div d=0, d=ϵ0e+p in Ω$
(6)

where ϵ0 is the vacuum permittivity. The curl, the divergence, and the gradient operators in the current configuration are denoted by “curl,” “div,” and “grad,” respectively. In contrast, the corresponding operators in the reference configuration are denoted by “Curl,” “Div,” and “$∇$.” The equality, $curl e=0$ in Eq. (6), indicates that there exists a scalar potential (voltage) ξ such that $e=−grad ξ$.

To represent e, d, and p in the undeformed domain ΩR, we use the composition maps [42]
$E=e°x, D=d°x, P=p°x$
(7)
In the undeformed domain ΩR, the nominal electric field is denoted by $Ẽ$, the nominal electric displacement by $D̃$, and the nominal polarization by $P̃$. We define the relations of the true fields in Eq. (7) and the nominal fields in ΩR as [42,43]
$Ẽ=FTE, D̃=JF−1D, P̃=JP$
(8)
where $FT$ and $F−1$ are the transpose and the inverse of the deformation gradient F, respectively, and $J=detF$ is the Jacobian. Maxwell's equations (Eq. (6)) in ΩR are
$Curl Ẽ=0, Div D̃=0, D̃=F−1(ϵ0JF−TẼ+P̃)$
(9)

where $Ẽ=−∇ξ$.

Combining Eq. (8), an alternative form of Eq. (9) in ΩR can be written as
$FTE=−∇ξ, Div D̃=0, FD̃=ϵ0JE+P̃$
(10)

### Free Energy of the System.

Electromechanics of deformable dielectrics can be formulated in a variety of ways cf. Refs. [36,42,43] for just a few examples. In this paper, we follow the energy formulation of continuum magneto-electro-elasticity as described by Liu [42]. The notion of invoking a minimum energy principle with Maxwell's equations as a constraint has roots in an earlier work on micromagnetism [44] and ferroelectrics [45].

Subject to both mechanical and electrical loads, the total free energy (see, for example, Refs. [42,46]) of the system in Fig. 1 is given by
$F[x,P̃]=U[x,P̃]+Eelect[x,P̃]$
(11)
Here, $U[x,P̃]$ is the internal energy
$U[x,P̃]=∫ΩRΨ(F,P̃)$
(12)
where $F=∇x$ and $Ψ(F,P̃)$ is the internal energy density. The electric energy, $Eelect[x,P̃]$, in Eq. (11) is
$Eelect[x,P̃]=ϵ02∫ΩRJ|E|2+∫Su∪SbξD̃·N$
(13)

where $J=det∇x$ is the Jacobian and N is the unit normal to the surfaces $Su∪Sb$. The relations among E, $D̃$, and $P̃$ in Eq. (13) are given by Eq. (10). Note that the mechanical work done by the loading device on the side surfaces $Sl∪Sr$ is not included into the total free energy due to the nominal displacement-controlled boundary condition (Eq. (4)1).

### First Variation of the Free Energy.

When the aforementioned electromechanical system is in equilibrium at a deformation x and a polarization $P̃$, the first variation of the energy functional $F[x,P̃]$ must vanish (subject to the constraint imposed by Maxwell's equations). Since there exist two functions $x:ΩR→Ω$ and $P̃:ΩR→ℝ2$ in $F[x,P̃]$, the vanishing of the first variation requires that both the first variations with respect to x and $P̃$ must be zero (see Appendix  A for details).

#### Variation of Polarization.

The first variation of $F[x,P̃]$ in Eq. (11) with respect to the polarization $P̃$ gives
$∂Ψ∂P̃−E=0 in ΩR$
(14)

where $E=−F−T∇ξ$. The detailed derivation of Eq. (14) is given in Appendix A1.

#### Variation of Deformation.

Vanishing of the first variation of $F[x,P̃]$ in Eq. (11) with respect to the deformation x yields the Euler–Lagrange equation (see Appendix A2 for details)
$Div(∂Ψ∂F+Σ̃−qF−T)=0 in ΩR$
(15)
and the natural boundary conditions
$(∂Ψ∂F+Σ̃−qF−T)e1=se1 on Sl∪Sr$
(16)
$(∂Ψ∂F+Σ̃−qF−T)e2=0 on Su∪Sb$
(17)
where q is the hydrostatic pressure required by the incompressibility constraint (Eq. (3)), s is the normal stress on the side surfaces $Sl∪Sr$, and $Σ̃$ is the so-called Piola–Maxwell stress defined by
$Σ̃=E⊗D̃−ϵ0J2|E|2F−T$
(18)
Equations (10), (14)(17), along with the constraint of incompressibility (Eq. (3)) and the boundary conditions (Eqs. (4) and (5)), form a boundary-value problem, whose solution includes all the possible equilibrium solutions for a finite block of soft dielectric subject to mechanical and electrical loads. The aforementioned boundary-value problem can be compactly summarized as:
$DivT=0, ∂Ψ∂P̃−E=0, detF=1FTE=−∇ξ, DivD̃=0 FD̃=ϵ0JE+P̃} in ΩR$
(19)
$x1=λX1, D̃·e1=0, Te1=se1 on Sl∪Sr$
(20)
$ξ=ξb, Te2=0 on Su∪Sb$
(21)
where the total nominal stress T is
$T=∂Ψ∂F+Σ̃−qF−T$
(22)

## Onset of Electromechanical Buckling

Of interest here is the condition for the onset of buckling of the dielectric block. Mathematically, buckling is governed by the onset of bifurcation in the trivial solution to the boundary-value problem (Eqs. (19)(21)). Based on the implicit function theorem [40,41], the equilibrium equations have a nontrivial solution bifurcating from its trivial solution only if the linearized equations of equilibrium possess a nonzero solution. It is obvious that the onset of bifurcation depends on the applied mechanical and electrical loads. The linearized equations describe the response of the dielectric block, in a state of equilibrium, to infinitesimal increments of the deformation and the polarization.

### Linearization With Respect to Both the Deformation and the Polarization.

Let $x*$ and $P̃*$ be the infinitesimal increments of the deformation x and the polarization $P̃$, respectively. Other increments depend on $x*$ and $P̃*$ at $(x,P̃)$. We denote other linearized increments (omitting higher terms) by taking advantage of the superscripts. For example, the total linearized increment of a general field Θ is denoted by $Θ*$ (see Appendix  B for details). $Θ*$ is usually the sum of two increments $Θ†$ and $Θ‡$, which denote, respectively, the increments related to the deformation and the polarization. Since the deformation gradient F and the Jacobian J are independent of the polarization, we have the linearized increments
$F*=∇x*, J*=JF−T·∇x*,(FT)*=(∇x*)T, (F−T)*=−F−T(∇x*)TF−T$
(23)
In contrast, the total linearized increments of the electric field E, the nominal electric displacement $D̃$, and the Piola–Maxwell stress $Σ̃$ consist of two parts:
$E*=E†+E‡, D̃*=D̃†+D̃‡, Σ̃*=Σ̃†+Σ̃‡$
(24)
Combining Eq. (24) and $Σ̃$ in Eq. (18), we have the relation between the linearized increments (see Appendix B2)
$Σ̃*=E*⊗D̃+E⊗D̃* −ϵ02{2J(E·E*)F−T+|E|2[J(F−T)*+J*F−T]}$
(25)
Similarly, the linearization of the boundary-value problem (Eqs. (19)(21)) can be written as
$DivT*=0, ∂2Ψ∂P̃∂F[F*]+∂2Ψ∂P̃2P̃*−E*=0F−T·F*=0, (FT)*E+FTE*=−∇ξ*DivD̃*=0, FD̃*+F*D̃=ϵ0JE*+ϵ0J*E+P̃*}in ΩR$
(26)
$x*·e1=0, D̃*·e1=0, T*e1=s*e1 on Sl∪Sr$
(27)
$ξ*=0, T*e2=0 on Su∪Sb$
(28)
where $T*$ is the linearized increment of the total nominal stress T in Eq. (22), given by
$T*=∂2Ψ∂F2[F*]+∂2Ψ∂F∂P̃P̃*+Σ̃*−q*F−T+qF−T(FT)*F−T$
(29)

We remark that the linearized boundary-value problem (Eqs. (26)(28)) considers the total increments including both the incremental deformation and the incremental polarization. The condition of the nonzero solution of $(x*,P̃*)$ in Eqs. (26)(28) determines the onset of the electrical buckling—and more precisely the onset of bifurcation from the solution $(x,P̃)$ —of a finite block of dielectric elastomer subject to electromechanical loads.

### Linearization With Respect to Only the Deformation.

A further simplification of Eqs. (26)(28) may be made by considering the linearization with respect to only the deformation. That is, we introduce an infinitesimal increment $x*$ of the deformation but a zero increment of the polarization when we linearize the boundary-value problem (Eqs. (19)(21)) at $(x,P̃)$. This simplification, of course, will yield a narrower solution space that is a subspace of the solution space considering both the incremental deformation and polarization; however, it significantly simplifies the analysis and also provides important results for electromechanical buckling. Therefore, by letting $P̃*=0$ the total increments (with superscript “*”) in Eqs. (24) and (25) reduce to the increments (with superscript “”) with respect to the deformation, the linearized boundary-value problem (Eqs. (26)(28)) can be reduced to
$DivT†=0, ∂2Ψ∂P̃∂F[F*]−E†=0F−T·F*=0, (FT)*E+FTE†=−∇ξ†DivD̃†=0, FD̃†+F*D̃=ϵ0JE†+ϵ0J*E} in ΩR$
(30)
$x*·e1=0, D̃†·e1=0, T†e1=s†e1 on Sl∪Sr$
(31)
$ξ†|SbSu=0, T†e2=0 on Su∪Sb$
(32)
where
$T†=∂2Ψ∂F2[F*]+Σ̃†−q*F−T+qF−T(FT)*F−T$
(33)
and
$Σ̃†=E†⊗D̃+E⊗D̃†−ϵ02{2J(E·E†)F−T+|E|2[J(F−T)*+J*F−T]}$
(34)

We remark here that the total incremental boundary-value problem (Eqs. (26)(28)) and the reduced incremental boundary-value problem (Eqs. (30)(32)) are valid for all incompressible elastic soft dielectrics. In the following, we will adopt the neo-Hookean constitutive law to generate specific results.

## Neo-Hookean Dielectrics

In the following, we consider incompressible neo-Hookean dielectrics, whose strain energy function [42,46] under the plane strain assumption is given by
$Ψ(F,P̃)=μ2(|F|2−2)+|P̃|22J(ϵ−ϵ0)$
(35)

where μ is the shear modulus, and ϵ and ϵ0 are, respectively, the permittivities of the dielectric elastomer and the vacuum. Note that the second term on the right-hand side of Eq. (35) reflects the usual linear dielectric behavior, that is, the permittivity ϵ of the dielectric elastomer is independent of the deformation.

The derivatives of the strain energy function (Eq. (35)) are given by
$∂Ψ∂F=μF−|P̃|22J(ϵ−ϵ0)F−T, ∂Ψ∂P̃=P̃J(ϵ−ϵ0)∂2Ψ∂F2=μI4+|P̃|22J(ϵ−ϵ0)(F−T⊗F−T−∂F−T∂F)∂2Ψ∂P̃2=I2J(ϵ−ϵ0),∂2Ψ∂F∂P̃=−F−T⊗P̃J(ϵ−ϵ0), ∂2Ψ∂P̃∂F=−P̃⊗F−TJ(ϵ−ϵ0)}$
(36)

where $I4$ and $I2$ are, respectively, the fourth- and second-order identity tensors in two dimensions.

### Homogeneous Deformation.

Substituting Eq. (36) into the boundary-value problem (Eqs. (19)(21)), a trivial solution that corresponds to homogeneous deformation is given by
$x0(X)=λX1e1+λ−1X2e2, P̃0(X)=−(ϵ−ϵ0)λẼ0e2$
(37)
where $Ẽ0=V/l2$ and other corresponding quantities are
$F0(X)=∇x0(X)=λe1⊗e1+λ−1e2⊗e2E0(X)=−λẼ0e2, D̃0(X)=−ϵλ2Ẽ0e2}$
(38)
The Piola–Maxwell stress tensor (Eq. (18)) is
$Σ̃0:=[−ϵ02λẼ0200(ϵ−ϵ02)λ3Ẽ02]$
(39)
the hydrostatic pressure (i.e., the Lagrange multiplier) is
$q0=μλ−2+ϵ2λ2Ẽ02$
(40)
the total stress tensor (Eq. (22)) is
$T0:=[μ(λ−λ−3)−ϵλẼ02000]$
(41)
and the compressive/tensile stress on $Sl∪Sr$ is
$s0=μ(λ−λ−3)−ϵλẼ02$
(42)

Note that a negative/positive s0 in Eq. (42) represents the nominal compressive/tensile stress on the side surfaces $Sl∪Sr$ in the reference configuration. In contrast, the true compressive/tensile stress in the current configuration is given by $λs0$.

Figure 2 shows how the electric field affects the mechanical behavior of the homogeneous deformation of a finite block. The dimensionless compressive/tensive stress $s0/μ$ and electric field $Ẽ0ϵ/μ$ are used. In the absence of the electric field (i.e., $Ẽ0ϵ/μ=0$) in Fig. 2(a) (or in Eq. (42)), for example, a prescribed stretch $λ>1$ corresponds to a tensile stress $s0>0$, while a stretch $λ<1$ corresponds to a compressive stress $s0<0$ on the side surfaces $Sl∪Sr$ of the block.

Fig. 2
Fig. 2
Close modal

The electric field in Eq. (42) will decrease the nominal stress s0 on the side surfaces. This is because the Maxwell stress in Eq. (39) will make the block decrease its height l2 (due to a positive component of the Maxwell stress in the X2 direction) and increase its length l1 (due to a negative component of the Maxwell stress in the X1 direction). However, the two lubricated rigid plates exert an additional compressive stress on the side surfaces to hinder the extension of the block. Therefore, at a prescribed stretch λ in Fig. 2(a), the electric field $Ẽ0ϵ/μ$ can decrease the nominal stress vector $s0/μ$ (or the true stress vector $λs0/μ$). For example, at a prescribed stretch $λ>1$, the increase of the electric field $Ẽ0ϵ/μ$ can cause the nominal stress $s0/μ$ in Fig. 2(a) (or the true stress $λs0/μ$ in Fig. 2(b)) decrease from positive (tensile stress) to negative (compressive stress). If we were to ignore electrical breakdown, the continually increasing compressive stress will eventually force the block to buckle.

To further illustrate the effects of the electric field on deformation, we consider two special cases: a prescribed stretch λ = 1 and a zero nominal stress $s0=0$.

In the first case, the block is undeformed prior to electromechanical buckling. This is because of the constraint of incompressibility and the plane strain assumption in our model, leading to the stretch ratios $λ1=λ=1, λ2=λ−1=1$, and $λ3=1$. Although the block is undeformed under the electric field, it is no longer a stress-free state. The nominal compressive stress s0 in Eq. (42) is $s0=−ϵẼ02$, which is a quadratic function of the nominal electric field (see Fig. 3(a)). Under zero electric field, the block is stress-free, corresponding to the origin $(Ẽ0ϵ/μ,s0/μ)=(0,0)$. It is clear from Fig. 3(a) that the parabola opens downward and the axis of symmetry is $Ẽ0ϵ/μ=0$. In this case, the electric field always induces a compressive state in the block. With a continuously increasing electric field, the block eventually will buckle. Note that only electromechanical buckling is considered in this paper for the constrained deformation—other instabilities including the electrical breakdown, the electro-creasing to cratering instabilities, and the electrocavitation instability [14,1820] are beyond the scope of this paper.

Fig. 3
Fig. 3
Close modal

In the second case, the homogeneously deformed block is stress-free. With $s0=0$ in Eq. (42), the relation between the stretch and the electric field becomes $λ=(1−ϵẼ02/μ)−1/4$ (see Fig. 3(b)). Without considering the electrical breakdown and the pull-in instability [22], the stretch λ, mathematically, can increase from one to infinity as the dimensionless electric field $Ẽ0ϵ/μ$ increases from zero to one.

In the following, we will analyze electromechanical buckling by studying the solution of the incremental boundary-value problem (Eqs. (30)(32)) at the homogeneous solution (Eqs. (37)(42)).

### Incremental Boundary-Value Problem.

Let us first address Eq. (30) in ΩR. The increments of the Jacobian $J*$ and the electric field $E†$ in Eq. (30) must vanish due to the constraint of incompressibility, $F−T·∇x*=0$, namely
$J*=JF−T·∇x*=0$
(43)
$E†=−P̃⊗F−TJ(ϵ−ϵ0)[∇x*]=−(F−T·∇x*)P̃J(ϵ−ϵ0)=0$
(44)
and the two relations, $(FT)*E+FTE†=−∇ξ†$ and $FD̃†+F*D̃=ϵ0JE†+ϵ0J*E$ in Eq. (30), reduce to
$∇ξ†=−(∇x*)TE, D̃†=−F−1F*D̃$
(45)
Then, the incremental Piola-Maxwell stress (Eq. (34)) reduces to
$Σ̃†=E⊗D̃†−ϵ0J|E|22(F−T)*$
(46)
Note that Eqs. (43)(46) hold for all kinematically admissible deformations of incompressible neo-Hookean dielectrics, including the homogeneous deformation. Moreover, the constraint makes the divergence of $D̃†$ in Eq. (45) vanish automatically in the case of homogeneous deformation, such that
$DivD̃†=−Div(F0−1F*D̃0)=−D̃0·∇(F0−T·F*)=0$
(47)
Based on Eqs. (43)(47) and the homogeneous solution (Eqs. (37)(42)), the linearized boundary-value problem (Eqs. (30)(32)) reduces to
$DivT†=0, F0−T·∇x*=0(∇x*)TE0=−∇ξ†, D̃†=−F0−1(∇x*)D̃0} in ΩR$
(48)
$x*·e1=0, D̃†·e1=0, T†e1=s†e1 on Sl∪Sr$
(49)
$ξ†|SbSu=0, T†e2=0 on Su∪Sb$
(50)
where the incremental voltage $ξ†$ and the incremental nominal electric displacement $D̃†$ are
$ξ†=λẼ0x2*+ξ0, D̃†=ϵλẼ0(x1,2*e1+λ2x2,2*e2)$
(51)
and the incremental nominal stress $T†$ is
$T†:=μ[T11†T12†T21†T22†]$
(52)
with
$T11†=[1+λ−4+(Ẽ0ϵ/μ)2]x1,1*−λ−1μ−1q*T12†=x1,2*+[λ−2+λ2(Ẽ0ϵ/μ)2]x2,1*T21†=x2,1*+λ−2x1,2*, T22†=2x2,2*−λμ−1q*$

#### Governing Equation.

From the constraint of incompressibility, $F0−T·∇x*=0$, in Eq. (48), we introduce a stream function $ϕ(X1,X2)$
$x1*=λϕ,2(X1,X2), x2*=−λ−1ϕ,1(X1,X2)$
(53)
where the subscript denotes the partial derivative. With Eqs. (52) and (53), $DivT†=0$ in Eq. (48) gives
$λ2(ϕ,112+ϕ,222)−μ−1q,1*=0$
(54a)
$λ−2(ϕ,111+ϕ,122)+μ−1q,2*=0$
(54b)
Eliminating $q*$ in Eq. (54) yields
$ϕ,1111+(λ4+1)ϕ,1122+λ4ϕ,2222=0$
(55)

#### Boundary Conditions on $Sl∪Sr$ and on $Su∪Sb$⁠.

Substituting Eq. (53) into the boundary conditions on $Sl∪Sr$ in Eq. (49), we have
$ϕ,2=0, ϕ,22=0, ϕ,22−ϕ,11=0 on Sl∪Sr$
(56)

where the first and the third equations come from the mechanical boundary conditions (i.e., the controlled-nominal displacement $x*·e1=0$ and the free shear stress $T†e1=s†e1$), while the second equation corresponds to the electrostatic boundary condition (i.e., $D̃†·e1=0$ in Eq. (49)).

Similarly, substituting Eq. (53) into the electrostatic boundary condition on $Su∪Sb$ in Eq. (50), we have
$ϕ,1(X1,l22)=ϕ,1(X1,−l22)$
(57)
Moreover, substituting Eq. (53) into the mechanical boundary conditions on $Su∪Sb$ in Eq. (50) yields
$[λ−4+(Ẽ0ϵ/μ)2]ϕ,11−ϕ,22=0 2μϕ,12+λ2q*=0} on Su∪Sb$
(58)

#### Solution of the Incremental Boundary-Value Problem.

Due to the boundary conditions on $Sl∪Sr$ in Eq. (56), the solution of $ϕ(X1,X2)$ in Eq. (55) admits the series form
$ϕ(X1,X2)=∑m=1∞Ym(X2)sin(kmX1)+BX1+C,$
(59)
where $km=(mπ/l1), m=1,2,3,…$, B and C are constants. C is physically irrelevant and may be chosen to be zero. B is related to rigid body motion of the deformed dielectric block. Substituting Eq. (59) into Eq. (53), we have
$x1*=λ∑m=1∞Ym′(X2)sin(kmX1)$
(60a)
$x2*=−λ−1{∑m=1∞kmYm(X2)cos(kmX1)+B}$
(60b)
where the prime denotes the derivative with respect to X2. Substituting Eq. (59) into Eq. (54) and performing the integration with respect to X1, we obtain
$q*=−μλ2∑m=1∞km−1[Ym‴(X2)−km2Ym′(X2)]cos(kmX1)$
(61)
Then the increment $s†$ of the nominal stress in Eq. (49) is
$s†=μλ{∑m=1∞[km−2Ym‴(X2) +(λ−4+(Ẽ0ϵ/μ)2)Ym′(X2)]km cos(kmX1)}$
(62)
Substituting Eq. (59) into the electrostatic boundary conditions on $Su∪Sb$ in Eq. (57) and using the orthogonality relation of Fourier series yields
$Ym(l22)=Ym(−l22)$
(63)
Moreover, substituting Eqs. (53), (59), and (61) into the mechanical boundary conditions on $Su∪Sb$ in Eq. (58) and using again the orthogonality relation of Fourier series, we obtain
$Ym″(X2)+[λ−4+(Ẽ0ϵ/μ)2]km2Ym(X2)=0Ym‴(X2)−(1+2λ−4)km2Ym′(X2)=0} on Su∪Sb$
(64)
Finally, substituting Eq. (59) into the governing equation (55), we find that $Ym(X2)$ yields the following fourth-order ordinary differential equation:
$Ym(4)(X2)−km2(1+λ−4)Ym″(X2)+km4λ−4Ym(X2)=0$
(65)
The general solution of $Ym(X2)$ in Eq. (65) is
$Ym(X2)={C1mcosh(kmλ−2X2)+C2msinh(kmλ−2X2) +C3mcosh(kmX2)+C4msinh(kmX2) for λ≠1(C¯1m+C¯2mX2)cosh(kmX2) +(C¯3m+C¯4mX2)sinh(kmX2) for λ=1$
(66)

where $Cim$ and $C¯im, i=1,2,3,4$, are constant coefficients.

### Bifurcation at Varying Stretch $λ≠1$⁠.

Substituting the general solution to $Ym(X2)$ in Eq. (66)$1$ for $λ≠1$ into Eq. (64), we obtain a system of four linear equations in four unknowns $Cim, i=1,2,3,4$. The system of four equations can be rewritten in a matrix form of $MCm=0$, where M is the 4 × 4 coefficient matrix and $Cm=(C1m,C2m,C3m,C4m)T$. The nonzero solution of $Cm$ requires a zero determinant of M, namely
$detM=|M110M1300M220M240M320M34M410M430|=0$
(67)
which can be decomposed into a product of two 2 × 2 determinants, such that
$|M11M13M41M43|•|M22M24M32M34|=0$
(68)
where
$M11=[2λ−4+(Ẽ0ϵ/μ)2]coshmπl22λ2l1M13=[1+λ−4+(Ẽ0ϵ/μ)2]coshmπl22l1M22=(1+λ−4)coshmπl22λ2l1, M24=2λ−2coshmπl22l1M32=[2λ−4+(Ẽ0ϵ/μ)2]sinhmπl22λ2l1M34=[1+λ−4+(Ẽ0ϵ/μ)2]sinhmπl22l1M41=(1+λ−4)sinhmπl22λ2l1, M43=2λ−2sinhmπl22l1$

Equation (68) holds if either of the 2 × 2 determinant vanishes, indicating the possibility of two types of buckling.

For the first type, the vanishing of the left 2 × 2 determinant in Eq. (68) admits nonzero $C1m$ and $C3m$ but zero $C2m$ and $C4m$, leaving two hyperbolic cosine functions of $X2∈[−l2/2,l2/2]$ in $Ym(X2)$ in Eq. (66)1 and making it become an even function of X2. This type of electrical buckling, of course, satisfies the electrostatic boundary conditions, Eq. (63), since $Ym(X2)$ is an even function. Moreover, the even function, $Ym(X2)$ in Eq. (66)1, makes the perturbed displacement $x1*(X1,X2)$ in Eq. (60) become an odd function of X2 and $x2*(X1,X2)$ become an even function of X2, such as $x1*(X1,X2)=−x1*(X1,−X2)$ and $x2*(X1,X2)=x2*(X1,−X2)$. It is assumed that the constant B in Eq. (60) for the coordinates is appropriately chosen to make $x2*(0,0)=0$. Then, the buckling modes of the first type are antisymmetric with respect to the X1 axis. This type of buckling is called an antisymmetric buckling about the X1 axis. For instance, Figs. 4(a) and 5(a) are antisymmetric bifurcation modes with m = 1 and m = 2, respectively.

Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal

For the second type, the right 2 × 2 determinant in Eq. (68) vanishes and then $Ym(X2)$ in Eq. (66)1 has nonzero $C2m$ and $C4m$ but zero $C1m$ and $C3m$. Thus, $Ym(X2)$ only contains two hyperbolic sine functions of X2 (i.e., $Ym(X2)$ becomes an odd function of X2). The perturbed displacement $x1*(X1,X2)$ in Eq. (60) is an even function of X2, such that $x1*(X1,X2)=x1*(X1,−X2)$, while the perturbed displacement $x2*(X1,X2)$ in Eq. (60) has the property $x2*(X1,X2)+x2*(X1,−X2)=−2λ−1B$. If the constant B in Eq. (60) is chosen to be zero for an appropriate fixity condition of coordinates, the perturbed displacement $x2*(X1,X2)$ in Eq. (60) becomes an odd function of X2, such as $x2*(X1,X2)=−x2*(X1,−X2)$. This type of electrical buckling satisfies the mechanical boundary conditions (Eq. (64)), however, it does not satisfy the electrostatic boundary conditions (Eq. (63)) since $Ym(X2)$ is an odd function now. For the mechanical compression, the buckling modes of the second type are symmetric with respect to the X1 axis. This type of buckling is called a symmetric buckling about the X1 axis. The symmetric bifurcation modes with m = 1 and m = 2 are shown, respectively, in Figs. 4(b) and 5(b).

The critical conditions for the two types of buckling can be explicitly written as

Type (i): Antisymmetric
$(1+λ−4)[1+λ−4+(Ẽ0ϵ/μ)2]tanhmπl22λ2l1 −2λ−2[2λ−4+(Ẽ0ϵ/μ)2]tanhmπl22l1=0$
(69)
Type (ii): Symmetric
$(1+λ−4)[1+λ−4+(Ẽ0ϵ/μ)2]tanhmπl22l1 −2λ−2[2λ−4+(Ẽ0ϵ/μ)2]tanhmπl22λ2l1=0$
(70)

### Bifurcation at Fixed Stretch λ = 1.

Similarly, substituting the general solution to $Ym(X2)$ in Eq. (66)2 for λ = 1 into Eq. (64), we obtain a system of four linear equations in four unknowns $C¯im, i=1,2,3,4$. The system of four equations can be written in a matrix form as $M¯C¯m=0$, where $M¯$ is the 4 × 4 coefficient matrix and $C¯m=(C¯1m,C¯2m,C¯3m,C¯4m)T$. A nonzero solution of $C¯m$ requires $detM¯=0$, which, similar to Eq. (68), can be reduced to a product of two 2 × 2 determinants, such that
$|M¯11M¯14M¯41M¯44|•|M¯22M¯23M¯32M¯33|=0$
(71)
where
$M¯11=[2+(Ẽ0ϵ/μ)2]coshmπl22l1M¯14=2l1mπcoshmπl22l1+l22[2+(Ẽ0ϵ/μ)2]sinhmπl22l1M¯22=l22sinhmπl22l1, M¯23=coshmπl22l1M¯32=2l1mπsinhmπl22l1+l22[2+(Ẽ0ϵ/μ)2]coshmπl22l1M¯33=[2+(Ẽ0ϵ/μ)2]sinhmπl22l1M¯41=sinhmπl22l1, M¯44=l22coshmπl22l1$
Similar to the analysis of Eq. (68), the left 2 × 2 determinant in Eq. (71) might correspond to the type of an antisymmetric buckling, while the right one might corresponds to the type of a symmetric buckling. However, the right 2 × 2 determinant in Eq. (71) is always nonzero, indicating the nonexistence of the type of symmetric buckling at λ = 1 regardless of how large the electric field is. While both symmetric and antisymmetric buckling satisfy all the mechanical boundary conditions for purely mechanical compression, only antisymmetric buckling satisfies both the mechanical and electrostatic boundary conditions for the combined electromechanical loading. The critical condition for the antisymmetric buckling of a dielectric block under electric field at λ = 1 in plane strain comes from the vanishing of the left 2 × 2 determinant in Eq. (71), and yields
$[1+12(Ẽ0ϵ/μ)2]mπl2l1−sinhmπl2l1=0$
(72)

## Discussion and Conclusions

### Comparison With Euler's Prediction for the Mechanical Compression.

The buckling of Euler's column studied by Leonhard Euler in 1757 is one of the classical problems in engineering. The formula derived by Euler gives the critical load at which a long, slender, ideal column is in a state of unstable equilibrium (i.e., even an infinitesimal lateral force will make the column buckle).

Considering the conditions of end support of the column, Euler's formula can be expressed as
$F=π2EeIe(Kl1)2$
(73)

where F is the critical force, Ee is the plane strain elastic Young's modulus, Ie is the area moment of inertia of the cross section, l1 is the length of the column, and K is column effective length factor that depends on the conditions of end support. For example, the factor K is 0.5 for both fixed ends while it is 1 for both pinned ends.

The effective Young's modulus under plane strain is $Ee=4μ$ in Eq. (73) for incompressible neo-Hookean materials with shear modulus μ, and the area moment of inertia is $Ie=l23/12$ for a block with height l2 and unit width. Thus, the critical nominal stress sc from Eq. (73) is
$sc=Fl2×1=μπ23(Kl1/l2)2$
(74)

In contrast to Euler's formula for slender structures, our buckling analysis is valid for a finite compressed elastic block with any aspect ratio $l1/l2$. In the absence of an electric field, Eq. (69) gives the critical stretch λc for the buckling of an incompressible neo-Hookean block subject to a purely mechanical compression. With the relation between the stretch and the nominal stress in Eq. (42), we can obtain the critical nominal stress that corresponds to the critical stretch λc. We remark here that Eq. (69) determines the critical stretch of antisymmetric buckling and Eq. (42) is used to transform the critical stretch into the critical nominal stress for the purpose of a direct comparison with Euler's prediction (Eq. (74)). Note that only the results of the antisymmetric buckling are used to compare with Euler's prediction (Eq. (74)) since the antisymmetric buckling always occurs prior to symmetric buckling. Moreover, only antisymmetric buckling satisfies all the mechanical and electrostatic boundary conditions, while the symmetric buckling satisfies only the mechanical boundary conditions. The detailed discussion of the difference between antisymmetric and symmetric buckling is given in Secs. 5.3 and 5.4.

Figure 6(b) shows the critical nominal stress for the antisymmetric buckling mode m = 2 in Eq. (69) of a finite block, whose buckling pattern is shown schematically in Fig. 5(a). The buckling pattern and the boundary condition on the left and the right surfaces in Fig. 5(a) are very similar to that of the buckling of Euler's column with fixed-fixed ends. Compared with Euler's prediction, the two predicted critical loads of buckling agree well with each other only at a sufficiently large aspect ratio (i.e., $l1/l2>5$). The obvious discrepancy at small aspect ratios is because Euler's analysis is only valid for a slender column.

Fig. 6
Fig. 6
Close modal

### Comparison of Euler's Prediction for Electroelastic Buckling at a Fixed Stretch λ = 1.

Compared with the mechanical compressive stress, the electrostatic Maxwell stress can also make the dielectric block buckle in our model. The special case of a prescribed stretch λ = 1 under an electric field in our model corresponds to zero strain but nonzero stress in the homogeneous solution. The magnitude of the compressive stress in Eq. (42) at λ = 1 is
$|s0|=ϵẼ02$
(75)
From Euler's prediction (Eq. (74)), when $|s0|$ in Eq. (75) increases to the critical value sc in Eq. (74), the dielectric block begins to buckle and the critical nominal electric field for the electromechanical buckling is given by
$Ẽc=π(Kl1/l2)μ3ϵ$
(76)

where the factor K = 0.5 is for both fixed ends, while K = 1 is for both pinned ends of the Euler column.

In contrast to the approximation (Eq. (76)) from Euler's formula, our analytical prediction of the critical nominal electric field from Eq. (72) is obtained as
$Ẽc=(l1mπl2sinhmπl2l1−1)2μϵ$
(77)

where the modes m = 1 and m = 2 are related to two different boundary conditions corresponding K = 1 and K = 0.5 in Eq. (76).

In Fig. 7, we plot the variation of the dimensionless critical electric field $Ẽcϵ/μ$ with respect to the aspect ratio $l1/l2$ from both Euler's approximation (Eq. (76)) and our analytical prediction (Eq. (77)). The critical electric field decreases monotonously with the increase of the aspect ratio $l1/l2$. This trend agrees with intuition that a more slender block (i.e., a larger aspect ratio $l1/l2$) is more likely to become unstable under external stimuli such as an electric field. In the limiting case $l1/l2→∞$, the critical electric field approaches zero and an exceedingly small electric field can make the block buckle.

Fig. 7
Fig. 7
Close modal

### Buckling of a Mechanically Compressed Block.

In contrast to the critical force for Euler's column, the critical stretch (or strain) is often used to define the critical conditions for the buckling of finite blocks or surface instability of soft materials [1,33,3638,]. In Biot's half-space problem [1], the critical stretch for surface instability of a homogeneous neo-Hookean half-space under plane strain compression is 0.544 at which all the wavelengths become unstable. Later, Levinson [37] studied the stability of a compressed block in the current configuration. Recently, Dorfmann and Ogden [36] studied the surface instability of the homogeneous deformation of a half-space subject to both mechanical and electrical loads by solving the incremental boundary-value problem.

In our work, the neo-Hookean block is compressed under plane strain by changing the stretch λ. The critical condition of the buckling is determined by either Eq. (69) for antisymmetric buckling or Eq. (70) for symmetric buckling in the absence of electric fields. The critical stretch λc for the mechanical buckling of the compressed block with different aspect ratios $l1/l2$ is plotted in Fig. 8. The critical stretches for antisymmetric/symmetric buckling with different modes $m=1,2,3,5$ are plotted in solid/dashed lines. In particular, the critical stretches for all modes approach 0.544 when the aspect ratio $l1/l2$ decreases to zero (i.e., $l2/l1$ increases to infinity). The critical stretch $λc=0.544$ of this limiting case ($l2/l1→∞$) coincides with Biot's prediction [1] since the limiting case ($l2/l1→∞$) of block is that of a half-space.

Fig. 8
Fig. 8
Close modal

Figure 8 also shows that the critical stretch for antisymmetric buckling is always larger than that of symmetric buckling. This means that the antisymmetric buckling in a compressed block occurs prior to symmetric buckling. Indeed, symmetric buckling cannot occur unless the passive constraints [37] are considered to be acting until $λ<0.544$. Therefore, only the antisymmetric buckling is compared with Euler's prediction in the preceding discussion.

### Buckling of an Electromechanically Compressed Block.

In Secs. 5.1, 5.2, and 5.3, we have shown that either the mechanical compression in Figs. 6 and 8 or the electric field in Fig. 7 can make the dielectric block buckle. An obvious extension to these notions is that the combined electromechanical loading ought to make buckling of the block yet easier. For purely mechanical loading, the block in our problem can only buckle under compression ($λ<1$) rather than extension ($λ>1$). Since the electric field can make the block buckle at λ = 1 in Fig. 7, the block may also become buckle in extension (i.e., $λ>1$) under an electric field.

Using Eqs. (69) and (70), we plot Fig. 9 the critical stretch λc as a function of the aspect ratio $l1/l2$ for the buckling mode m = 2 under different electric fields $Ẽ0ϵ/μ$. The solid lines denote antisymmetric buckling, while the dashed lines represent symmetric buckling. Note that antisymmetric buckling satisfies all the boundary conditions, while the symmetric buckling satisfies all the boundary conditions other than the electrostatic boundary conditions of the perturbed voltage on the upper and lower surfaces. Furthermore, the critical stretches for the antisymmetric buckling (solid lines) rather than the symmetric buckling (dashed lines) in Fig. 9 are very sensitive to the electric fields. When the aspect ratio $l1/l2$ is larger than five, for example, the differences of the critical stretches for the symmetric buckling between the mechanical compression and the electromechanical loading are negligible. On the other hand, since the occurrence of the symmetric buckling is always later than the onset of antisymmetric buckling, in practice only the effects of the electric fields on antisymmetric buckling are of interest.

Fig. 9
Fig. 9
Close modal

Compared with the critical stretch for buckling of a mechanically compressed block, the critical stretch that accounts for the electric field is shifted upward for a small aspect ratio $l1/l2$. For example, the critical stretch for the buckling of a mechanically compressed block in the limiting case $l1/l2→0$ is 0.544 while it increases to 0.628 at an electric field $Ẽ0ϵ/μ=2$ in Fig. 9.

It is clear from Fig. 9 that the electric field can cause the block to buckle more easily in a compressed state ($λ<1$). Moreover, the electric field can make the block buckle even if the block is in extension ($λ>1$).

We know that both mechanical compression and the electric field can make the block buckle. For antisymmetric buckling with mode m = 2, the variation of the critical stretch λc with respect to the critical electric field $Ẽcϵ/μ$ is plotted in Fig. 10. It is obvious that a slender block (i.e., with high aspect ratio $l1/l2$) is more likely to buckle when it is subject to a combined loading. For example, at a zero electric field (i.e., $Ẽ0ϵ/μ=0$), the critical stretch is slightly less than one in the case of a large aspect ratio $l1/l2=10$, while it approaches 0.544 at an aspect ratio $l1/l2=1$. For each aspect ratio $l1/l2$ in Fig. 10, the critical stretch λc increases monotonically with the increase of $Ẽcϵ/μ$. It clearly shows how the electric field makes the block buckle in an extended state (i.e., $λc>1$). We finally remark that for actual applications, electric breakdown should also be considered and the comparison of the critical electric fields between the electric breakdown and the electrical buckling is needed for a safe design of electrical devices.

Fig. 10
Fig. 10
Close modal

In summary, for a mechanical compression without electric field, the block mathematically exhibits two types of buckling modes, i.e., antisymmetric and the symmetric buckling, however, the antisymmetric buckling will always precede the other. Our results, in the asymptotic limit of large aspect ratio, agree well with Euler's prediction for the buckling of a slender block and, furthermore, at a zero aspect ratio are the same as Biot's critical strain of surface instability of a compressed homogeneous half-space of a neo-Hookean material. For the case where electric fields are included, aside from similar interesting asymptotic connection to Euler's formula, we find that the electric field can cause the block to buckle more easily in a compressed state, and the electric field can even cause the block to buckle in a state of tension.

## Acknowledgment

Financial support from the M.D. Anderson Professorship, NSF CMMI Grant No. 1463339 and NPRP award [NPRP 6-282-2-119] from the Qatar National Research Fund (a member of The Qatar Foundation).

### Appendix A: First Variation of the Energy Functional

The infinitesimal variations of the deformation $x=x(X)$ and the polarization $P̃=P̃(X)$ are denoted, respectively, by $δx$ and $δP̃$, such that
$δx=η1xd, δP̃=η2P̃p$
(A1)

where $η1,η2∈ℝ$ and $max{|η1|,|η2|}≪1$, and $xd$ and $P̃p$ are two smooth variations.

The deformation $x=x(X)$ and the polarization $P̃=P̃(X)$ may not appear directly in the energy functional including the deformation gradient $F=∇x$, the Jacobian $J=det∇x$, the voltage ξ, and the nominal electric displacement $D̃$, among others. Thus, we perform their first variations implicitly and only keep the first-order terms of η1 and η2, such that
$F→F+η1Fd, J→J+η1Jd, ξ→ξ+η1ξd+η2ξp,D̃→D̃+η1D̃d+η2D̃p, E→E+η1Ed+η2Ep$
(A2)
where the subscripts “d” and “p” denote, respectively, the variations related to the deformation and the polarization. For example, $Fd$ and Jd in Eq. (A2) are
$Fd=ddη1∇(x+η1xd)|η1=0=∇xd,Jd=ddη1det(F+η1Fd)|η1=0=JF−T·∇xd$
(A3)
Substituting Eq. (A2) into the Maxwell equation (Eq. (10)) and taking partial derivatives with respect to η1 and η2 at $η1=η2=0$, respectively, we have
$DivD̃d=DivD̃p=0 in ΩR$
(A4)
and the relations
$FTEd+(FT)dE=−∇ξd, FTEp=−∇ξp,FD̃d+FdD̃=ϵ0JEd+ϵ0JdE, FD̃p=ϵ0JEp+P̃p$
(A5)
With Eqs. (A1) and (A2), the variations of Eqs. (4) and (5) are
$xd·e1=0, D̃d·e1=D̃p·e1=0 on Sl∪Sr$
(A6)
$ξd=ξp=0 on Su∪Sb$
(A7)

#### First Variation With Respect to the Polarization

The first variation of the energy functional (Eq. (11)) with respect to the polarization $P̃$ is
$ddη2F[x,P̃+η2P̃p]|η2=0=ddη2U[x,P̃+η2P̃p]|η2=0+ddη2Eelect[x,P̃+η2P̃p]|η2=0=∫ΩR∂Ψ∂P̃·P̃p+ϵ0∫ΩRJE·Ep +∫Su∪Sb(ξpD̃·N+ξD̃p·N)$
(A8)
With Eqs. (A4)(A7) and the divergence theorem, Eq. (A8) becomes
$ddη2F[x,P̃+η2P̃p]|η2=0=∫ΩR(∂Ψ∂P̃·P̃p+ϵ0JE·Ep)+∫∂ΩRξD̃p·N=∫ΩR(∂Ψ∂P̃·P̃p+ϵ0JE·Ep)+∫ΩR(ξDivD̃p+D̃p·∇ξ)=∫ΩR(∂Ψ∂P̃·P̃p+ϵ0JE·Ep)−∫ΩRE·FD̃p=∫ΩR(∂Ψ∂P̃·P̃p+E·(ϵ0JEp−FD̃p))=∫ΩR(∂Ψ∂P̃−E)·P̃p$
(A9)

Based on the basic lemma of calculus of variations, vanishing of Eq. (A9) gives Eq. (14).

#### First Variation With Respect to the Deformation

We introduce a Lagrange multiplier function $q:ΩR→ℝ2$ to address the variation of a constrained problem such that the deformation x is subject to the constraint of incompressibility $J=det∇x=1$. The modified energy functional of Eq. (11), including the Lagrangian multiplier q to enforce incompressibility, is
$F̂[x,P̃]=∫ΩR(Ψ(∇x,P̃)+ϵ02J|E|2−q(J−1))+∫Su∪SbξD̃·N$
(A10)
The first variation of Eq. (A10) with respect to the deformation x is
$ddη1F̂[x+η1xd,P̃]|η1=0=∫ΩR(∂Ψ∂F·∇xd+ϵ02Jd|E|2+ϵ0JE·Ed−qJd) +∫Su∪Sb(ξdD̃·N+ξD̃d·N)$
(A11)
With Eqs. (A3)(A7) and the divergence theorem, Eq. (A11) becomes
$ddη1F̂[x+η1xd,P̃]|η1=0=∫ΩR(∂Ψ∂F·∇xd+ϵ02Jd|E|2+ϵ0JE·Ed−qJd) +∫∂ΩRξD̃d·N=∫ΩR(∂Ψ∂F·∇xd+ϵ02Jd|E|2+ϵ0JE·Ed−qJd) +∫ΩR(ξDivD̃d+D̃d·∇ξ)=∫ΩR(∂Ψ∂F·∇xd+ϵ02Jd|E|2−qJd+E·(ϵ0JEd−FD̃d))=∫ΩR(∂Ψ∂F·∇xd+ϵ02Jd|E|2−qJd+E·(FdD̃−ϵ0JdE))=∫ΩR(∂Ψ∂F+E⊗D̃−ϵ0J2|E|2F−T−qJF−T)·∇xd=∫ΩR(∂Ψ∂F+Σ̃−qJF−T)·∇xd=∫∂ΩRxd·(∂Ψ∂F+Σ̃−qJF−T)N −∫ΩRxd·Div(∂Ψ∂F+Σ̃−qJF−T)$
(A12)

where $Σ̃$ is the Piola–Maxwell stress defined by Eq. (18). Similar derivations of the Piola–Maxwell stress can also be found in the work [42,46] and many other references. With the boundary condition of $xd$ in Eq. (A6), the vanishing of Eq. (A12) gives Eqs. (15)(17).

### Appendix B: Linearized Analysis

Suppose that a deformation x and a polarization $P̃$ have infinitesimal increments $x*$ and $P̃*, ||x*||,||P̃*||≪1$. For a general field $Θ(x,P̃)$ that is (Fréchet-) differentiable at $(x,P̃)$, we have the expansion in the neighborhood of $(x,P̃)$, such that
$Θ(x+x*,P̃+P̃*)=Θ(x,P̃)+∂Θ∂x·x*+∂Θ∂P̃·P̃*+o(||x*||,||P̃*||)$
(B1)
We define
$Θ*=Θ†+Θ‡$
(B2)
where
$Θ†=∂Θ∂x·x*, Θ‡=∂Θ∂P̃·P̃*$
(B3)

Here, $Θ*$ denotes the total linearized increment, and $Θ†$ and $Θ‡$ denote the linearized increments with respect to the deformation and the polarization.

With Eqs. (B1)(B3) and the chain-rule, the linearized increments of the deformation gradient $F=∇x$ and the Jacobian $J=detF=det∇x$ are
$F*=∇x*, (F−T)*=−F−T(∇x*)TF−T,J*=∂J∂F·F*=JF−T·∇x*$
(B4)
Similarly, the linearized increments of other fields at $(x,P̃)$ can be written implicitly as
$(qsξ⋮ED̃Σ̃)Linearization→(q*=q†s*=s†+s‡ξ*=ξ†+ξ‡⋮E*=E†+E‡D̃*=D̃†+D̃‡Σ̃*=Σ̃†+Σ̃‡)$
(B5)

We remark that the linearized increments of the deformation gradient F, the Jacobian J and the Lagrange multiplier q only depend on the increment $x*$ at $(x,P̃)$.

#### Linearized Relation

Consider the linearization of the Maxwell equation (Eq. (10)1) as an example. Substituting the sum of the fields and their linearized increments defined in Eqs. (B4) and (B5) into Eq. (10)1, and ignoring the higher order terms, we obtain
$FTE→[FT+(FT)*](E+E*)→FTE+FTE*+(FT)*E$
(B6a)
$−∇ξ→−∇(ξ+ξ*)=−∇ξ−∇ξ*$
(B6b)
then we have
$FTE*+(FT)*E=−∇ξ*$
(B7)

Other linearized relations can also be obtained in a similar manner.

#### Linearized Piola–Maxwell Stress

Substituting the sum of the fields and their linearized increments defined in Eqs. (B4) and (B5) into the Piola–Maxwell stress (Eq. (18)), and ignoring higher order terms, such that
$→(E+E*)⊗(D̃+D̃*) −ϵ0(J+J*)2|E+E*|2[F−T+(F−T)*]→E⊗D̃−ϵ0J2|E|2F−T+(E⊗D̃*+E*⊗D̃) −ϵ02{2J(E·E*)F−T+|E|2[J(F−T)*+J*F−T]}=Σ̃+Σ̃*$
(B8)

then we have the linearized increment of the Piola–Maxwell stress $Σ̃*$ in Eq. (25).

2

The assumption of homogeneous deformation restricts their analysis essentially to tensile loading to avoid buckling instability.

3

By “physically reasonable,” we imply conditions that are easily realizable in an experimental setup.

4

We explicitly allow for inhomogeneous deformation modes to study buckling under compression.

## References

1.
Biot
,
M.
,
1963
, “
Surface Instability of Rubber in Compression
,”
Appl. Sci. Res.
,
12
(
2
), pp.
168
182
.
2.
Yang
,
S.
,
Khare
,
K.
, and
Lin
,
P.-C.
,
2010
, “
Harnessing Surface Wrinkle Patterns in Soft Matter
,”
,
20
(
16
), pp.
2550
2564
.
3.
Gent
,
A.
, and
Cho
,
I.
,
1999
, “
Surface Instabilities in Compressed or Bent Rubber Blocks
,”
Rubber Chem. Technol.
,
72
(
2
), pp.
253
262
.
4.
Hong
,
W.
,
Zhao
,
X.
, and
Suo
,
Z.
,
2009
, “
Formation of Creases on the Surfaces of Elastomers and Gels
,”
Appl. Phys. Lett.
,
95
(
11
), p.
111901
.
5.
Hohlfeld
,
E.
, and
,
L.
,
2011
, “
Unfolding the Sulcus
,”
Phys. Rev. Lett.
,
106
(
10
), p.
105702
.
6.
Lu
,
N.
, and
Kim
,
D.-H.
,
2014
, “
Flexible and Stretchable Electronics Paving the Way for Soft Robotics
,”
Soft Rob.
,
1
(
1
), pp.
53
62
.
7.
Shian
,
S.
,
Bertoldi
,
K.
, and
Clarke
,
D. R.
,
2015
, “
Dielectric Elastomer Based ‘Grippers’ for Soft Robotics
,”
,
27
(
43
), pp.
6814
6819
.
8.
Rogers
,
J. A.
,
Someya
,
T.
, and
Huang
,
Y.
,
2010
, “
Materials and Mechanics for Stretchable Electronics
,”
Science
,
327
(
5973
), pp.
1603
1607
.
9.
Shankar
,
R.
,
Ghosh
,
T. K.
, and
Spontak
,
R. J.
,
2007
, “
Dielectric Elastomers as Next-Generation Polymeric Actuators
,”
Soft Matter
,
3
(
9
), pp.
1116
1129
.
10.
Moscardo
,
M.
,
Zhao
,
X.
,
Suo
,
Z.
, and
Lapusta
,
Y.
,
2008
, “
On Designing Dielectric Elastomer Actuators
,”
J. Appl. Phys.
,
104
(
9
), p.
093503
.
11.
Keplinger
,
C.
,
Kaltenbrunner
,
M.
,
Arnold
,
N.
, and
Bauer
,
S.
,
2010
, “
Röntgen's Electrode-Free Elastomer Actuators Without Electromechanical Pull-in Instability
,”
,
107
(
10
), pp.
4505
4510
.
12.
Koh
,
S. J. A.
,
Zhao
,
X.
, and
Suo
,
Z.
,
2009
, “
Maximal Energy That Can Be Converted by a Dielectric Elastomer Generator
,”
Appl. Phys. Lett.
,
94
(
26
), p.
262902
.
13.
Bauer
,
S.
,
Bauer-Gogonea
,
S.
,
Graz
,
I.
,
Kaltenbrunner
,
M.
,
Keplinger
,
C.
, and
Schwödiauer
,
R.
,
2014
, “
25th Anniversary Article: A Soft Future: From Robots and Sensor Skin to Energy Harvesters
,”
,
26
(
1
), pp.
149
162
.
14.
Zhao
,
X.
, and
Wang
,
Q.
,
2014
, “
Harnessing Large Deformation and Instabilities of Soft Dielectrics: Theory, Experiment, and Application
,”
Appl. Phys. Rev.
,
1
(
2
), p.
021304
.
15.
Stark
,
K.
, and
Garton
,
C.
,
1955
, “
,”
Nature
,
176
(
4495
), pp.
1225
1226
.
16.
Plante
,
J.-S.
, and
Dubowsky
,
S.
,
2006
, “
Large-Scale Failure Modes of Dielectric Elastomer Actuators
,”
Int. J. Solids Struct.
,
43
(
25
), pp.
7727
7751
.
17.
Zhao
,
X.
, and
Suo
,
Z.
,
2009
, “
Electromechanical Instability in Semicrystalline Polymers
,”
Appl. Phys. Lett.
,
95
(
3
), p.
031904
.
18.
Wang
,
Q.
, and
Zhao
,
X.
,
2013
, “
Creasing-Wrinkling Transition in Elastomer Films Under Electric Fields
,”
Phys. Rev. E
,
88
(
4
), p.
042403
.
19.
Wang
,
Q.
,
Zhang
,
L.
, and
Zhao
,
X.
,
2011
, “
Creasing to Cratering Instability in Polymers Under Ultrahigh Electric Fields
,”
Phys. Rev. Lett.
,
106
(
11
), p.
118301
.
20.
Wang
,
Q.
,
Suo
,
Z.
, and
Zhao
,
X.
,
2012
, “
Bursting Drops in Solid Dielectrics Caused by High Voltages
,”
Nat. Commun.
,
3
, p.
1157
.
21.
Ha
,
S. M.
,
Yuan
,
W.
,
Pei
,
Q.
,
Pelrine
,
R.
, and
Stanford
,
S.
,
2006
, “
Interpenetrating Polymer Networks for High-Performance Electroelastomer Artificial Muscles
,”
,
18
(
7
), pp.
887
891
.
22.
Zhao
,
X.
, and
Suo
,
Z.
,
2007
, “
Method to Analyze Electromechanical Stability of Dielectric Elastomers
,”
Appl. Phys. Lett.
,
91
(
6
), p.
061921
.
23.
Kofod
,
G.
,
2008
, “
The Static Actuation of Dielectric Elastomer Actuators: How Does Pre-Stretch Improve Actuation?
,”
J. Phys. D: Appl. Phys.
,
41
(
21
), p.
215405
.
24.
Li
,
B.
,
Zhou
,
J.
, and
Chen
,
H.
,
2011
, “
Electromechanical Stability in Charge-Controlled Dielectric Elastomer Actuation
,”
Appl. Phys. Lett.
,
99
(
24
), p.
244101
.
25.
Akbari
,
S.
,
Rosset
,
S.
, and
Shea
,
H. R.
,
2013
, “
Improved Electromechanical Behavior in Castable Dielectric Elastomer Actuators
,”
Appl. Phys. Lett.
,
102
(
7
), p.
071906
.
26.
Niu
,
X.
,
Stoyanov
,
H.
,
Hu
,
W.
,
Leo
,
R.
,
Brochu
,
P.
, and
Pei
,
Q.
,
2013
, “
Synthesizing a New Dielectric Elastomer Exhibiting Large Actuation Strain and Suppressed Electromechanical Instability Without Prestretching
,”
J. Polym. Sci. Part B: Polym. Phys.
,
51
(
3
), pp.
197
206
.
27.
Jiang
,
L.
,
Betts
,
A.
,
Kennedy
,
D.
, and
Jerrams
,
S.
,
2016
, “
Eliminating Electromechanical Instability in Dielectric Elastomers by Employing Pre-Stretch
,”
J. Phys. D: Appl. Phys.
,
49
(
26
), p.
265401
.
28.
Keplinger
,
C.
,
Li
,
T.
,
Baumgartner
,
R.
,
Suo
,
Z.
, and
Bauer
,
S.
,
2012
, “
Harnessing Snap-Through Instability in Soft Dielectrics to Achieve Giant Voltage-Triggered Deformation
,”
Soft Matter
,
8
(
2
), pp.
285
288
.
29.
Shivapooja
,
P.
,
Wang
,
Q.
,
Orihuela
,
B.
,
Rittschof
,
D.
,
López
,
G. P.
, and
Zhao
,
X.
,
2013
, “
Bioinspired Surfaces With Dynamic Topography for Active Control of Biofouling
,”
,
25
(
10
), pp.
1430
1434
.
30.
Díaz-Calleja
,
R.
,
Riande
,
E.
, and
Sanchis
,
M.
,
2008
, “
On Electromechanical Stability of Dielectric Elastomers
,”
Appl. Phys. Lett.
,
93
(
10
), p.
101902
.
31.
Xu
,
B.-X.
,
Mueller
,
R.
,
Klassen
,
M.
, and
Gross
,
D.
,
2010
, “
On Electromechanical Stability Analysis of Dielectric Elastomer Actuators
,”
Appl. Phys. Lett.
,
97
(
16
), p.
162908
.
32.
Bertoldi
,
K.
, and
Gei
,
M.
,
2011
, “
Instabilities in Multilayered Soft Dielectrics
,”
J. Mech. Phys. Solids
,
59
(
1
), pp.
18
42
.
33.
Dorfmann
,
L.
, and
Ogden
,
R. W.
,
2014
, “
Instabilities of an Electroelastic Plate
,”
Int. J. Eng. Sci.
,
77
, pp.
79
101
.
34.
Leng
,
J.
,
Liu
,
L.
,
Liu
,
Y.
,
Yu
,
K.
, and
Sun
,
S.
,
2009
, “
Electromechanical Stability of Dielectric Elastomer
,”
Appl. Phys. Lett.
,
94
(
21
), p.
211901
.
35.
Suo
,
Z.
,
2010
, “
Theory of Dielectric Elastomers
,”
Acta Mech. Solida Sin.
,
23
(
6
), pp.
549
578
.
36.
Dorfmann
,
A.
, and
Ogden
,
R.
,
2010
, “
Nonlinear Electroelastostatics: Incremental Equations and Stability
,”
Int. J. Eng. Sci.
,
48
(
1
), pp.
1
14
.
37.
Levinson
,
M.
,
1968
, “
Stability of a Compressed Neo-Hookean Rectangular Parallelepiped
,”
J. Mech. Phys. Solids
,
16
(
6
), pp.
403
408
.
38.
Triantafyllidis
,
N.
,
Scherzinger
,
W.
, and
Huang
,
H.-J.
,
2007
, “
Post-Bifurcation Equilibria in the Plane-Strain Test of a Hyperelastic Rectangular Block
,”
Int. J. Solids Struct.
,
44
(
11
), pp.
3700
3719
.
39.
Kankanala
,
S.
, and
Triantafyllidis
,
N.
,
2008
, “
Magnetoelastic Buckling of a Rectangular Block in Plane Strain
,”
J. Mech. Phys. Solids
,
56
(
4
), pp.
1147
1169
.
40.
Golubitsky
,
M.
, and
Schaeffer
,
D. G.
,
1985
,
Singularities and Groups in Bifurcation Theory
, Vol. 1,
Springer
,
Berlin
.
41.
Chen
,
Y.-C.
,
2001
, “
Singularity Theory and Nonlinear Bifurcation Analysis
,”
Nonlinear Elasticity: Theory and Applications
, Y. B. Fu and R. W. Ogden, eds.,
Cambridge University Press
,
Cambridge. UK
.
42.
Liu
,
L.
,
2014
, “
An Energy Formulation of Continuum Magneto-Electro-Elasticity With Applications
,”
J. Mech. Phys. Solids
,
63
, pp.
451
480
.
43.
Suo
,
Z.
,
Zhao
,
X.
, and
Greene
,
W. H.
,
2008
, “
A Nonlinear Field Theory of Deformable Dielectrics
,”
J. Mech. Phys. Solids
,
56
(
2
), pp.
467
486
.
44.
James
,
R. D.
, and
Kinderlehrer
,
D.
,
1990
, “
Frustration in Ferromagnetic Materials
,”
Continuum Mech. Thermodyn.
,
2
(
3
), pp.
215
239
.
45.
Shu
,
Y. C.
, and
Bhattacharya
,
K.
,
2001
, “
Domain Patterns and Macroscopic Behaviour of Ferroelectric Materials
,”
Philos. Mag. Part B
,
81
(
12
), pp.
2021
2054
.
46.
Deng
,
Q.
,
Liu
,
L.
, and
Sharma
,
P.
,
2014
, “
Flexoelectricity in Soft Materials and Biological Membranes
,”
J. Mech. Phys. Solids
,
62
, pp.
209
227
.