Branching dynamics in electrohydrodynamic instabilities of viscoelastic soft gels

Abstract

An electric field imposed on a bilayer of fluids that are stably stratified in the presence of gravity leads to an instability manifested by interfacial deflections. The layer of perfect conductor is simulated using a linear viscoelastic model and the perfect dielectric is considered to be a layer of air. Under neutral conditions, the key dimensionless groups are the dimensionless electric potential, Bond number and the Weissenberg number. The branching behavior upon instability to sinusoidal disturbances is determined by weak nonlinear analysis with the dimensionless potential advanced from its critical value at neutral stability. An analytical expression obtained from weak nonlinear analysis leads to the unintuitive result that sinusoidal deflections can either lead to supercritical saturated waves or lead to subcritical breakup depending on the elasticity of the perfect conductor. The analytical expression also indicates that there is a transition wave number below which supercritical saturation ought to occur, it can be shown that such wave numbers can be geometrically accessed, thus permitting any supercritical saturation to steady waves. In contrast, our results demonstrate that when the perfect conductor is modeled as an Oldroyd-B fluid, the branching remains subcritical in nature, ultimately leading to interface rupture—mirroring the behavior observed in the Newtonian fluid case (as demonstrated by B. Dinesh and R. Narayanan, Phys. Rev. Fluids, 2021, 6, 054001).

---

1 Introduction and physics

Electrohydrodynamics concerns the interaction between fluids and electric fields. The foundational studies by Taylor and Melcher^2,3 established the theoretical basis for analyzing such phenomena. Subsequent research has examined how electric fields influence the morphology and stability of interfaces in systems such as liquid bridges, drops, jets, and thin films. When an electric field is applied perpendicular to an initially flat interface, it typically induces destabilization.

Vertical electric fields have been applied in quiescent multilayer flows for the controlled self-assembly of microscale hierarchical structures in polymer melts and processes like soft lithography.^4–8 Other applications include inkjet printing,⁹ enhancement of mixing and heat transfer, as well as the development of soft devices with adjustable shapes.^10,11 For an overview of electrohydrodynamic instability in thin fluid films, see the review by Papageorgiou.¹²

In this work, we investigate the branching behavior of electrohydrodynamic instability in a system comprising a perfect conductor–dielectric fluid pair. The perfect conductor is modeled in two ways: first, using a linear viscoelastic model, and second, using the Oldroyd-B fluid model. The dielectric fluid is assumed to be a hydrodynamically passive air layer. As shown in Fig. 1, the fluids are confined between two rigid plates, with the bottom plate held at a constant voltage D, while the top plate is grounded. The applied voltage counteracts gravity and drives an instability that leads to the formation of an array of pillars in the conducting fluid.^1,12–15 However, when the elastic stresses in the conducting fluid overcome the destabilizing effect of the constant voltage, the interface undergoes saturated deformation without branching.

Fig. 1 Schematic of electrostatic instability. The bottom viscoelastic fluid is a perfect conductor and the top fluid is a perfect dielectric. The bottom plate is maintained a constant voltage D and the top plate is grounded.

The instability arises due to competition between the applied potential and stabilizing factors such as gravity, elastic stresses, and surface tension. These competing effects can produce a minimum in the voltage–wave number curve at the onset of instability.^1,16

At neutral stability, if the perfect conductor is modeled as a Newtonian fluid, the velocity perturbations vanish. For a viscoelastic perfect conductor, at low wave numbers, the destabilizing voltage is balanced by gravity and elastic stresses, while at high wave numbers, surface tension opposes the destabilizing effect.^1,16,17 The high-wave number behavior is similar to that seen in the Rayleigh–Taylor instability, which is typically subcritical in nature. In contrast, the low-wave number balance resembles that of thermocapillary-driven problems, such as the Bénard–Marangoni instability, which generally exhibits supercritical behavior.¹⁸

The objective of this study is to determine whether a transition occurs between supercritical and subcritical behavior at a critical wave number, or if the instability consistently exhibits one of these two behaviors across all wave numbers.

We review previous studies relevant to this work. Wu and Chou¹⁹ showed that, for Maxwell fluids subjected to DC electric fields, elasticity increases the growth rate of pillar formation. They identified a critical Deborah number above which the growth rate for certain wavenumbers diverges.

Tomar et al.²⁰ examined the surface instability of a confined viscoelastic liquid film under an applied electric field, considering both Maxwell and Jeffreys models. Their analysis showed that incorporating fluid inertia removes the singularity, resulting in large but finite growth rates for all Deborah numbers and identifying the dominant instability wavelength. They also observed that the small inertia limit does not coincide with the zero-inertia limit when describing the instability dynamics and wavelength in polymer melts.

Espín et al.⁷ studied the effects of viscoelasticity on instabilities under AC and DC electric fields using the Jeffreys model for both perfect and leaky dielectrics. In the DC case, asymptotic methods were applied to resolve a singularity occurring when solvent viscosity is neglected, corresponding to the Maxwell limit. Their results show that elasticity increases both the maximum growth rate and the wavenumber of the instability.

Dinesh et al.¹ investigated the electrohydrodynamic instability between a perfect conductor and a dielectric, modeling both as Newtonian fluids. Their study found that the branching behavior is always subcritical, resulting in rupture of the interface. In contrast, the present study demonstrates that both subcritical and supercritical branching behavior can occur when the perfect conductor is described using a viscoelastic model.

Observations from previous studies indicate that pillar formation due to instability spans the gap between the plates, which is characteristic of subcritical interface behavior. However, elasticity in the viscoelastic fluid can also result in supercritical saturation of the interface. The objective of this work is to demonstrate this behavior through mathematical analysis and to identify the physical conditions that lead to either subcritical or supercritical instability.

The analysis presented here follows the methodology outlined in the works of Dinesh and Narayanan.^1,21 We consider a system comprising a perfect conductor–dielectric fluid pair and derive an analytical expression to characterize the branching behavior. The perfect conductor is modeled using a linear viscoelastic framework. Additionally, we examine a case where the perfect conductor is represented by an Oldroyd-B fluid.

A perturbation expansion is carried out around the neutral stability state, with the applied voltage slightly exceeding its critical threshold. This approach allows us to investigate the nonlinear evolution of the interface and distinguish between subcritical and supercritical regimes.

2 Mathematical model

The model assumes a hydrodynamically active viscoelastic fluid with constant properties whose free surface is exposed to a passive gas in a stabilizing gravitational field, as depicted in Fig. 1. The gas layer is considered to be air and hydrodynamically inactive. The viscoelastic fluid and the air layers are between two rigid electrically conducting plates located at −h* and h, across which a constant voltage difference, D, is applied (cf.Fig. 1). The potential field in the top fluid is given by

∇²ψ = 0 (1)

The bottom viscoelastic fluid, represented by an asterisk, is assumed to be a perfect conductor while the top fluid is taken to be a perfect dielectric. The bottom fluid is taken to be linearly viscoelastic for algebraic simplicity whilst retaining essential physics in this study. Now, the stress tensor in a viscoelastic fluid may be expressed in terms of the displacement field. The displacement vector, R, is the displacement of the position vector, x, in the current configuration from the position vector, ζ, in the reference configuration. In other words,

x = ζ + R(x) (2)

For a linear viscoelastic fluid (cf., ref. 22–24) the stress tensor, T, is given by

T = −pI + G(∇R + ∇R^T) + μ_g(∇v + ∇v^T) (3)

Here G and μ_g are the shear modulus and viscosity of the fluid material and v is the velocity field. Now the velocity field in the fluid is itself expressed in terms of the displacement field (cf., ref. 25 for a detailed explanation). This expression is given by

In a linear viscoelastic soft-gel layer, the elastic field is described by the displacement field vector (cf.eqn (2)), which represents the deformation of material points, while the velocity field vector in eqn (4), describes the rate of change of these displacements over time. The relaxation mechanism occurs as the internal elastic stresses, generated by deformation, gradually dissipate through viscous dissipation in the soft-gel. This means the displacement field changes slowly, and the velocity field captures how the material returns toward its undeformed state by relaxing the elastic stresses. Over time, energy stored in elastic deformation is released or converted to viscous dissipation, resulting in a decrease of internal stresses and a gradual return to equilibrium. This relaxation process is characterized by time-dependent decay of stresses and strains governed by the viscoelastic properties of the gel.^{21,23,24,26–29} The equations of motion in the viscoelastic medium are thus

where T and v are given by eqn (3) and (4), where in eqn (4) the superscript T on ∇R is the transpose of the tensor ∇R and I is the identity tensor and where i_z is the base vector in the positive z-direction in eqn (5). The horizontal and vertical components of the displacement vector, R, are denoted by X and Z.

The viscoelastic fluid is taken to be incompressible and the mass conservation demands (ref. 25 and 30)

det(F) = 1 (6)

where F is the deformation tensor given by . Here ζ_i represent the components of the position vector in the reference configuration. Upon expansion of eqn (6) and using eqn (2) we geteqn (5) and (7) along with the representations for T and v constitute the domain equations. Other models, such as the neo-Hookean^25,29,31 and Oldroyd-B models,³² can also describe the soft-gel layer. In this work, we use the linear viscoelastic model to obtain a mathematically tractable solution and focus on the core physics of the electrohydrodynamic instability. A stability analysis using the Oldroyd-B model is included in the Appendix D. The governing equations are complemented by boundary conditions at the rigid wall and at the soft-gel air interface.

At the wall, the displacement fields are taken to be zero. At the interface, z = h(x, t), the normal and the tangential components of the momentum balance hold. They are

n·T·t = 0  and  n·T·n − D(T^M·n) = −γ∇·n (8)

where the unit normal vector (n) and the unit tangent vector (t) are given byIn addition we have an impermeable interface along with its kinematic relationIn eqn (8), the dimensional form of the Maxwell stress tensor, i.e., T^M, (cf. ref. 33), is given by where E = −∇ψ is the electric field, ε is the relative permittivity of the fluid and ε₀ is the permittivity of free space.

The governing equations are made dimensionless by using the following scales denoted by the subscript ‘c’:

Here U is a characteristic velocity scale. The scaled potential field in the top fluid is given by

∇²ψ = 0 (12)

The potential field is subject to a constant value of unity at the bottom plate, i.e., at z = −1, while the top plate at z = [script letter H] is assumed to be grounded, where [script letter H] = h/h*.

The dimensionless governing equations for the bottom fluid are given by

where , and and . The mass conservation eqn (7) remains unchanged and does not induce any dimensionless groups. The interfacial force balance conditions at z = h(x, t) becomeThe normal component of the Maxwell stress in dimensionless form can be written asHere, [scr D, script letter D] = εε₀D²/h*μ*U. For the perfect conductor–dielectric model, the tangential components of the Maxwell stress tensor vanish (cf. ref. 33). From the scaling of the equations and boundary conditions, five key dimensionless groups emerge: Re, Ca, Bo, [scr D, script letter D], and Wi, which are evaluated using representative physical parameters (see Tables 1 and 2). The base state R = v = T = 0 is subject to a stability analysis to explore the onset and nature of the bifurcation. We seek the neutral stability condition and characterize steady-state branching, governed primarily by Wi, Ca, Bo, and [scr D, script letter D]. Notably, under steady-state conditions, Ca appears only as Ca/Wi, allowing U to be chosen such that Ca = 1.

Table 1 Physical properties of polyacrylamide viscoelastic fluid^34,35

Physical property	Range
Density (ρ)	1000 kg m^−3
Viscosity (μ)	1.5 kg m^−1 s^−1
Shear modulus (G)	10–45 Pa
Surface tension (γ)	73 × 10^−3 N m^−1
Thickness (h*)	1 × 10^−3–20 × 10^−3 m

---

Table 2 Range of key dimensionless quantities for polyacrylamide viscoelastic fluid. The velocity scale U is based on the capillary number such that

Dimensionless parameter	Definition	Range
Re	ρh*U/μ	1.5
Wi	μU/Gh*	0.1–10
Bo	ρgh*^2/γ	0.1–50

---

3 Linear stability analysis

A linear stability analysis is carried out by introducing small perturbations ψ₁, R₁, p₁, and ζ₁ about the quiescent base state. These are expressed using the expansion

ψ(x, z, t) = ψ₀(z) + ψ₁(z)cos(kx)e^σt,

with similar forms assumed for R, p, and ζ. Here, k denotes the wavenumber of the perturbation and the subscript “0” refers to the base-state quantities. Neutral stability corresponds to σ = 0 in the above expansion.

Assuming constant fluid volume imposes the constraint

is the disturbance wavelength.

Since the base-state displacement R₀ = 0, the base-state pressure gradient is balanced by the gravitational body force:

here [capital G, script] = ρgh*²/μU and δ^j = ρ^j/ρ, takes values ρ_dieltric/ρ_conductor and 1. The base-state potential field is given by ψ₀ = 1 − z/[script letter H]. The perturbation fields, being non-zero, are determined by solving the linearized governing equations for the soft gel. Substituting the assumed forms of the disturbances into the governing equations yields:

Linearizing the boundary conditions at the reference interface, z = 0

Z₁ = ζ₁ (19)

The no displacement fields X₁ = 0 and Z₁ = 0, boundary conditions, at the wall i.e., z = −1 are imposed. Upon eliminating p₁ and X₁ from the governing equations and using the boundary conditions we get an expression for Z₁, i.e.,The linearized equation for the potential field in the top fluid are given bywhere at z = 0 and z = [script letter H], we have . Upon solving eqn (22), we getTo determine the neutral stability criteria, we turn toward the perturbed normal component of the momentum balance along the interface, z = 0. First, we take the derivative of this equation with respect to x to eliminate pressure, p₁ using the x-momentum equation. Second, using the continuity equation, we eliminate X₁ in terms of Z₁, i.e.,Upon substituting the solution for ψ₁ and Z₁, in eqn (24), we getHere, [capital G, script]Ca (also known as the Bond number, Bo), Wi/Ca and [scr D, script letter D]Ca, hereafter termed [scr D, script letter D]₀Ca at the neutral point, can be grouped as pairs and are observed to be independent of the characteristic velocity scale, U, and the viscosities. In addition, here δ is the ratio of the density of the perfect dielectric to the density of the perfect conductor.Observe that for low wavenumbers and shallow film thicknesses, term I is of and term II is of [scr O, script letter O](1). This is obtained by expanding cosh(2k), sinh(k), and cosh(k) in term I in series of the wavenumber k and neglecting higher order terms of wavenumbers. For example, coth(k[script letter H]) when expressed in a series of wavenumbers can be written as 1/k[script letter H]. Upon multiplying eqn (27) by k², the Bo term is of [scr O, script letter O](k²), the term consisting of is of [scr O, script letter O](1), the curvature term is of [scr O, script letter O](k⁴), and the term containing [scr D, script letter D]Ca is of [scr O, script letter O](k²). This indicates that for low wavenumber, the stabilizing effect of elastic stresses is dominant over the stabilizing effects of gravity, surface tension, and the destabilizing effect of the potential.

It is noteworthy that the term associated with the elasticity (Wi) of the perfect conductor has an opposite sign compared to the term associated with the critical potential ([scr D, script letter D]). This indicates that the elastic stresses, like gravity, cause a stabilizing effect, whereas the potential leads to a destabilizing effect on the interfacial deformations. It is noteworthy that when a similar linear stability analysis is performed by employing the Oldroyd-B model to capture the behavior of the soft-gel layer, the following expression relating Bo and [scr D, script letter D]Ca is obtained (cf. Appendix D for the derivation), i.e.,

This expression matches the critical potential for a Newtonian fluid (perfect conductor) and a hydrodynamically passive gas (perfect dielectric).¹ The elasticity of the Oldroyd-B fluid does not change the critical potential for the onset of instability. As a result, the Oldroyd-B model does not affect the outcome of the linear stability analysis under neutral conditions. This result contrasts with our current findings, where decreasing the Weissenberg number, Wi, raises the critical potential required for instability. This is the main reason for using the linear viscoelastic model to represent the soft-gel layer in both the linear stability and weakly nonlinear analyses.

Typical neutral stability curves obtained from eqn (27) are shown in Fig. 2 drawn for the example of air on top of polyacrylamide gel. Observe that on the falling branch of the neutral stability curve, i.e., for low wave numbers the electrostatic potential driving the instability is balanced by gravity and elastic stresses, while on the rising branch, i.e., for high wave numbers the voltage applied across the horizontal conducting walls is stabilized by interfacial tension. Interfacial tension is a force acting tangentially along the interface but creates a net restoring force normal (perpendicular) to the interface through curvature effects. The force pulls the interface inward to minimize its surface area and restore stability. Interfacial tension stabilizes the interface by acting as a restoring force that resists deformation, minimizing the surface area between two phases and suppressing disturbances that would increase irregularities. This helps maintain a smooth and stable boundary. Interfacial tension enters the mathematical formulation through the normal stress balance equation in eqn (8) (dimensional form) and via the Capillary number Ca in eqn (14). This stabilizing effect of interfacial tension for high wave numbers is analogous to Rayleigh–Taylor instability of a fluid, wherein the instability is always subcritical in nature. The balance between potential, elastic stresses and gravity for low wave numbers is reminiscent of the Bénard instability which is supercritical in nature. Note that for high Weissenberg numbers i.e., Wi = 1 and Wi = 10, the neutral curve tends towards the critical curve corresponding to that of a Newtonian perfect conductor, i.e., the curve with the solid dots (cf.Fig. 2(a) of Dinesh et al.¹). For low Weissenberg numbers, i.e., Wi = 0.1, the critical potential [scr D, script letter D]₀Ca for the onset of an instability increases. This increment in the critical potential is attributed to the stabilizing elastic stresses in the soft-gel layer.

Fig. 2 [scr D, script letter D]₀Ca vs. k² obtained from eqn (26). The case of air on top of water is assumed. Figure (a) is drawn for [script letter H] = 1, Bo = 5 and Wi = 0.1, 1 and 10. Note that the solid dots in the figure correspond to the case where the perfect conductor is a Newtonian fluid. Figure (b) is drawn for [script letter H] = 1, Bo = 5, and Wi = 0.1 taking , with n being the horizontal mode index.

Not all wave numbers in Fig. 2a are accessible. When a minimum appears in the curves, the accessible wave numbers, k are related to the horizontal modal index, n, as shown in Fig. 2b for a one-dimensional system of width, w, where . where n takes integer values 1, 2, 3,…, with n = 1 representing one wave at the interface, n = 2 representing two waves, and so on. Here, w is the width of the container confining the perfect conductor and perfect dielectric between flat plates. To obtain the curve of the critical dimensionless potential versus the wavenumber, the waveform (i.e., the value of n) is fixed while the width “w” of the container is varied. This variation results in a continuous and smooth curve for the dimensionless potential as a function of the wavenumber. There are specific values of “w” where two consecutive modes can coexist, and around these points, the curves are non-monotonic. These points and the wave numbers between them are crucial to our study's conclusions. Our aim is to determine whether the instability is supercritical or subcritical for each accessible wave number by performing a weak nonlinear analysis around the critical state.

4 Weak nonlinear analysis and discussion: linear viscoelastic model

A weak nonlinear analysis is used to study how an interface behaves just beyond the onset of an instability, whether it stabilizes at a finite amplitude (supercritical saturation) or grows toward rupture (subcritical rupture). For the case of the perfect conductor behaving like a linear viscoelastic fluid, the ensuing analysis follows the procedure outlined in the work of Dinesh and Narayanan.^1,21 We start with a linear stability analysis to find the critical conditions, the governing equations are then expanded in a small parameter (ε, given below in eqn (29)) that measures deviation from this threshold. Next, we solve the equations order by order, i.e., [scr O, script letter O](ε²) and [scr O, script letter O](ε³). This leads to an equation in terms of the amplitude of the interfacial deformation (which we refer as ζ₁ or [scr A, script letter A]). To determine the nature of the bifurcation, the dimensionless potential, [scr D, script letter D], is advanced from its critical value, [scr D, script letter D]₀ by an amount ε, such that ε is defined byThe response of the potential to an increase of the control variable [scr D, script letter D] is given byThe displacement fields, pressure field and the interface deflection change from their base state byandAt the interface, the interior displacement field, X, is expressed as ref. 36and likewise for Z(x, z) and p(x, z). The equations at various orders can thus be obtained. The equations at [scr O, script letter O](ε) are given byThe boundary conditions at [scr O, script letter O](ε) at the rigid wall, i.e., z = −1 are

Z₁ = 0 and X₁ = 0 (39)

The interfacial conditions, i.e., at z = 0 at [scr O, script letter O](ε) become

Z₁ − ζ₁ = 0 (40)

andObserve that the [scr O, script letter O](ε) equations are precisely the neutral stability equations. Hence ζ₁(x) becomes ζ₁(x) = [scr A, script letter A] cos(kx). Our job now is to determine the sign of [scr A, script letter A]², noting that a positive value of [scr A, script letter A]² implies a supercritical bifurcation and a negative value implies a subcritical bifurcation. When the system consisting of the perfect conductor and perfect dielectric becomes unstable, interfacial deformations start to grow. If the square of the amplitude, i.e., [scr A, script letter A]², is positive, the interfacial deformations grow gently and the amplitude of these deformations increases smoothly as we go past the critical potential—this is a supercritical bifurcation. This finally leads to the formation of saturated interfacial deformation. If the square of the amplitude, i.e., [scr A, script letter A]², is negative, the interface resists disturbances but then suddenly jumps to a large amplitude deformation once pushed away from the critical potential, often showing hysteresis—this is a subcritical bifurcation. In short, supercritical bifurcation means a smooth transition to a saturated interfacial deformation, while subcritical bifurcation means an abrupt jump leading to the rupture of the interface, which ultimately leads to pillars of the perfect conductor spanning the gap between the plates.¹² To determine [scr A, script letter A]², we proceed to the next order.

The governing equations at are given by

and

The boundary conditions at z = −1 are

Z₂ = 0  and  X₂ = 0 (46)

and the interfacial conditions, i.e., at z = 0, are nowandThe potential field at is given by

∇²ψ₂ = 0 (50)

subject to the following conditionsObserve the pattern of eqn in (36)–(42) and compare them with eqn (43)–(51). The right-hand sides are comprised of forcing terms that are quadratic combinations of [scr O, script letter O](ε) terms and directly proportional to [scr A, script letter A]². We call these quadratic combinations (1, 1) terms due to their bi-linear combinations of first order variables. The solution to the second order problem must therefore be the sum of several [scr A, script letter A]² dependent terms. In addition, the boxed term in eqn (51) evolves due to the expansion at . This term leads to the determination of [scr A, script letter A]² at . Upon observing that eqn (36)–(42) are homogeneous and employing solvability conditions on eqn (36)–(42), we see that solvability of eqn (43)–(51) is automatically satisfied. Thus we cannot determine [scr A, script letter A]² at this order and must advance to the next order. In weak nonlinear analysis, the solvability condition ensures that the perturbation expansions remain bounded by preventing unphysical growth of perturbations. Applied to the normal stress balance equation at the interface, this condition requires that higher-order nonlinear terms satisfy a compatibility criterion, correlating factors like disturbance amplitude, surface tension, elasticity, and electrostatic potential. This balance leads to an amplitude equation that governs the evolution and saturation of interfacial deformations beyond the linear stability threshold, ensuring a physically consistent description of the interface evolution.

The governing equations at are given by

andThe boundary conditions at z = −1 are

Z₃ = 0 and  X₃ = 0 (55)

The interfacial conditions at z = 0 are

and

At the third order the potential field is governed by

∇²ψ₃ = 0 (59)

subject to ψ₃ = 0 at z = [script letter H] and the following condition at z = 0

In the absence of electrostatic potential and upon reversing the direction of gravity, the governing equations for the weakly nonlinear regime reduce to the classical Rayleigh–Taylor instability problem. In this limiting case, the Bond number advances beyond its critical value, i.e., Bo = Bo₀ + ε²/2. In Rayleigh–Taylor instability, gravity destabilizes the system when a denser fluid overlies a lighter fluid (cf. ref. 21 for weak nonlinear analysis of Rayleigh–Taylor instability of a soft-gel layer). In our current study, the destabilization arises from Maxwell stresses at the soft-gel fluid interface, which causes interfacial deformations. However, the stabilization of the interfacial deformations arises from elastic stresses, gravity and surface tension. It is noteworthy that from eqn (27) term II arises from the Maxwell stresses in the normal force balance equation. This term is solely responsible for driving an instability when a critical potential is applied across the perfect conductor and perfect dielectric fluid layers. However, this term is balanced by the Bond number term, term I, and surface tension. The Bond number term indicates that gravity opposes the potential, term I shows that elastic stresses act against the potential, and the surface tension term also acts in opposition to the potential.

If [scr A, script letter A]² is positive, the branch is supercritical, but if it is negative, [scr D, script letter D] would not be advanced as indicated by eqn (29) but reduced instead and the branch would become a subcritical pitchfork. The calculations that involve solvability require the use of symbolic manipulation, which was carried out in Mathematica®. The algebraic complications can be reduced if the nonlinear terms of eqn (7) are dropped.²¹ It has been observed by the present authors that the results do not change qualitatively.

4.1 Solution of the nonlinear equations and tracing the cause of the transition from super to subcritical branching

The solution to the first order governing equations is obtained from the linear stability analysis discussed in the Section 3. We now proceed to seek the solution for the governing equation at . At this second order, the displacement fields, pressure field and interface deformation must be expressed asHere the displacement fields and the pressure field consist of an x-dependent part and an x-independent part. The domain equations for the x-dependent part of the problem areandAs before, upon eliminating from the governing eqn (63) and (64), by using the continuity eqn (62), we getAt z = 0, we have the tangential stress and kinematic conditions. Eliminating from the x-dependent part of the tangential stress condition, eqn (48), by taking the horizontal derivative and using the continuity equation, (62), we getwhile the x-dependent part of kinematic condition, (47), at z = 0 givesSolving eqn (65) and applying the tangential stress condition, (66), along with the kinematic condition, (67), we getEmploying eqn (68) in the continuity eqn (62), and the x-momentum eqn (63), yieldsand

At the second order we are left with obtaining the solutions for the x-independent parts of X₂, Z₂ and p₂. To do this, we eliminate pressure from the momentum eqn (44) and (45), and further eliminate X₂ by using the continuity eqn (43), resulting in an equation for Z₂. The x-independent part of this equation gives

At z = 0, we have the tangential stress condition and the kinematic conditions. The tangential stress condition for Z₂₀ is obtained by eliminating X₂ from the eqn (48) by taking the horizontal derivative and using the continuity equation, (43). The x-independent part of the resulting equation givesand at z = 0 the x-independent part of the kinematic condition, (47), yieldsThe solution to eqn (71) using the conditions (72) and (73) along with the no displacements conditions at z = −1, gives

Z₂₀ = 0 (74)

Applying this solution to the z-momentum eqn (45), yieldsNoting this, we turn toward the remaining equations at . The potential field is given by

∇²ψ₂ = 0 (76)

subject to the following conditionsThe normal component of the momentum balance at the interface, z = 0, iswhere T^Maxwell_2f consists of the forcing terms which are bilinear, also called (1, 1) terms, because they are products of terms with subscript 1 (cf. Appendix A for the expression of T^Maxwell_2f). The terms associated with [scr D, script letter D]₀ arrive from the expansion of the Maxwell stresses at the second order. We note that all of the forcing terms at this order are bilinear combinations of first order terms, i.e., (1,1) terms with the exception of the boxed term in eqn (77). This implies that the forcing terms are superposition of second harmonics, i.e., cos(2kx) terms and x-independent terms. Their projection on to the eigenspace, i.e., cos(kx) is zero and thus solvability at second order is automatically satisfied. The boxed term in eqn (77) is so identified as it is the sole reason for us to determine [scr A, script letter A]² at the third order. At the second order, ψ₂ is expressed as,where is given byandand where ψ₂₀(z) in eqn (79) is given byThe reason for us to be able to calculate [scr A, script letter A]² at the third order is due to the term ψ^f₂₀. This term arises from the correction of the base state potential at the second order.

Having obtained the solutions to the first and second order equations in terms of [scr A, script letter A], we now turn to the equations at the third order, i.e. to determine [scr A, script letter A]². As noted earlier, the forcing terms in eqn (52)–(60) are bi-linear combinations of second order and first order terms i.e., (2, 1) terms and trilinear combinations of first order terms i.e., (1, 1, 1) terms. From these combinations we can infer that the displacement fields, pressure field and the interface deformation at this order can be expressed as

andAt this order only the cos(kx) part of the displacement, pressure and the interface deformation fields (the terms underlined in eqn (83) and (84)) play a role in determining the amplitude, [scr A, script letter A]. Therefore, we see that the domain equations for the cos(kx) part of the variables areandAgain as before, eliminating X̂₃ from (86) and (87) and using the continuity equation, (85), we getAt z = 0 we have the tangential stress and the kinematic conditions. The tangential stress is modified by eliminating X̂₃ from the cos(kx) part of the eqn (57). This is done by taking the horizontal derivative of the cos(kx) part of the eqn (57) and by using the continuity equation, (52). This givesand at z = 0, the cos(kx) part of the kinematic condition, (56), yieldsSolving eqn (85)–(90) givesWe split into two parts, one that is free of ĥ₃ and the other that is homogeneous in ĥ₃. Thus,where and .

At the third order the potential field is governed by

∇²ψ₃ = 0 (93)

subject to ψ₃ = 0 at z = [script letter H] and the following condition at z = 0The normal component of the momentum balance along the interface at z = 0 giveswhere T^Maxwell_3f consists of the (1, 2) and (1, 1, 1) forcing terms (cf. Appendix A for the expression of T^Maxwell_3f). The boxed term in eqn (95) gives rise to a term in the expression for [scr A, script letter A] which is always negative indicating the subcritical nature of the bifurcation and plays a key role in the high wave number regime, leading to subcritical branching upon instability. It is analogous to a similar term in the Rayleigh–Taylor problem.

At this order, ψ₃ is expressed as

Here, the cos(kx) part of the third order problem alone plays a role in the determination of [scr A, script letter A]² due to the requirement of the solvability condition at the third order (cf. Appendix B for the solution of ). At the third order, solvability requires that the inhomogeneous terms reside in the null space of the homogeneous problem. This yields an expression for [scr A, script letter A]², the details of which may be found in the Appendix C. In the limit of Wi → ∞, [scr A, script letter A]² is given bywith β = [k(cosh(2k[script letter H]) + 2)csch²(k[script letter H])]/[(2 tanh(k[script letter H])(Bo + k²) − 6k²coth(k[script letter H]))/(Bo + k²)]. The above expression for 1/[scr A, script letter A]² is exactly the same that was reported in Dinesh et al.¹ Note that in the limit of Wi → ∞, the bottom fluid mimics the behavior of a Newtonian perfect conductor. We now focus in the expression for 1/[scr A, script letter A]² when the bottom fluid is considered to be a viscoelastic perfect conductor, i.e.,It is noteworthy that the above expression in eqn (98) for 1/[scr A, script letter A]² helps us glean the physics of electrostatic potential acting on a viscoelastic perfect conductor-dielectric fluid pair. A positive sign for 1/[scr A, script letter A]² indicates that the branching behavior for the interface deformation will lead to a supercritical saturation and a negative sign signifies that the branching behavior leads to a subcritical rupture of the interface. Also note that this expression consists of terms containing Wi, these terms play a key role in stabilization of interfacial deformations driven by the electrostatic potential. This stabilization is primarily attributed to the elastic stresses in the perfect conductor that counteract the destabilizing potential. This can be seen from the expression for [scr A, script letter A]² in eqn (99), which is obtained by expanding [scr A, script letter A]² from eqn (98) as a series in the wavenumber k. At low wavenumbers, the [scr O, script letter O](k²) term is dominant. This term is positive and includes the Weissenberg number (Wi), indicating that the interface undergoes supercritical saturation. The presence of Wi shows that elastic stresses contribute to the saturation of interfacial deformations. This term arises from the tangential and normal elastic stresses in eqn (41), (42), (48), (49), (57) and (58).

To further glean the physics of electrostatic instability, we calculate the sign of 1/[scr A, script letter A]² from eqn (98), along the neutral stability curve for Wi = 0.1. In addition, we also calculate the wavenumber beyond which the sign of 1/[scr A, script letter A]² transits from a positive to negative sign. This transition wavenumber is obtained from the expression for 1/[scr A, script letter A]² by equating [scr A, script letter A] ² = 0. The neutral curve and the curve for the transition wavenumber (straight line) are shown in Fig. 3. Note that the region on the neutral stability curve above the straight line which is the transition wavenumber always leads to supercritical saturation of the interface deformation and the region below the straight line leads to subcritical rupture of the interface. The intersection of the neutral curve and the straight line corresponding to the transition wavenumber is a co-dimension 2 point. Observe that in the region for the supercritical behavior, the elastic stresses along with gravity counter act the Maxwell stresses due to the electrostatic potential. These stresses balance out leading to a saturated interface deformation. However, if the fluid were to be Newtonian the Maxwell stresses dominate over gravity and always lead to a subcritical behavior. The region on the neutral curve below the straight line depicts the subcritical branching behavior of the interface. This behavior is attributed to both the Maxwell stresses and the surface tension which ultimately lead to the rupture of the interface. The surface tension term arises from the diminished curvature, i.e., the boxed term in eqn (95) which is of [scr O, script letter O](k³). In short, if the wavenumber k, is chosen such a way that the critical potential lies to the left of the intersection point of the neutral curve and the transition wavenumber, this will lead to a supercritical saturation of the interface. On the contrary, if the critical potential lies to the right of the intersection point, the interfacial deformations lead to a subcritical rupture.

Fig. 3 The critical potential versus wavenumber given by the neutral stability curve and the transition wavenumber is represented by the straight line. The critical potential versus wavenumber to the left of the intersection point lead to a supercritical saturation of the interface and the region to the right of the intersection point represents the subcritical rupture of the interface. The results are presented for Wi = 0.1.

4.2 Comments on the use of the Oldroyd-B fluid model for the perfect conductor: linear stability and weakly nonlinear analysis

A similar stability analysis is performed for the case where the perfect conductor is simulated using the Oldroyd-B fluid model (the details of the linear stability analysis are relegated to Appendix D). For an Oldroyd-B fluid modeled as a perfect conductor, the neutral stability relationship involving the critical dimensionless potential, Bond number, and curvature terms remains identical to the relationship observed when the perfect conductor is modeled as a Newtonian fluid. In addition, the expression for [scr A, script letter A]² in the case of an Oldroyd-B fluid is identical to that obtained when the perfect conductor behaves as a Newtonian fluid (the derivation of the weak nonlinear analysis is provided in Appendix D). This equivalence arises because, in both cases, the perturbed velocities and pressure vanish at orders [scr O, script letter O](ε), , and . As a result, the governing equations remain the same for both Newtonian and Oldroyd-B fluids. Consequently, the bifurcation behavior is also identical, with both fluids exhibiting subcritical rupture of the interface.¹

5 Conclusions

Using weakly nonlinear analysis, it is shown that an electric field imposed on a linear viscoelastic fluid, modeled as a perfect conductor, adjacent to a dielectric air layer in the presence of gravity can lead to either supercritical or subcritical instability of the interface. The analysis is based on a regular perturbation expansion, where the perturbation parameter is defined in terms of the deviation of the applied potential from its critical value.

The results indicate that for low wavenumbers (i.e., wide containers), the interface deforms smoothly and exhibits supercritical saturation. In contrast, at high wavenumbers (i.e., narrow containers), subcritical branching occurs. An analytical expression shows that supercritical behavior arises due to the stabilizing influence of elastic stresses and gravity, which dominate over the destabilizing Maxwell stresses. On the other hand, the subcritical nature of the instability at high wavenumbers is primarily driven by electrostatic forcing and the reduced influence of curvature, which scales as [scr O, script letter O](k³). These effects are sufficient to overcome the stabilizing contributions from gravity and elasticity.

Additionally, when the perfect conductor is modeled as an Oldroyd-B fluid, the analysis reveals that the branching remains subcritical, resulting in interface rupture. This behavior is consistent with that observed in Newtonian fluids.

Author contributions

G. B.: conceptualization, formal analysis, investigation, methodology, software, visualization, writing; D. B.: conceptualization, methodology, investigation, resources, supervision, writing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that support the findings of this study are available within the article.

Appendices

A Maxwell stresses at and at

T ^Maxwell_2f is given by the following expressionT^Maxwell_3f is given by the following expression

B Solution to at

We only need to determine the solution of the cos(kx) part of the third order problem, because this part alone plays a role in the determination of [scr A, script letter A]² at this order. The governing equation for is given byThe corresponding boundary conditions for at z = 0 and at z = [script letter H], areandThe solution to the governing equation for along with its boundary conditions is,

C Solvability condition at

The normal force balance at the third order, yieldsRecall that the cos(kx) part of the normal force balance at the first order is given byMultiplying eqn (C1) with [small zeta, Greek, circumflex]₁ cos(kx) and eqn (2) with [small zeta, Greek, circumflex]₃ cos(kx) and upon subtracting the resulting equations, we get an expression for [scr A, script letter A]².

D Linear stability and weak nonlinear analysis using the Oldroyd-B model for the perfect conductor

We now perform a weakly nonlinear analysis of the electrohydrodynamic instability when the perfect conductor is represented by an Oldroyd-B model. The governing equations are expanded and solved at successive orders: [scr O, script letter O](ε), , and . The ensuing analysis follows the procedure outlined in the work of Dinesh and Narayanan.^1,21 In addition, the solution methodology is outlined in Section 4. For a perfect conductor, the dimensionless governing equations using the Oldroyd-B model are given as follows, beginning with the continuity equation:the x-momentum equationand the z-momentum equationwhere τ_xx, τ_xz, τ_zx, τ_zz are components of the polymeric stress tensor. The xx-component of shear stress is given byThe zz-component of shear stress yieldsThe xz-component of shear stress is represented byHere , and η = μ_p/μ* are the Reynolds number, Deborah number and the viscosity ratio, i.e., the viscosity of the polymer divided by the viscosity of the solvent.

We now perform a linear stability analysis of the governing equations following a similar procedure presented in Section 3. The linearized form of the above equations under neutral conditions, yield,

the x-momentum equationand the z-momentum equationThe xx-component of shear stress isThe zz-component of shear stress is given byThe xz-component of shear stress is represented byIt is noteworthy that the Deborah number, De, is absent in the governing equations corresponding to the shear stress components of the Oldroyd-B model. The Deborah number which takes into the elastic nature of the fluid an important characteristic of the viscoelastic fluid is absent under the neutral conditions, indicating that the elastic nature of the fluid does not alter the critical potential required for an instability. The Deborah number measures how fast stresses relax compared to the flow time scale. In neutral stability, perturbations are marginal with no net growth or decay, so relaxation does not influence the neutral boundary but only affects the rate at which instabilities grow or decay away from it. Consequently, De appears in transient growth rate problems but not in neutral stability curves. It can also be shown that under neutral stability conditions, the perturbed velocity and pressure are zero. Toward this end, we calculate the perturbed velocity in the z-direction by eliminating pressure and the x-component velocity in the governing equations, eqn (C9)–(C14). This yields,The solution for ŵ can be written as

ŵ = c₁e^kz + c₂e^−kz + c₃ze^kz + c₄ze^−kz (C16)

The corresponding boundary conditions, i.e., no slip, no penetration, kinematic and tangential stress in terms of ŵ are given by,Upon solving the governing equations along with the boundary conditions, we get ŵ = 0. Since the perturbed velocities ŵ and û are zero, the perturbed pressure p̂ is also zero. Therefore, upon substituting the perturbed quantities in the linearized normal stress balance yields an expression for the critical potential, i.e.,Note that this expression is similar to the critical potential for a Newtonian fluid (recall eqn (28) in Section 3). This is attributed to the fact that the perturbed velocities and pressure are zero at [scr O, script letter O](ε). Therefore, the Oldroyd-B model for the soft-gel layer results in a critical potential that is independent of the Deborah number. We now focus on the solution of the governing equations at .

D.1 Governing equations and solution. The governing equations at are given by the continuity equation, i.e.,the x-momentum equationand the z-momentum equationwhere (τ_xx)₂, (τ_xz)₂, (τ_zx)₂, (τ_zz)₂ are components of the polymeric stress tensor at . The xx-component of shear stress is given byThe zz-component of shear stress isand the xz-component of shear stress can be written asIt is noteworthy that in the above equations [scr O, script letter O](ε) terms are zero, i.e., the first order velocity components in the viscoelastic fluid layer are zero. In addition, all the time derivative terms are also zero since we are conducting a weak nonlinear analysis around the neutral stability point. Therefore, [scr O, script letter O](ε) terms and the time derivative terms vanish from the momentum and in the polymeric stress tensor equations.

We now calculate the perturbed velocity in the z-direction by eliminating the pressure and the x-component velocity in the momentum equations. This yields,

At this order, velocities, pressure and interface formation are written asThe solution for can be written asThe corresponding boundary conditions, i.e., no slip, no penetration, kinematic and tangential stress in terms of ŵ are given by,Upon solving the governing equations along with the boundary conditions, we get . Note that will play a key role in the calculations at [scr O, script letter O](ε³/6).The solution for can be written as

w₂₀ = c₁z³ + c₂z² + c₃z + c₄ (C31)

The corresponding boundary conditions, i.e., no slip, no penetration, kinematic and tangential stress in terms of w₂₀ are given by,Upon solving the governing equations along with the boundary conditions, we get w₂₀ = 0. It is noteworthy that w₂₀ = 0, will play a key role in the calculations at [scr O, script letter O](ε³/6). The potential field is given by

∇²ψ₂ = 0 (C33)

subject to the following conditionsThe normal component of the momentum balance at the interface, z = 0, is(cf. Appendix A for the expression of T^Maxwell_2f) where is given byandand where ψ₂₀(z) is given byA detailed solution for the potential is provided in Dinesh et al.¹ We now proceed to the solution of the governing equations at .

D.2 Governing equations and solution. The governing equations at are given bythe x-momentum equationand the z-momentum equationwhere (τ_xx)₃, (τ_xz)₃, (τ_zx)₃, (τ_zz)₃ are components of the polymeric stress tensor at . The xx-component of shear stress is expressed asThe zz-component of shear stress can be written asand the xz-component of shear stress is represented byIt is noteworthy that in the above equations [scr O, script letter O](ε) and terms are zero, i.e., the first and second order velocity components and pressure in the viscoelastic fluid layer are zero. Therefore, we do not see the product of these [scr O, script letter O](ε) terms and product of [scr O, script letter O](ε) term with in the momentum equations and the polymeric stress tensor equations. Toward this end, we calculate the perturbed velocity in the z-direction by eliminating pressure and the x-component velocity in the momentum equations. This yields,At this order, velocities, pressure and interface deformation are written asThe solution for ŵ₃ can be written as

ŵ₃ = c₁e^kz + c₂e^−kz + c₃ze^kz + c₄ze^−kz (C48)

The corresponding boundary conditions, i.e., no slip, no penetration, kinematic and tangential stress in terms of ŵ are given by,Upon solving the governing equations along with the boundary conditions, we get ŵ₃ = 0. At the third order the potential field is governed by

∇²ψ₃ = 0 (C50)

subject to ψ₃ = 0 at z = [script letter H] and the following condition at z = 0The normal component of the momentum balance along the interface at z = 0 giveswhere T^Maxwell_3f consists of the (1, 2) and (1, 1, 1) forcing terms (cf. Appendix A for the expression of T^Maxwell_3f). Therefore, at this order, ψ₃ and ζ₃ are expressed asHere, the cos(kx) part of the third order problem alone plays a role in the determination of [scr A, script letter A]² due to the requirement of the solvability condition at the third order. The solvability condition yields an expression for [scr A, script letter A]², i.e.,with β = [k(cosh(2k[script letter H]) + 2)csch²(k[script letter H])]/[(2 tanh(k[script letter H])([scr B, script letter B] + k²) − 6k² coth(k[script letter H]))/([scr B, script letter B] + k²)]. The expression for [scr A, script letter A]² in the case of an Oldroyd-B fluid is identical to that obtained when the perfect conductor behaves as a Newtonian fluid. This equivalence arises because, in both cases, the perturbed velocities and pressure vanish at orders [scr O, script letter O](ε), , and . As a result, the governing equations remain the same for both Newtonian and Oldroyd-B fluids. Consequently, the bifurcation behavior is also identical, with both fluids exhibiting subcritical rupture of the interface.¹ The Deborah number characterizes the relative importance of elastic relaxation compared to the characteristic time scale of the flow, indicating whether the material responds more like a solid (large De) or more like a fluid (small De). In the case of neutral stability, perturbations are marginal and neither grow nor decay. At this threshold, stress relaxation does not contribute to shifting the onset of instability, since the system is balanced precisely between stable and unstable behavior. Instead, the Deborah number influences the dynamics away from neutrality by determining how rapidly instabilities amplify or decay. Consequently, De plays a central role in governing transient growth rates of perturbations, but it does not explicitly appear in the neutral stability curves that define the onset of instability.

Acknowledgements

The authors gratefully acknowledge funding from Ministry of Higher Education (MHRD), Anusandhan National Research Foundation (ANRF) via the grant numbers ITS/2024/005066 and ANRF/ECRG/2024/001265/ENS.

References

1. B. Dinesh and R. Narayanan, Phys. Rev. Fluids, 2021, 6, 054001.
2. G. I. Taylor, Proc. R. Soc. London, Ser. A, 1966, 291, 159–166.
3. J. R. Melcher and G. I. Taylor, Annu. Rev. Fluid Mech., 1969, 1, 111–146.
4. E. SchaÈffer, T. Thurn-Albrecht, T. P. Russell and U. Steiner, Nature, 2000, 403, 874–877.
5. V. Shankar and A. Sharma, J. Colloid Interface Sci., 2004, 274, 294–308.
6. R. M. Thaokar and V. Kumaran, Phys. Fluids, 2005, 17, 084104.
7. L. Espn, A. Corbett and S. Kumar, J. Non-Newtonian Fluid Mech., 2013, 196, 102–111.
8. K. Mondal, P. Kumar and D. Bandyopadhyay, J. Chem. Phys., 2013, 138, 1–15.
9. J.-U. Park, M. Hardy, S. J. Kang, K. Barton, K. Adair, D. K. Mukhopadhyay, C. Y. Lee, M. S. Strano, A. G. Alleyne and J. G. Georgiadis, et al. , Nat. Mater., 2007, 6, 782–789.
10. D. Riccobelli and P. Ciarletta, Philos. Trans. R. Soc., A, 2017, 375, 20160421.
11. J. Marthelot, E. F. Strong, P. M. Reis and P. T. Brun, Nat. Commun., 2018, 9, 1–7.
12. D. T. Papageorgiou, Annu. Rev. Fluid Mech., 2019, 51, 155–187.
13. S. Y. Chou and L. Zhuang, J. Vac. Sci. Technol., B: Microelectron. Process. Phenom., 1999, 17, 3197–3202.
14. S. A. Roberts and S. Kumar, J. Fluid Mech., 2009, 631, 255.
15. P. Gambhire and R. M. Thaokar, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 036301.
16. K. Ward, S. Matsumoto and R. Narayanan, J. Fluid Mech., 2019, 862, 696–731.
17. D. Koulova-Nenova, J. Electrost., 1997, 40, 185–190.
18. P. Colinet, J. C. Legros and M. G. Velarde, Nonlinear dynamics of surface-tension-driven instabilities, Wiley-VCH, Berlin, 2001.
19. L. Wu and S. Y. Chou, J. Non-Newtonian Fluid Mech., 2005, 125, 91–99.
20. G. Tomar, V. Shankar, A. Sharma and G. Biswas, J. Non-Newtonian Fluid Mech., 2007, 143, 120–130.
21. B. Dinesh and R. Narayanan, J. Fluid Mech., 2021, 915, 1–27.
22. L. D. Landau and E. M. Lifshitz, Theory of Elasticity, Pergamon Press, 1989.
23. V. Shankar and V. Kumaran, J. Fluid Mech., 2000, 407, 291–314.
24. B. Dinesh and S. Pushpavanam, Phys. Rev. E, 2017, 96, 013119.
25. R. Patne, D. Giribabu and V. Shankar, J. Fluid Mech., 2017, 827, 31–66.
26. V. Kumaran, J. Fluid Mech., 1995, 294, 259–281.
27. V. Kumaran, J. Fluid Mech., 1995, 302, 117–139.
28. V. Kumaran, J. Fluid Mech., 1998, 357, 123–140.
29. B. Dinesh and S. Pushpavanam, Phys. Fluids, 2018, 30, 104103.
30. P. Howell, G. Kozyreff and J. Ockendon, Applied solid mechanics, Cambridge University Press, 2009.
31. P. Ramkarn and V. Shankar, J. Fluid Mech., 2018, 860, 837–885.
32. J. Snoeijer, A. Pandey, M. Herrada and J. Eggers, Proc. R. Soc. A, 2020, 476, 20200419.
33. D. A. Saville, Annu. Rev. Fluid Mech., 1997, 29, 27–64.
34. S. Mora, T. Phou, J.-M. Fromental and Y. Pomeau, Phys. Rev. Lett., 2014, 113, 178301.
35. H. W. Müller and W. Zimmermann, EPL, 1999, 45, 169.
36. L. E. Johns and R. Narayanan, Interfacial instability, Springer Science & Business Media, 2002, pp. 255–304.

---