Active elastocapillarity in soft solids with negative surface tension.

Active solids consume energy to allow for actuation, shape change, and wave propagation not possible in equilibrium. Whereas active interfaces have been realized across many experimental systems, control of three-dimensional (3D) bulk materials remains a challenge. Here, we develop continuum theory and microscopic simulations that describe a 3D soft solid whose boundary experiences active surface stresses. The competition between active boundary and elastic bulk yields a broad range of previously unexplored phenomena, which are demonstrations of so-called active elastocapillarity. In contrast to thin shells and vesicles, we discover that bulk 3D elasticity controls snap-through transitions between different anisotropic shapes. These transitions meet at a critical point, allowing a universal classification via Landau theory. In addition, the active surface modifies elastic wave propagation to allow zero, or even negative, group velocities. These phenomena offer robust principles for programming shape change and functionality into active solids, from robotic metamaterials down to shape-shifting nanoparticles.

INTRODUCTION

Embedding stress-generating active elements into a passive solid powers functionality inaccessible in thermal equilibrium. These active metamaterials occupy the space between materials and machines, giving rise to exotic phenomena from actuation and shape change (–) to overdamped wave propagation (–) and nonreciprocal interactions (–). For two-dimensional (2D) active interfaces, powerful design principles exist for the distribution and control of stress-generating elements to achieve a target behavior (, , –). A key challenge is to develop these principles for the control of bulk 3D solids. Realizations of far-from-equilibrium solids range from macroscopic mechanical materials (, ) and hydrogels (–) down to the microscale (, , , ). The challenge of spanning these systems and length scales requires principles based on a continuum approach.

In thermal equilibrium, the shape and structure of a soft solid is determined not only by 3D bulk elasticity but also by surface stresses on its 2D boundary. This competition, termed (passive) elastocapillarity (, ), has been used to stiffen composites (), self-assemble micro-objects (, ), and drive the coiling of nanoparticle helices (). These phenomena all originate in the minimization of surface area due to the isotropic surface stress tensor , where δ is the Kronecker delta. Because passive elastocapillary solids are in equilibrium, their surface tension γ is constrained to be positive.

In this work, we show that elastocapillarity offers a distinct set of design principles when considering 3D soft solids driven out of equilibrium. The interplay between bulk and surface stresses translates to this far-from-equilibrium context, but now, each stress tensor may acquire additional, active terms (, ). We refer to the resulting phenomenology as active elastocapillarity. This scenario stands in contrast to the buckling of elastic shells (, –) or morphing of thin programmable sheets (, , , ). Rather than a thin layer, we consider a fully 3D solid (Fig. 1A), which is sufficiently soft for active surface stresses to strain the entire bulk. We focus on a minimal change to the surface stress tensor due to activity. This is the addition of an isotropic dilational stress , where γ < 0. Sufficiently strong dilations will overpower passive tension γ, giving an overall surface tension γ, which is effectively negative, γ ≡ γa + γp < 0 (, , ).

Fig. 1.

Active elastocapillary phenomena in soft solids with negative surface tension.

(A) At the boundary of a bulk 3D soft solid, an active layer (blue) exerts dilational surface stresses (red arrows), which compete against bulk elasticity (orange). Left: A cutaway view of the entire solid. Right: A zoom in near the boundary. (B and C) Mechanisms leading to active surface stresses, from the nanoscale to the macroscale. (B) Insertion of complementary polymers into a nanoparticle surface causes steric crowding (). (C) An active metamaterial in which mechanical actuators are embedded in the surface of an elastic medium (, , ). (D to G) Continuum phenomenology. (D) Switching between spheres and complex shapes of higher surface area. (E) Negative group velocity of elastocapillary Rayleigh waves. (F) Sharpening of corners and edges, with a cone of tip angle θ deforming into a cusp of the form ∣Z(R)∣∼ θ−1R1/4. (G) Promotion of crack propagation.

Intuitively, positive surface tension rounds any 3D object into a sphere. For negative surface tension, does a unique favored shape exist? Developing a continuum theory of active elastocapillarity, here, we show that final shape is not determined by∣γ∣alone. Instead, the shape can be selected by varying either surface stresses or bulk elastic moduli. As active driving increases, we find that distinct geometries discontinuously snap between one another, as realized in a simple particle-based numerical model. As well as shape, active surface stresses also control dynamic phenomena. For example, we find that negative surface tension softens Rayleigh wave propagation, leading to zero (or even negative) group velocity. Together, our results form a toolkit for programming the functionality of 3D active solids.

Experimental realizations

In Fig. 1 (B and C), we give examples of experimental mechanisms for how active surface stresses can be designed. At the microscale, Fig. 1B shows a nanoparticle made of a diblock copolymer with a hydrophilic head (blue) and a hydrophobic tail (orange). The long tails form a melt, leading to bulk elasticity inside the nanoparticle (). Insertion of a complementary polymer (blue head and green tail) into the nanoparticle surface drives area growth. As an example, Hua et al. () achieved this insertion using binding between paired DNA nucleobases (orange to green). To accommodate this growth, the nanoparticles deform away from their equilibrium spherical shape.

On the macroscale, a simpler realization is to embed dilational motors into the surface of a solid rubber ball (or another object), as shown in Fig. 1C. Alternatively, macroscale surface dilation can be designed using self-folding origami sheets or other deployable structures that spontaneously grow their effective surface area as they unfold (, ). The incompressibility of the underlying rubber ball prohibits isotropic expansion at fixed spherical shape. Instead the ball shears, with the coupling between rubber shear modulus and dilational surface stresses in the active sheet determining the final shape of the composite object.

We focus on physical realizations using active elements such as independent motors. Similar to all examples of active matter, these motors locally consume energy to exert mechanical stresses. When these decentralized stresses couple together, shape change can become an emergent phenomenon via spontaneous symmetry breaking. An alternative mechanism for generating negative surface tension commonly used to make wrinkling patterns (, , , –) is to globally prestress a planar sheet, which is then glued onto a flat substrate (). By contrast, active elements can lead to a negative surface tension even when embedded into the surface of an arbitrarily shaped object. These dilational elements are readily found across many length scales (, ). Here, we take a continuum approach that is independent of the microscopic origin of negative surface tension or system scale. For example, in Fig. 1B, control over the continuum-level surface tension γ can be implemented by varying the concentration of inserted polymer. In Fig. 1C, such control can be achieved by varying the forces exerted by the macroscopic motors.

RESULTS

Figure 1 (D to G) illustrates a selection of the active elastocapillary phenomena that arise because of spontaneous growth of surface area. Below, we focus on exact solutions for deformations of a sphere (Fig. 1D) and surface waves and instabilities (Fig. 1E). First, we give an example that typifies the phenomenology of active elastocapillarity: the sharpening of a cone to a cusp (Fig. 1F).

Passive elastocapillarity smooths out a vertical cone of tip angle θ (–). By contrast, active elastocapillarity will sharpen this feature to a power law cusp. The solid’s boundary is then described by the curve ∣Z(R)∣∼θ−1R1/4, with the deformed conical height Z converging sharply to the origin as a function of its deformed radius R. To derive this power law, we consider the elongation of an incompressible elastic cylinder, of elastic modulus μ and undeformed radius ρ, under the action of negative surface tension on its curved boundary. The stretch factor λ > 1 is found by balancing the elastic deformation energy μρ2λ2 against the surface energy to give λ ∼ (∣γ∣/μρ)2/3. The deformed height is then Z = λz, with a corresponding radial contraction of . We then take the cone to be a stack of infinitesimally thin cylinders of progressively decreasing radius ρ ∼ θz. Each cylinder experiences a z-dependent elongation λ(z) ∼ (θz)−2/3. Integrating these extensions and expressing the result in terms of the deformed radius, R ∼ θzλ(z)−1/2 ∼ (θz)4/3, yields the R1/4 power law. In the Supplementary Materials, we describe this theoretical approach in detail and provide in fig. S1 a numerical demonstration of the cone-to-cusp transition, as well as confirmation of the R1/4 power law. Sharpening of corners corresponds to the case of a small conical angle θ. Taking a large θ instead corresponds to crack proliferation through stress concentration, shown in Fig. 1G.

Continuum theory

We now proceed to the general framework of our continuum description. Active elastocapillarity is defined by two intrinsic length scales. The first, so-called elastocapillary length lγ ≡ ∣γ∣/μ, measures the ratio of effective surface tension ∣γ∣ to shear elastic modulus μ of the solid. Intuitively, at length scales larger than the elastocapillary length lγ, 3D elasticity stabilizes the surface. What happens at ever smaller scales? In passive elastocapillarity, stability is provided by positive surface tension. By contrast, within active elastocapillarity, the destabilizing contribution of negative surface tension γ must be regularized by higher-gradient stresses. We focus here on surface effects and introduce a second bendoelastic length lκ ≡ (κ/μ)1/3. This length scale can arise, for example, from a surface with bending modulus κ () or from the length dependence of active stresses γa. A distinct stabilization mechanism is bulk dispersion. In the Supplementary Materials, we consider the effects of a shear modulus μ(q) = μ0 + μ1q2, dependent now on the planar wave number q, which instead leads to a stabilizing length scale . We contrast this effect with viscous dissipation, from which no such stabilization is possible.

An object’s shape results from the competition between bulk elasticity and boundary conditions containing active surface stresses. We solve the equations of linear elastodynamics with a stress-matching condition at the surface (–)for slow variations in initial curvature H, where σ is the component of the 3D elastic stress tensor normal to the surface, ∇2 is the covariant (surface) Laplacian, and δH is the variation of mean curvature (see Fig. 2A). In “Linear instability and shape change in an active elastocapillary droplet” section, we discuss the derivation of Eq. 1 and details of our solutions. Notably, our approach accounts for 3D elastic coupling between active surface elements. This allows us to describe deformations that are inaccessible when treating bulk elasticity as an external 2D potential coupled to a surface energy [cf. (, –, –)]. We contrast these approaches further in the Supplementary Materials.

Fig. 2.

Shape instability and wave propagation within active elastocapillarity.

(A) Schematic of an elastic solid with shear modulus μ and surface bending rigidity κ deformed by a negative surface tension −∣γ∣, which gives excess mean curvature H → H + δH. (B) Phase diagram for the stability of an active elastocapillary sphere of radius R, controlled by rescaled surface tension and bending modulus [≡κ/(μR3)]. The border between the white and the colored regions corresponds to activity at which the sphere goes unstable. Color indicates the dominant azimuthal mode number l, with insets showing unstable modes at l = 2,3, and 4. (C) Dispersion of spherical oscillations corresponding to points marked by square, triangle, and circle in (B). (D and E) The limit l → ∞ describes a half-space, with 1D phase diagram (E) controlled by rescaled bending modulus . Colored curves (red, orange, black, and blue) and regions (green and gray) of (D) map to corresponding points and regions of (E). At , the half-space destabilizes at rescaled wave number [(E), green region]. For , , we find negative group velocity surface waves [(E), gray region]. (F) Surface wave dispersions corresponding to colored lines in (D) and (E).

For any shape of size R, the solutions that we find must be characterized by the length scale triplet (lγ, lκ, R). We then define two independent dimensionless ratios as the surface tension and bending modulus rescaled by the size R: and . When size R is too large to be relevant (such as for an infinite half-space), the solutions depend only on the quantity . Here, describes the ratio of stabilizing elasticity (μ and κ) to destabilizing activity ∣γ∣. We conclude that under overall rescaling, the phenomenology remains unchanged and that continuum elastocapillarity remains valid across any scale, from nanoparticles to macroscopic metamaterials.

Tunable instabilities and shape selection

Spheres minimize area at fixed volume, and positive surface tension drives every initial shape toward that of a sphere. By contrast, negative surface tension drives transitions away from a sphere into a variety of shapes. Our exact results demonstrate how to select between these shapes using the elasticity of the underlying solid. The phase diagram in Fig. 2B shows that for low active driving, is small and spherical shapes are stable. As activity increases, spheres spontaneously destabilize. The threshold activity for instability, encoded in , and the azimuthal wave number l of the dominant unstable mode are both determined by the balance of surface and 3D moduli (∼κ/μ).

When the bending modulus is small compared to 3D elasticity (small ), we find an instability at wavelength ∼lκ, allowing for control over fine structure and surface texture. The threshold and wave number of this instability map to those of the wrinkling instability of a stiff, prestressed layer bonded onto a substrate (, , , –, ). To make this mapping, we identify the bendoelastic length lκ with the physical thickness of the stiff layer and the negative surface tension γ with a uniform dilational stress within the layer, as described in the Supplementary Materials.

However, as bulk elasticity weakens (large ), both the wavelength and penetration depth of these wrinkles become comparable to the system size. In this limit, we cannot describe the instability as that of a nearly planar surface, and our approach via fully 3D elasticity becomes essential (see the Supplementary Materials). Here, we find destabilization at the level of global shape change (Fig. 2B, right-hand side). In particular, at the largest , we find shapes with uniaxial anisotropy, which are inaccessible, for example, via instabilities at fixed surface area (). These unstable modes originate in the dispersion relations shown in Fig. 2C; here, l is the azimuthal wave number, and frequency is rescaled by the sphere radius R and transverse speed of sound (with ρ as the bulk density). Figure 2C shows dispersions corresponding to points marked by square, triangle, and circle in Fig. 2B, just above the onset of instability. In each case, as threshold negative surface tension is crossed, a single mode is driven unstable (, dashed line). This instability gives a scale-free tool for designing global shape change in bulk 3D solids, distinct from the buckling of thin elastic shells ().

For a flat object, energy injection will soften surface waves and drive wrinkling instabilities [see Fig. 2 (D and E)]. Figure 2D shows a zoom-in of the l → ∞ (equivalently R → ∞) limit of Fig. 2B. In this limit, the phase diagram collapses to the 1D diagram shown in Fig. 2E, controlled solely by the rescaled bending modulus : Each curve in Fig. 2D maps to a value of in Fig. 2E (colored bars). At a threshold , the half-space destabilizes (Fig. 2E, green region), which, in unscaled variables, occurs when active driving ∣γ*∣ = 3(κμ2)1/3. The length scale for this instability is characterized by a rescaled planar wave number (Fig. 2E). At instability, the critical wave number or in unscaled units .

Even below instability, a vestige of surface activity can be measured via the negative group velocity of surface elastocapillary waves (Fig. 2E, gray region). In Fig. 2F, we show the planar dispersion relation in which the frequency ω is now rescaled by both the elastocapillary length lγ and the speed of sound c. For weak active driving (large ), the half-space is stable and plane waves stiffen at high wave numbers. As active driving increases ( decreases), energy injection softens high wave numbers, leading first to negative group velocity and then to full-blown surface instability. Intuitively, this behavior stems from an effective shift of the shear modulus by negative surface tension, μ → μ − ∣γ∣ q (see the “Waves and instabilities in an active elastocapillary half-space” section). This rescaling causes μ, restoring elastic forces, and phase velocities all to vanish on a scale set by the elastocapillary length lγ.

The instability thresholds for wave number q* and active driving ∣γ*∣ can both be tuned using the surface modulus κ and 3D shear modulus μ. In other words, by selecting the material parameters of the passive solid, we can select the first mode that goes unstable once activity is turned on.

Nonlinear elasticity and universality

Linear analysis can select only mode number, not mode amplitude. This amplitude is typically selected by nonlinear terms, which break symmetry. For example, geometric curvature of a surface will break the up-down symmetry of normal displacements, leading to the selection of dimples or labyrinthine structure in wrinkling (–). Here, we explore a distinct symmetry breaking mechanism: bulk material nonlinearity. Nonlinear elastic response will break the symmetry between positive and negative amplitude deformations present in our linear analysis. Thereby, nonlinearity selects global shape.

We focus on the l = 2 mode of uniaxial deformations. The amplitude of this uniaxial strain can be approximated by a homogeneous deformation with principal stretch factor λ (see Fig. 3A). Stretch factor λ > 1 corresponds to a worm-like (prolate) shape, and λ < 1 corresponds to a pancake-like (oblate) one. A general nonlinear elasticity introduces an infinite set of materials parameters, making inaccessible an exact solution such as the one that we obtained in the linear regime. In keeping with the minimal approach taken thus far, we first consider a simple model of nonlinear elasticity that accounts for the effects of material nonlinearity: the Mooney-Rivlin model, often used for rubbers and polymer gels such as those shown in Fig. 1 (B and C) (). We will return to the importance of this choice below. An effective energy for this far-from-equilibrium solid is given by

Fig. 3.

Nonlinear active elastocapillarity selects mode amplitude.

(A) Homogeneous deformations are described by stretch factor λ, with λ > 1 giving a worm and λ < 1 giving a pancake. (B) Worm/pancake phase diagram, in activity and material nonlinearity . Spheres destabilize as increases, but the resulting shape (worm/pancake) depends on . At fixed , worms suffer a second snap-through transition at large . These discontinuous transitions meet at a critical point. (C) Landscapes for effective energy describing discontinuous transitions along the arrows marked in (B). (D and E) Simulations of minimal ball-spring model demonstrate continuum theory predictions (see movie S1). By tuning spring-level nonlinearity, we realize worms (top), pancakes, and snap-through (bottom). (E) Theory-simulation comparison for dependence of λ on active driving . For numerical data, orange squares and blue triangles correspond to data shown in (D). Theoretical fit shows all (meta)-stable minima for the Mooney-Rivlin continuum theory, Eq. 2, with parameters matched to the orange squares (). (F) An anisotropic shape, with initial stretch factor λ0, splits the critical point from part (B). For both worm (λ0 > 1, top right) and pancake (λ0 < 1, bottom right), initial anisotropy either grows smoothly or snaps through, depending on the path in parameter space. A cut in space through the critical point (bottom left) gives an Ising-like transition described by Landau theory, Eq. 3.

For the active surface contribution , we take the area of a uniaxial ellipsoid (Fig. 3A). This is balanced by a Helfrich bending energy () and equilibrium bulk elasticity, composed of the neo-Hookean contribution and the Mooney-Rivlin term . We give expressions for the surface and bending energies in terms of λ in the “Nonlinear theory” section. Crucially, this model includes a single parameter to continuously tune material nonlinearity, which ranges between the neo-Hookean limit of only geometric nonlinearity, , and maximal nonlinearity, .

Minimizing Eq. 2 yields the phase diagram shown in Fig. 3B. The discontinuous shape transitions indicated by solid lines correspond to bistable configurations in the effective energy (Fig. 3C). Increasing active driving destabilizes spheres, but now, the elastic nonlinearity controls whether the resulting shape is a worm or a pancake. In the neo-Hookean limit , the preferred shape is a compressed pancake. Intuitively, this corresponds to maximizing surface area without any elastic effects (Fig. 3C, left). For sufficiently large , the preferred shape is instead an elongated worm (Fig. 3C, middle), reflecting the bias toward uniaxial elongation over compression encoded in the Mooney-Rivlin theory (). These worms can undergo a second snap-through transition (Fig. 3C, right), morphing to pancakes as active driving is further increased or as the nonlinearity is tuned.

In Fig. 3 (D and E), we demonstrate these continuum predictions in a microscopic model of active elastocapillarity, as shown in movie S1. We simulate the deformation of a spherical mesh of bulk nonlinear springs, coupled to surface springs exerting active stresses, as shown in Fig. 3D (see the “Numerics” section for details). At a critical active stress, the meshed solid destabilizes, with the resulting shape tunable via spring-level nonlinearity. By tuning the springs such that crosses the critical value , we realize worms (Fig. 3D, top), pancakes, and snap-through transitions (Fig. 3D, bottom), in agreement with analytical predictions. In Fig. 3E, we compare the continuum theory Eq. 2 with numerical data for the stretch factor λ. Significantly, the quantitative theory-simulation agreement near the transition hints at universality.

The three transition lines in Fig. 3B meet at a critical point. This critical point is a direct consequence of the symmetry of the initial shape. In contrast to a sphere, an initial uniaxial anisotropy λ0 causes the critical point to split (Fig. 3F). For example, the phase diagram shown in Fig. 3F (top right) demonstrates how an initially elongated shape with λ0 > 1 can smoothly extend as active driving is increased. However, to compress into a pancake-like shape, the object must still cross a snap-through transition. The converse is true for an initially compressed shape (Fig. 3F, bottom right). In activity-anisotropy space, a cut through the critical point reveals an Ising-like transition (Fig. 3F, bottom left), with playing the role of temperature and λ0 playing the role of an external field. This observation motivates a universal characterization of shape transitions near the critical point using Landau theory.

In the above discussion, we introduced the Mooney-Rivlin model as a minimal choice. Other nonlinear elasticities will yield quantitatively different phase diagrams. However, near the critical point, the complete behavior of the active solid can be understood using symmetry-based arguments. For any initial shape, the Landau expansion guarantees that the effective free energy has the formwhere the linearized strain ϵ(=λ − 1) plays the role of an order parameter, the control parameter probes the distance to linear instability, within the Mooney-Rivlin model) is the lowest-order nonlinear term, and u > 0 guarantees stability. The linear term h ∼ λ0 − 1 captures the effects of either shape anisotropy or external uniaxial stresses and is absent for cubes, spheres, and other spherical tops (i.e., shapes with an isotropic moment-of-inertia tensor). As a result, a critical point is generically expected for these symmetric shapes, with three weakly discontinuous transitions emanating from it, as in Fig. 3B. Although Landau theory breaks down at higher strains, this critical point controls the entire phase diagram shape.

We derive expressions for parameters r, w, u, and h within the Mooney-Rivlin model in the “Nonlinear theory” section. However, we emphasize that the form of Eq. 3 is fully constrained by symmetry. Hence, we expect the qualitative structure of our results to hold for a general nonlinear elasticity; the Landau theory presents a universal classification across all elastocapillary materials and shapes. The essential feature is simply that nonlinear effects appear in w, as is generically the case. In the Supplementary Materials, we explore the phase diagram of the Gent model of rubber elasticity (, ), a singular example where nonlinearities only appear at quartic order.

DISCUSSION

Active elastocapillarity couples field theories of different dimensionalities toward new materials design principles. Here, we have adopted a continuum approach applicable across systems, using linear elasticity to describe shape selection and wave propagation analytically. We have found quantitative agreement between continuum nonlinear elasticity and particle-based numerics in predicting both the final shape of active elastocapillary solids and snap-through transitions between them. Landau theory explains this agreement in terms of universal behavior about a critical point, allowing a classification based solely on the initial symmetry of the solid.

We have focused on the minimal case of a passive elastic solid, coupled to active stresses in the form of an effectively negative surface tension ϒ = −∣γ∣δ. The concept of active elastocapillarity may include a broader range of phenomena in which either bulk or surface stresses contain active components. While our attention has been on synthetic realizations, some of these ideas may also be relevant across living systems where surface stresses affect shape. Examples of biological phenomena with surface growth and (2D or 3D) elasticity include cellular symmetry breaking (), buckling of thin actin shells (), cellular layers (), tissues (, , ), and modified wetting via differential growth (). In contrast to the complexity of these systems, here, we have highlighted how minimal ingredients such as isotropic surface stresses and bulk elastic nonlinearity can already be used to design unexpected 3D metamaterial functionality.

The instabilities that we have uncovered, from snap-through to smooth deformations, suggest active elastocapillarity as a portable mechanism to achieve complex reconfigurable shapes. At the macroscale, we envision the design of soft robotic arms composed of anelastic backbone covered in simple actuators (Fig. 1C). Scaling active elastocapillarity down to soft nanoparticles (Fig. 1B), for which no reliable shape-control mechanism exists, may prove useful for applications ranging from drug delivery to self-assembly of photonic crystals.

MATERIALS AND METHODS

Linear instability and shape change in an active elastocapillary droplet

In this section, we derive the dispersion relation and linear instabilities of an active elastocapillary sphere. Recently, Tamim and Bostwick () have given an analysis of the vibrations of the passive case, extending classical results for the purely elastic () and capillary limits (). Here, we instead consider the active case in which the surface terms consist of both a negative surface tension γ and bending modulus κ. The approach is to take a bulk ansatz satisfying the equations of linear elastodynamics () and impose a matching boundary condition between bulk and surface stresses. This boundary condition leads to a solvability criterion, whose solution gives the dispersion and regimes of linear instability.

We first derive the boundary condition, balancing active stresses with restoring elasticity. Given a (2D) surface with surface tension γ and bending rigidity κ, a variation of the Helfrich surface free energygives the shape equation for vesicles ()where P is (minus) the normal component of the bulk stress tensor, H is the mean curvature of the solid’s boundary, K is the Gaussian curvature, and c0 allows for a preferred nonzero mean curvature. We now expand H, K, and P to first order about a spherical shape, considering a normal perturbation ψ, where is the outward unit normal

Here, δH and δK are given by (, )

Substituting Eqs. 6 and 7 into the shape equation Eq. 5 yields the normal component of the stress matching conditionwith the tangential component στ = 0, where τ denotes the surface tangent. For an expansion in terms of spherical harmonics, we take , with as the associated Legendre polynomial and N as a normalization factor (). Then, ∇2ψ = − l(l + 1)ψ and Eq. 8 simplifies to

From Eqs. 8 and 9, we see the effect of the bending modulus and spontaneous curvature is to shift γ as

Two natural special cases of this result are c0 = 0 (no spontaneous curvature) and c0 = − 1/R (spontaneous curvature matching the initial mean curvature H0). Here, we focus on the case c0 = 0, for which Eq. 8 simplifies to Eq. 1. However, note that only the last term in Eq. 10 depends on l and the effect of c0 can be completely reabsorbed into the effective surface tension γ.

With the mapping of Eq. 10, the problem reduces to that considered in (). The bulk stress tensor σ is given by the general solution to the linear elastodynamic equations in spherical coordinates (). Substituting this solution into Eq. 9 gives the solvability condition, which we invert numerically to obtain the dispersion. Below, we state this solvability condition for the limit of an incompressible material, in terms of dimensionless variables coming from the sphere radius R and associated timescale

In terms of these dimensionless quantities, the dispersion is given by the solution ofwhere j is a spherical Bessel function of the lth kind.

Figure 2 (B and C) was found by solving Eq. 12 numerically. Note that the dispersion Eq. 12 has an infinite number of branches, corresponding to the roots of j, spheroidal modes ω(s, l) are indexed by a radial “quantum number” s and polar wave number l, being degenerate with respect to the azimuthal wave number m. Only the s = 1 branch couples to the instability described here, and it is this branch that is shown in Fig. 2 (B and C).

Waves and instabilities in an active elastocapillary half-space

In this section, we derive the spectrum of the linearized equations of motion for an active elastocapillary half-space. These results follow from the l → ∞ limit of the “Linear instability and shape change in an active elastocapillary droplet” section, in particular Eq. 12. However, an independent derivation in the planar case has the virtue of being much simpler than the spherical case, and we extend it to study the effects of viscosity and bulk dispersion in the Supplementary Materials. Passive elastocapillary waves have been studied from the perspective of a viscous fluid () or an elastic solid. Here, we take the elastic solids perspective () and consider the active case in which a negative surface tension γ is regularized in the high wave number limit by a bending modulus κ. Our approach applies equally to 2D or 3D materials.

Before giving a detailed argument, basic scaling considerations capture the main phenomenology. Consider a slab of material of undeformed surface area A with a deformed surface height h(x). The energy stored in surface deformations is , and the bulk energy E ∼ μh′(x)2Al, where l is the depth that surface deformations penetrate into the bulk. Assuming l ∼ 1/q, the total energy per unit volume is then , and we see that γ acts as a q-dependent shift to the shear modulus, . As q increases, for γ < 0, μ(q) softens, with higher wave numbers feeling progressively weaker elastic restoring forces. At q ∼ 1/lγ, restoring elasticity vanishes entirely, with κ regularizing high wave numbers. We thus expect the threshold wave number for instability to scale as 1/lγ.

We now give a detailed derivation of the spectrum. We consider a half-space z ≤ 0. The equations of linear elastodynamics in the bulk material are ()where ρ is the density and u is the displacement. For an isotropic material, the stress tensor , where is the linearized strain and d is the spatial dimension. Equation 13 supports longitudinal and transverse waves, of wave vector q and frequency ω, propagating along x and decaying as z → − ∞ (, )

Here, and are the inverse penetration depths along z, with the longitudinal modulus . A general displacement u = ue + ue is written as

We now take the bulk ansatz Eq. 21 and substitute it into a stress matching boundary condition at z = 0. The boundary has surface tension γ, bending modulus κ, and effective free energywhere h(x) is the height of the free surface above z = 0. Matching the z component of bulk displacement u to height h yields the stress matching condition

Substituting Eq. 15 into Eq. 17 gives

At this point, we take the incompressible limit: as B, M → ∞, α → q, and Eq. 18 simplifies to

The dispersion is obtained from the solvability condition that the determinant of Eq. 19 must vanish

Equation 20 reduces to the passive elastocapillary dispersion (in the incompressible limit) for γ > 0 and κ → 0 (), and the standard Rayleigh wave dispersion for κ, γ → 0 (). Intrinsic length and time scales are given bywhich we use to define nondimensionalized variables

In dimensionless form, Eq. 20 is then

Here, we focus on the case of negative surface tension, sgn(γ) = −1. The phase diagrams shown in Fig. 2 (D and E) and the dispersions shown in Fig. 2F are found by solving Eq. 23 numerically. However, the threshold at which instability occurs can be found analytically: letting in Eq. 23 gives the condition

We require the cubic in Eq. 24 to have degenerate roots, as must be the case at instability. This gives values for the wave number and dimensionless bending modulus at which instability sets in

We may compare the structure of the above derivation to the intuitive argument presented at the beginning of this section. To the extent that α ≈ q (strictly true in the Rayleigh wave limit q → 0), Eq. 19 shows that we can indeed formally absorb the effects of surface tension into a q-dependent shear modulus, . Furthermore, referring to Eq. 25, we find that the threshold wave number for instability occurs at q ∼ 1/lγ, as expected.

Nonlinear theory

Derivation of the Mooney-Rivlin energy

In this section, we detail the derivation of Eq. 2, giving expressions for , , , and for the case of a uniaxial ellipsoid. We emphasize at the outset that the details of the results depend on the ellipsoidal geometry that we have chosen, but their structure does not. One may repeat these calculations for other starting geometries, for example, a cube, and obtain similar results.

Elastic deformations are described by the deformation gradient tensor Λ (). For a 3D material, is a 3D tensor that gives the local mapping of material points from the undeformed state x to the deformed state X. In general, Λ depends on x, the position within the undeformed state. Here, we assume a homogeneous deformation, for which Λ is constant. Given Λ, the three lowest-order rotational invariants that can appear in the elastic free energy density felastic are I1 = TrC, , I3 = DetC, where C = ΛΛ is the right Cauchy-Green deformation tensor (). I1 is the neo-Hookean term, and I3 describes volumetric deformations, i.e., I3 = 1 for incompressible solids, as we consider here. I2 is the Mooney-Rivlin term, often used to phenomenologically account for material nonlinearity in rubbers (). The elastic part of the free energy density can be written as

We consider a uniaxial deformation, , and let [μ is indeed the linear elastic shear modulus ()], . Equation 26 is then

For a uniaxial ellipsoid of radii the total elastic free energy iswhere we identify the first term in Eq. 28 as FNH and the second as FMR.

The surface energy iswhere is the eccentricity. For the bending energy Fbend, we use the Helfrich form discussed in the “Linear instability and shape change in an activeelastocapillary droplet” section

We now evaluate Eq. 30 for the case of a uniaxial ellipsoid. Finite c0 does not qualitatively change the structure of our results, and we consider c0 = 0 for simplicity. The area element iswhere u and v are the azimuthal and polar angles on the ellipsoid. The mean curvature H is

Equation 30 then simplifies towhich may be evaluated exactly; the result is

Combining Eqs. 28, 29, and 34, all divided by μR3, gives the total free energy Eq. 2.

Landau theory coefficients for the Mooney-Rivlin model

Here, we argued that the behavior of an active elastocapillary sphere near the critical point λ = 1 can be understood only on the basis of symmetries using the Landau expansion Eq. 3. We now derive the coefficients r, w, and u of Eq. 3 for the case of the Mooney-Rivlin free energy Eq. 2. Expanding Eqs. 28, 29, and 34 in the strain ϵ (=λ − 1), we obtain

Equation 35 has a critical point at , , , where

Expanding as , , and , we obtain the structure of the free energy about this critical pointwhere we omit terms such as Δγϵ4.

In Fig. 4, we show phase diagrams obtained from minimizing the exact free energy Eq. 2. We can interpret their structure in light of Eqs. 36 and 37, with a focus on the interplay between bulk elasticity and surface effects. The Landau theory of Eq. 37 corresponds to three weakly discontinuous transitions in the plane, meeting at a critical point (Fig. 4A). This structure is unchanged by varying the bending modulus . However, increasing the bending modulus drives the critical point to higher values of active driving and lower material nonlinearity (Fig 4B), enlarging the region of phase space in which worms are favored over pancakes. In this sense, both material nonlinearity and bending rigidity conspire to produce worm-like structures, as opposed to the more intuitively obvious pancake.

Fig. 4.

Structure of the worm/pancake phase diagram.

(A) Three weakly discontinuous transitions, described by the Landau theory Eq. 37, meet at a critical point given by Eq. 36. (B) As the bending modulus increases, the critical point is driven to larger active driving and lower material nonlinearity , enlarging the worm-like region of parameter space.

Numerics

In this section, we describe the microscopic ball-spring model and numerical methods used to realize the predictions of the continuum theory shown in Fig. 3.

Microscopic model

We first construct a disordered tetrahedral meshing of the ball, as shown in Fig. 5. We label the vertices i, edges ij, triangles α, and tetrahedra t. A microscopic energy for deformations of this mesh is given by a spring energy Fspring, a surface bending energy Fbend, and an approximate volume constraint Fvol

Fig. 5.

Microscopic ball-spring model.

(A) Schematic of the microscopic model at the surface of the meshed ball, showing vertices i, j, k…, edges ij, jk…, surface triangular plaquettes α, and volume tetrahedra t. (B) Cut-through of the meshed ball used in simulation, showing prestressed surface springs (blue) and unstressed bulk springs (orange).

For Fspring, we place a spring along every edge ij of the mesh. These springs have length r, and rest length , from which we define the extension . The spring energy is theni.e., each spring acts as an incompressible Mooney-Rivlin solid with microscopic neo-Hookean constant k and material nonlinearity α. To implement dilational surface stresses, we prestress the springs on the surface of the ball, initializing them at an extension λ < 1. The bulk springs are initialized at their rest length, λ = 1. We vary the macroscopic nonlinearity in material response by varying α for the bulk springs, keeping the surface springs at α = 0. A bending energy Fbend is given by ()where the sum is over neighboring triangular plaquettes α and β on the surface of the ball, with nα as the normal to plaquette α. Last, we approximately enforce incompressibility with an additional energetic penalty on volume changes of the tetrahedra t of the mesh ()where V is the current volume of a tetrahedron and is its rest volume.

The microscopic energy Eq. 38 contains k, α κ, λ, and B as microscopic parameters. We now describe a mapping to the continuum shear modulus μ, bulk modulus B, material nonlinearity , bending rigidity κ, and surface tension γ. Given a typical mesh length scale a, dimensional analysis gives

For an analytical estimate of the relation between κ and κ, we may calibrate using the continuum limit of the discrete bending energy Eq. 40 for a sphere, (). Comparing this to the continuum energy 8πκ gives the relation

For an estimate of the mapping from λ to |γ|, one can show that the energy per unit area of a triangular spring mesh of side length λ, composed of neo-Hookean springs, is given by

Numerical methods

The data in Fig. 3 are generated by numerically minimizing Eq. 38 at fixed k, α, κ, and B, with λ progressively decreasing from λ = 1 (giving progressively stronger dilational surface stresses). The final state of each minimization is then used as an initialization condition for the next. The minimizer used is the SciPy implementation of Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm, with mesh vertex coordinates as input and gradient norm stopping threshold of 10−2. For the data shown in Fig. 3, a ball of radius R = 1 (which defines the arbitrary spatial unit) is meshed with typical edge spacing a = 0.2. Microscopic parameters k = 0.013, B = 50,000, and κ = 2.5 are fixed for all runs. The three curves in Fig. 3E correspond to microscopic nonlinearities α = 0 (green triangles), 0.3 (blue circles), and 0.4 (orange squares). From the numerical data, we first map λ to γ using Eq. 46. To obtain the values of λ shown in Fig. 3E, an ellipsoid is then least-squares fit to the boundary vertices of the numerically relaxed mesh. The fit returns three ellipsoid axes, two of which are of similar magnitude (δλ/λ < 0.1), the third of which defines λ.

Last, we use the location of the critical point to fit the continuum theory Eq. 2 to this data, with , as fitting parameters. The theoretical fit shown in Fig. 3 corresponds to , via Eq. 36. An independent assessment of may be made using Eqs. 45 and 44, from which we estimate κ ≈ 0.72, μ ≈ 5 (independent measurement of μ for the meshed sphere used in simulation places μ ≈ 4). These estimates give , consistent with the value that we obtain from our fit.