Snap-induced morphing: From a single bistable shell to the origin of shape bifurcation in interacting shells

The bistability of embedded elements provides a natural route through which to introduce reprogrammability to elastic meta-materials. One example of this is the soft morphable sheet, in which bistable elements that can be snapped up or down, are embedded within a soft sheet. The state of the sheet can then be programmed by snapping particular elements up or down, resulting in different global shapes. However, attempts to leverage this programmability have been limited by the tendency for the deformations induced by multiple elastic elements to cause large global shape bifurcations. We study the root cause of this bifurcation in the soft morphable sheet by developing a detailed understanding of the behaviour of a single bistable element attached to a flat 'skirt' region. We study the geometrical limitations on the bistability of this single element, and show that the structure of its deformation can be understood using a boundary layer analysis. Moreover, by studying the compressive strains that a single bistable element induces in the surrounding skirt we show that the shape bifurcation in the soft morphable sheet can be delayed by an appropriate design of the lattice on which bistable elements are placed.


Introduction
The ability to change shape is as important to an emerging class of engineering applications as it is to biological organisms: just as animals and plants morph in response to external stimuli, soft robots must be able to change shape to adapt to different environments and complete different tasks (Alapan et al., 2020;Shah et al., 2021;. In both cases, our understanding of different artificial mechanisms through which this shape change can be achieved has exploded in recent years with examples including pneumatic inflation (Pikul et al., 2017;Siéfert et al., 2019), multi-material 4D printing (Boley et al., 2019), magneto-responsive elastomers (Zhang et al., 2021) and 3D-printed composites (Kim et al., 2018), designed director fields within liquid crystal elastomers (Aharoni et al., 2018), designed cuts in planar sheets (Celli et al., 2018;Liu et al., 2020), as well as swelling of patterned hydrogels (Wang et al., 2017).
In each of the above artificial examples, techniques have been developed that allow a particular three-dimensional target shape to be achieved by the actuation of an initially flat sheet. However, these techniques in general are only capable of generating one, or perhaps two, designed shapes. Nature, however, is able to adapt form repeatedly and in a variety of ways. To achieve something similar in artificial systems requires a means of reprogramming the shape.
There are two key hurdles to achieve this reprogrammability in practice: firstly, a means of reversibly actuating different elements of an object between two different states (Faber et al., 2020;Alapan et al., 2020). Secondly, interactions between neighbouring elements may lead to unwanted global deformation mode that make it impossible to reach a well-controlled state (Moessner and Ramirez, 2006;Gilbert et al., 2016;Siéfert et al., 2021): the system is then generically stuck in a local minimum, with many similar states 'nearby' -what may be termed 'soft modes'.
In this paper, we study a soft morphable sheet, illustrating another system in which shape can be changed by controlling the state of individual elements; we show how understanding the behaviour of these individual elements may yield new understanding, allowing the appearance of soft modes to be postponed.

C B
A Figure 1: A silicone chocolate mould is an everyday instance of a 'Soft Morphing Sheet': by popping individual hemispherical elements, the mould morphs from its initial flat state (shown in (A)) to different stable curved states ((B) and (C)). In (B) the sheet curves around its long axis (highlighted in red in panel A) in response to elements being popped along this axis. In (C) the sheet curves about its short axis (highlighted in blue in panel A) in response to elements along the short axis being popped. (The diameter of each spherical cap is 3.5 cm.) To achieve a fundamental element with two different stable (activated and natural, or deactivated) states, elastic snap-through is a natural candidate: curved elastic structures such as arches and shells often have two stable states and can be switched between these two states by the application of a suitable external load (Holmes and Crosby, 2007;Taffetani et al., 2018;. It is therefore somewhat natural to consider a simple elastic sheet with multiple bistable elements embedded at various locations within it. Conveniently, a similar system can be found in the silicone moulds used for moulding chocolate (see Fig. 1). Popping individual hemispherical dimples (the moulds) from their natural state to inverted state induces deformations in the neighbouring sheet that is reminiscent of the deformation induced by a localized dilation (Oshri et al., 2019;Plummer and Nelson, 2020;Oshri et al., 2020;Hanakata et al., 2022). However, as more snapping elements are activated, we see that the shape of the whole sheet changes globally from the initially flat state (Fig. 1A) to different stable curved states ( Fig. 1B and C), depending on which dimples one chooses to pop. In this sense, the silicone mould represents a "Soft Morphing Sheet". A rigid version of this idea was first presented by Seffen (Seffen, 2006(Seffen, , 2007, but the use of an elastomer here allows deeper shells to pop reversibly (without plastic deformation), thereby creating larger local deformations and, as a result, greater global shape changes. We seek to understand this deformation, and so shall not focus on the means of switching the individual elements between different states, merely noting that this could be achieved using pneumatic pressure (Gorissen et al., 2020) or externally applied magnetic fields (Chen et al., 2021), for example.
A similar system has been studied recently (Udani and Arrieta, 2022). This work showed that the system becomes frustrated with many co-existing states with similar energies: geometrical frustration of neighbouring elements leads to degeneracy of the global shape, which are sometimes called 'soft modes' of deformation (Moessner and Ramirez, 2006). The system can be switched between neighbouring frustrated states by gently teasing the system (e.g. by hand); however, this frustration cannot be eliminated. Examples of the different global shapes that can be achieved as a result of this frustration for dimples on square and triangular lattices are shown in Fig. 2. In the application to shape change, however, the aim is to control the state of the individual elements to be buckled (or not), and hence to understand how the macroscopic behaviour (i.e. global shape) emerges from the microscopic behaviour (i.e. local deformation). In this paper, therefore, we seek to understand the buckling instability that leads to shape bifurcation, as well as illustrating one strategy to delay its onset.
The paper is organized as follows. We start by analysing the individual snapping shell as a single element in Section 2. Here, we make use of shallow shell theory (presented in Section 2.2) to study first the conditions for bistability (Section 2.3). To understand the inverted state further we then present a boundary layer analysis of the deformed shape of a single element (Section 2.4); in each case, we compare our theoretical results with those from both Finite Element Method (FEM) simulations and experiments. The boundary layer analysis of Section 2.4 lends insight into the root cause for the shape bifurcation that is observed and results in frustration. We therefore study, in Section 3, the effect of lattice design on shape bifurcation and hence on frustration, showing that shape bifurcation can be delayed by placing elements on a hexagonal lattice, rather than a triangular lattice. Finally, in Section 4, the paper is concluded with some remarks and suggestions for future work.

The fundamental element: A single snapping shell
As a first step to understand the properties of the collective dimpled sheet, we consider the fundamental element i.e. a single spherical cap embedded within a planar thin plate -a 'skirt'. (Alternatively, the sheet consists of 'inclusions', the shells, embedded in a 'matrix'. For a single element, we prefer to retain the term skirt as a reminder that it has a circular boundary, and to differentiate it from the matrix between multiple elements in a sheet.) We call this fundamental element a 'snappit' (analogous to the bits used in electronic devices) and assume that its material and geometrical properties (save for the natural curvature) are identical to those of the remainder of the collective sheet; moreover, we shall assume that all the snappits of a dimpled sheet have the same characteristics. It should be noted that in this section we shall focus on the axisymmetric behaviours of a single element in §2.1- §2.4; we will turn to consider the causes of asymmetry in §2.5. The simplest geometry for a single snappit is to have an axisymmetric flat plate beyond the spherical cap, as shown in Fig. 3A. The geometrical parameters of this system are then given by Figure 2: The frustration of inverted dimples placed on a square lattice (top) and triangular lattice (bottom) leads to different global modes of deformation that can be easily shifted by hand ('soft modes'). (A) The natural (undeformed) state of the sheet with natural dimples on a square lattice (with total size of 159 × 159 mm) leads to at least two deformation shapes when the dimples are inverted, (A-i) and (A-ii). (B) The natural (undeformed) state of the sheet with natural dimples on a triangular lattice (169 × 157 mm) leads to at least two deformation shapes when the dimples are inverted, (B-i) and (B-ii). In each case these sheets are 3D printed as described in Appendix A.
those of the spherical cap (radius of curvature R, half-width L, thickness t) and of the flat skirt region (outer radius L + ∆L and, again, thickness t).
A typical example of this is shown in Fig. 3B and D: a 3D printed snappit can be popped from its natural state (see Fig. 3B and C) (Taffetani et al., 2018) to an inverted state (see Fig. 3D and E). This inversion induces a deformation of the plate that tends to be localized in the region near the interface with the shell.
Two natural questions emerge from this first picture of an inverted snappit: firstly, when is the inverted state stable (i.e. when is the snappit bistable?) and secondly, what is the induced deformation in the sheet?
We consider each of these questions in turn and will use a combination of experiments, simulations and theory to tackle them: physical experiments using 3D printed samples, FEM simulations using commercial software (ABAQUS) and detailed theoretical analysis using shallow shell theory. With this we obtain a quantitative comparison between simulations and experiments, as well as analytical and asymptotic results.

Stability of a shell: A review
To get an indication of the sort of behaviour that might be expected for a single snappit, we briefly review the scenario of a single spherical cap (without any extra material around the edge), which has been considered previously (Taffetani et al., 2018). In this case the stability of the inverted state of a spherical cap is independent of the Young's modulus E of the material: since there is no other force scale in the problem (assuming that gravity can be neglected) the transition between bistability and monostability is determined purely by the Poisson's ratio, ν, When pushed in, the spherical cap is able to invert, deforming the outer plate in the process (yellow shape). (B)- (E) Photos and Scanned 3D profiles of the snappit in the natural configuration (B and C) and the inverted configuration (D and E). Here L = 15 mm, R = 18.75 mm, t = 1 mm and ∆L = 2L (corresponding to α = 0.93, λ d = 6.95). Here U3,max = 9.0 and 8.0 mm for C and E, respectively. Details of the experimental protocols are given in Appendix A. and the geometrical properties of the shell, namely its radius of curvature R, half-width L, and thickness t.
On dimensional grounds, it is clear that with three length scales in the problem, two dimensionless groups can be formed. One natural choice is the angle sub-tended between the pole and the edge of the shell, i.e. α (= sin −1 (L/R)). Shells with α 1 are shallow (Ventsel and Krauthammer, 2001), but with α = O(1) are deep. To determine a second dimensionless parameter, Taffetani et al. (2018) used energy arguments to consider the competition between bending and stretching energies: the bending energy density induced by inversion is E B ∼ B(1/R 2 ) (where B = Et 3 /12(1 − ν 2 ) is the bending stiffness of the shell) while the stretching energy density may be estimated to be E S ∼ Etα 4 . The relative importance of stretching to bending energy densities is, therefore, E S /E B = 12(1 − ν 2 )(R 2 /t 2 )α 4 . We codify the relative importance of bending and stretching energies via the fourth root of this parameter (Libai and Simmonds, 2005), namely Clearly the parameter λ d involves both the depth of the shell, measured via α, and its slenderness, t/R. However, Taffetani et al. (2018) found that the transition between monostability and bistability occurs at a critical value, λ d = λ crit that is approximately independent of α for α 1.0 (from their results, the variation in λ crit is less than 1% for α 1.3). Finally, Taffetani et al. (2018) showed that the value of λ crit can be accurately calculated using shallow shell theory, which is only formally valid in the limit α 1. For the snappit with an attached outer skirt (see Fig. 3), there is an additional parameter, namely the size of the skirt ∆L measured relative to the size of the shell itself, L; we therefore introduce the ratio ∆ = ∆L/L.
We now turn to characterize how the presence of a skirt affects the bistability of the snappit, as well as the overall shape in the inverted state. Ultimately, we will be interested in this variation as all three of the dimensionless parameters (λ d , α and ∆) vary, but we initially simplify the problem by considering the problem for shallow shells, so that α 1.

Shallow shell formulation
A first description of the shape of an inverted element was presented by Sobota (2020); this involved employing a Rayleigh-Ritz approach with up to four degrees of freedom to solve the geometrically nonlinear shell model numerically. Here we provide some analytical insight by employing instead "shallow-shell" theory (Calladine, 1989;Ventsel and Krauthammer, 2001).

Governing equations.
Shallow-shell theory is based on the equations of axisymmetric plate theory modified to incorporate the finite radius of curvature of the shell. As a result, the same equations can be used to describe the vertical (normal) deflection, w(r), and stress potential, ψ(r), in a cylindrical polar geometry for both the shell and the skirt region if we introduce a spatially varying natural curvature where H(·) is the Heaviside step function. We therefore have and where Y = Et is the stretching modulus of the material and the stress potential ψ is the derivative of the Airy stress function (defined such that σ rr = ψ/r and σ θθ = dψ/dr).

Non-dimensionalization.
We follow Taffetani et al. (2018) by introducing dimensionless variables so that the governing equations (3)-(4) become and is the indicator function for the shell region (i.e. S = 1 in the shell region, 0 < ρ < 1, and S = 0 in the skirt region, 1 < ρ < 1 + ∆). This non-dimensionalization introduces two dimensionless parameters: ∆ as expected, and (Note that λ s , which characterizes the geometry of the shell, is the shallow shell version of the parameter, i.e. λ s ≈ λ d for α 1.) To be able to solve equations (5)-(6), we must specify appropriate boundary conditions at the boundaries, ρ = 0 (shell centre), 1 + ∆ (skirt edge) and, crucially, the join between the shell and skirt regions, ρ = 1.

Boundary conditions.
The appropriate boundary conditions at the shell centre, ρ = 0, and skirt outer edge, ρ = 1+∆, are relatively straightforward: at the shell centre, we use the symmetry and no displacement conditions W (0) = W (0) = U r (0) = 0, where U r is the horizontal (radial) displacement. (Note that since the governing equations only involve W , not W itself, there is some freedom in the choice of W (0).) At the outer boundary of the skirt, ρ = 1 + ∆, we assume there is zero bending moment, zero shear force and no radial stress (σ rr = 0).
The interface between the shell and skirt region at ρ = 1 deserves further discussion. We have immediately that the horizontal and vertical displacements are continuous here, as is the radial stress, so that It is also natural to assume that the bending moment and out-of-plane shear force are continuous. (The in-plane shear force, σ rθ = 0 from our assumption of axisymmetry.) The slope of the deflection at the boundary may, in principle, be discontinuous (depending on how this join is implemented); however, we assume that this join is rigid so that the angle between the tangent to the shell at its edge and the tangent to the skirt remains constant at α. Given the difference in natural shapes, this means that the slope of the displacement is continuous, i.e. [W ] 1 + 1 − = 0. The boundary conditions discussed above are summarized in terms of W , ψ and their derivatives in Table C.1 of Appendix C.
The Eqs. (5)-(6) can be solved subject to these boundary conditions using, for example, the multipoint boundary value problem solver bvp4c in MATLAB. We shall discuss the behaviour of the inverted solutions of these equations in Section 2.4 but first focus on determining when this inverted shape exists numerically, i.e. the critical condition for the bistability of the snappit. Note that this numerical system is controlled by just two dimensionless parameters, ∆ and λ s , since the shallow shell formulation assumes that the shells are shallow, i.e. α 1.

Bistability
Our first goal is to understand the geometrical conditions under which a single snappit is bistable, i.e. when does the inverted shape actually exist without the application of any external loads? As already discussed, there remain only three dimensionless parameters in the problem (α, λ and ∆), together with the material property ν. We present two approaches to determine the critical condition for bistability here: firstly, we use FEM simulations (see Appendix B for details) with a range of values of α, λ d and ∆ and determine whether the inverted state is stable (and hence the snappit is bistable). This gives a rough indication of the behaviour of the critical value λ crit d (∆, α; ν); our results are shown in Fig. 4A, with the dependence of the threshold λ crit d as a function of ∆ shown in Fig. 4B. Note that here we determine the threshold, λ crit d , as the mean of the smallest λ d with bistable behaviour and the largest λ d with monostable behaviour for given ∆; error bars are then used to record the difference between these two values. Secondly, we use arc-length continuation in AUTO (Doedel et al., 2007) to solve the shallow shell system (5)-(6) and determine the critical value λ crit s (∆) at which the inverted state ceases to exist. Our continuation calculation involves three key steps: first, the perfectly everted solution is used as an initial guess to the solution of the system (5)-(6) everted solution for an arbitrary, but large, value of λ s and ∆L/L = 3. The true solution is 'close' to the everted solution for large λ s , and so the true solution can readily be found using a relaxation method. Secondly, we perform a one parameter continuation analysis, reducing the value of λ s , but holding ∆L/L = 3 constant. This allows the branch of stable everted solutions to be followed until a fold bifurcation in the system is reached. The position of the fold bifurcation in terms of λ s corresponds to the λ crit at the monostable-bistable threshold for the given ∆L/L = 3. Finally, we perform a two-parameter continuation, following the fold while varying ∆L/L and λ s accordingly. This allows us to extract the value of λ crit for each ∆L/L down to ∆L/L = 0. This calculation gives the solid curve in Fig. 4B, which agrees very well with the FEM results. On the whole, Fig. 4 shows that the presence of a skirt decreases the threshold value of λ d at which the shell becomes bistable: compared with the case of a spherical cap without the skirt (λ crit (∆ = 0) ≈ 5.75), the presence of a skirt region makes the element 'more' bistable because, for a given ∆, all the geometrical configurations with λ d in between λ crit (∆) and λ d ≈ 5.75 are now within the bistable region. An interesting feature of the behaviour seen in the phase diagram ( Fig. 4A) is that, while there is some dependence of the critical shell depth λ d for bistability on the size of the plate region, ∆, this dependence is (perhaps surprisingly) small: over a wide range of ∆, the critical value of λ d changes by less than 10%. We also see that, as observed previously (Taffetani et al., 2018), the effect of the angle α = sin −1 (L/R) on the bistability is quite limited (indeed, it is not detectable in the results presented in Fig. 4B). As a result, we shall follow Taffetani et al. (2018): we use shallow shell theory to understand the importance of the parameter λ = λ s , and apply these results to less shallow shells (with α 1, not α 1) simply by letting λ = λ d .

Shape of an axisymmetric deformed element
Having determined when the element is bistable, we now move on to understand the shape in the inverted state. Given the success of shallow shell theory in describing the transition from mono-to bi-stable just demonstrated, it would be convenient to be able to use shallow shell theory to study in detail the deformed shape. First, however, we should confirm whether the effect of nonlinear elasticity is significant experimentally. In Fig. 5A we compare results from FEM simulations (see Appendix B for details), both using hyperelastic and linearly elastic constitutive models with results from physical experiments (see Appendix A). We see first that FEM results with each constitutive law are indistinguishable at the scale of the plot and, further, that either constitutive law does an equally satisfactory job of describing the experimental data. We shall therefore continue to study the shape assuming a linearly elastic model and, since we have also seen that the parameter α has a limited effect on the transition from mono-to bistability if we identify λ d → λ s , we use shallow shell theory. The numerical solution of the shallow shell equations gives a good account of the simulation results (see Fig. 5B); crucially, however, we will be able to gain some analytical insight in the limit of thin shells λ s 1, even though shallow shell theory is formally only valid for α 1.
Boundary layer analysis for λ s 1. Both our FEM simulations and numerical solutions of the shallow shell equations suggest that the deformation from the perfectly inverted shape becomes increasingly localized to the boundary between skirt and shell, at ρ = 1, as λ s increases (see Fig. 5B). We see from (5) that as λ s → ∞ the bi-Laplacian term becomes less important, except possibly in small regions where the derivatives become large -this term may therefore be neglected except in spatial boundary layers, which we expect to be centred on ρ = 1. Moreover, the 'outer' solution in the majority of the shell corresponds simply to "perfect eversion" within the shell region, W ∼ ρ 2 −1, and a constant vertical displacement within the skirt region, W ∼ cst; this deformation also results in zero stress within the majority of the shell. In this limit, therefore, the deformation is localized within the boundary layer, and the size of the boundary layer gives the scale over which the snappit deformation is expected to decay. We shall see that this length scale is ∼ λ −1 , which in physical variables corresponds to the Pogorelov length scale, p ∼ (tR) 1/2 , as seen in many other problems involving the inversion of spherical caps (Libai and Simmonds, 2005;Pogorelov, 1988;Gomez et al., 2016).
Firstly, we can decompose the displacement of the snappit as a perfect inverted shape plus a small perturbation, i.e. W (ρ) = S(ρ)(ρ 2 − 1) + ζ(ρ) For the moment, however, we proceed generally by letting Examining the transformed (5) in the plate region, and requiring terms to balance, we find that: while from (6) we have: Combining this with (10), we therefore find that µ + β = 2. To make any further progress requires information from the shell region; since this information only enters through the ρζ term in (6), we require this to enter at the same order as the other terms so that µ = β, and hence µ = β = 1, γ = 2.
Plotting the results of FEM simulations in terms of the boundary layer variables suggested by this analysis (specifically, plotting Z(ξ) = ζ(ρ) × λ s as a function of ξ = λ s (ρ − 1)), shows a good collapse as λ s → ∞ (see Fig. 5C). To understand the limiting behaviour, i.e. with λ s = ∞, we note that (9) transforms (5) and (6) to and These equations are to be solved subject to the condition that displacement and stress should vanish away from the join region; within the shell region we take Z s , Z s , Ψ s → 0, as ξ → −∞ (since the natural shell curvature has already been subtracted). Within the skirt region we take Z p , Z p , Ψ p → 0, as ξ → ∞, setting zero vertical displacement at infinity to fix the object in space. At the interface between shell and skirt, we also require continuity of these fields as well as the horizontal displacement i.e. 0 = Z s −Z p = Z s −Z p = Z s −Z p = Z s −Z p and 0 = Ψ s −Ψ p = Ψ s −Ψ p at ξ = 0. Solving (12)-(13) numerically subject to these condition yields the solid black curve shown in Fig. 5C. Note that this boundary layer solution is in good agreement not only with the results of FEM simulations (Fig. 5C), but also agrees well with experiments when cast back into the shape variables (Fig. 5D). This analysis shows that the join between the skirt and shell is the cause of the deformation observed in an inverted snappit: there is an incompatibility between the perfectly inverted shape (an inverted spherical cap joined to a flat skirt region) and the (fixed) angle between the shell and the flat region at the join. This geometrical origin of the boundary layer is distinct to the boundary layer region in an inverted spherical cap (Libai and Simmonds, 2005;Taffetani et al., 2018), which is caused by a mismatch between the zero applied moment condition and the finite, but constant, bending moment required to maintain the inverted spherical shape.
As well as showing that the join region is what is responsible for the deformation of the skirt region, the boundary layer analysis also gives us insight into the structure of the stress profile within the deformed snappit. In particular, the boundary layer scaling gives σ θθ (= dΨ/dρ) = O(λ −1 s ) while σ rr (= Ψ/ρ) = O(λ −2 s ); as a result, we expect that σ θθ σ rr , which is consistent with the numerical results. We also observe that both of these stresses are negative at some point within the vicinity of the join between the sheet and shell regions: the snappit is under both radial and azimuthal compression. In general, such a compression may be expected to lead to a buckling instability, and so we turn now to consider this instability.

Origin of instability
The analyses of Sections 2.3 and 2.4 were predicated on the assumption that the deformation remains axisymmetric. However, in some scenarios, a visible azimuthal instability occurs (see Fig. 6A,B). It is natural to treat this as a buckling instability and the boundary layer theory just presented gives two insights: firstly, both the hoop and radial stresses are negative (compressive) in some region close to the join between sheet and shell, see Fig. 6C,D. Secondly, in the axisymmetric state, the hoop stress is an order λ s larger than the radial stress (from Fig. 6C, σ rr ∼ λ −2 s , while from Fig. 6D, σ θθ ∼ λ −1 s ); therefore we might expect relieving the hoop compression by azimuthal buckling to be more energetically favourable than relieving the radial compression by buckling in the radial direction.
While this qualitative discussion explains the origin of the instability, to go further would require a more detailed analysis, along the lines of standard post-buckling analyses. On the face of it, the fundamental problem we consider has some similarities with the Lamé problem (Coman and Bassom, 2007;Davidovitch et al., 2011) in which an annular sheet is subject to different internal and external tensions; in particular, here an annular sheet (the skirt) is subject to a radial displacement at its inner edge and suffers an azimuthal instability as a result. The distribution of the scaled (C) hoop stresses (σ θθ = σ θθ · λ, blue symbols) and (D) radial stresses (σrr = σrr · λ 2 , red symbols) stresses within the snapped elements with λ d = 8.63 (squares) and λ d = 17.26 (circles) as calculated from FEM simulations, plotted as a function of ξ = (ρ − 1)λ. The corresponding predictions of the boundary layer theory (solid curves in both (C) and (D)) show that the compressive hoop stress |σ θθ | = O(λ −1 ), and hence is significantly larger than the typical compressive radial stress |σrr| = O(λ −2 ); this explains the azimuthal, rather than radial, buckling instability observed here.
However, there are some important differences with the standard Lamé problem. Most notably, in the Lamé problem, there are two dimensionless parameters: one corresponding to the ratio of the applied tensions (or imposed displacements) and the other related to the ease with which the sheet bends (the 'bendability' defined by Davidovitch et al. (2011)). As a result, for a fixed bendability one can consider gradually increasing the tension ratio until a small, buckled perturbation with a preferred wavenumber appears (as an eigenfunction of the problem). This is a 'Near threshold' analysis. In contrast, in the inversion-induced instability we report here, there is only a single dimensionless parameter (λ), which encodes both how easily bent the sheet is and the magnitude of the inward displacement that occurs at the inner edge. In some sense, the difference arises because there is no external tension applied at the outer boundary. This is similar to the case of an elastic ring subject to an internal tension, but no external tension, which has been shown to buckle with a wavelength λ ∼ (B/ T R 2 in ) 1/2 (Box et al., 2020), which would correspond to a mode number of instability n ∼ λ. However, the instability studied by Box et al. (2020) is driven by the dynamics of buckling and cannot occur statically -this is therefore not relevant to the static problem considered here. Instead, it seems that the mode number of instability is determined in some other way, as discussed in the tensionless Lamé problem by Pal et al. (2022). We do not consider this problem further here, leaving it for the future, but focus instead on how the system collectively chooses to relieve compression.

Many elements: Global response to distributed compression
The response of a single skirt region to the azimuthal compression caused by the presence of an inverted dimple is naturally to buckle azimuthally. When multiple elements are combined in a single sheet, i.e. the skirt region is shared between multiple elements to become a matrix, it is not clear what the global response to compression should be. We therefore begin by studying experimentally and numerically the response of the smallest symmetric system that contains multiple elements embedded within a circular sheet: a single snappit surrounded by six (hexagonally packed) other snappits. In this case, the snappits form a triangular lattice.

A triangular lattice
The fundamental geometry is shown in Fig. 7E. This geometry is characterized by the properties of each, identical, snappit (i.e. its thickness, radius of curvature and base width) but also the separation between the centres of each snappit, d. We focus here on the effect of the snappit separation, d, on the behaviour of the whole sheet. Heuristically, we might expect that if the separation between snappits is 'large' then they are effectively isolated with no significant interaction between them. Moreover, our analysis of the single snappit suggests that the boundary layer thickness BL = L/λ is the natural horizontal length scale over which the deformation induced by each snappit decays. Fixing λ, we therefore focus on the effect of the dimensionless gap width between snappit edges:δ Defined in this way,δ measures the half-spacing between two snappits in terms of the boundary layer width surrounding each snappit; we therefore expect that forδ sufficiently small the interaction between snappits should be strong, while for sufficiently largeδ the snappits are effectively isolated. Figs. 7A-D show comparisons between the shapes observed experimentally and in FEM simulations as the distance between snappits changes. The shapes predicted by FEM simulations agree well with those observed experimentally. As expected, both approaches show that as the separation between snappits increases the collective deformation of the sheet decreases. More interestingly, however, we see that it is not just the magnitude of deformation observed that changes with δ, but also the type of deformation: for smallδ the sheet deforms cylindrically (forming a taco-like shape, as in Fig. 7A) while for large δ the shape is (approximately) axisymmetric, see Fig. 7D for example. More specifically, in Fig. 7F we plot profiles, taken along the dashed magenta line, of snapped dimpled sheets with dimple spacing d varying in the range 32 mm ≤ d ≤ 56 mm . It is clearly seen that decreasing d is associated with a significant change in shape. We claim that the appearance of this buckling transition is crucial in the appearance of soft modes in the soft shapeable sheet: on a triangular lattice, the cylindrical mode has a multiplicity of three (the fold can happen along each of the three axes) and so multiple global shapes can be formed in a larger sheet.
To understand this global transition in shape, it is natural to draw analogies to other systems in which out-of-plane deformation occurs as a result of in-plane strain. While this is the fundamental origin of Euler buckling, two more specific examples are worthy of some discussion: deformation caused by inhomogeneous growth or shrinkage in a thin, single-layer material (Klein et al., 2007;Dervaux and Ben Amar, 2008;Sharon and Efrati, 2010) and that caused by differential thermal expansion in a bilayer material (Freund and Suresh, 2006). In both cases, deforming out of plane relaxes some of the in-plane strain when the system is flat, but the means of doing so, and conditions under which it occurs, are quite different. When the strain profile through the sheet thickness is uniform a radial inhomogeneity in growth/shrinkage must be introduced to obtain buckled shapes; moreover, bowl-like shapes (positive Gaussian curvature) are obtained from increasing shrinkage at the edge while azimuthally oscillating shapes (negative Gaussian curvature) are observed with increased shrinkage at the centre (Klein et al., 2007). In contrast, a bilayer with differential strains in the two layers adopts a bowl-like shape (positive Gaussian curvature) before transitioning to a more cylindrical shape (though still with positive Gaussian curvature) above a critical strain (Freund and Suresh, 2006). It is natural to assume that reducing the spacing between identical snapping elements (i.e. reducingδ) corresponds to increasing the level of strain in the material. Given this, the qualitative behaviour observed asδ varies in a dimpled sheet is essentially the same as that observed in the bilayer sheet: Fig. 7 shows that sheets with largestδ appear to be rotationally symmetric in the inverted state, while those with smallestδ clearly break the underlying rotational symmetry. This similarity to a bilayer sheet might be rationalized by the observation that the contraction applied by the individual dimples lies below the neutral plane of the unstrained parts of the matrix region between them.
To quantify the transition from rotational symmetry to two-fold symmetry, we measure the horizontal displacement perpendicular to the fold that ultimately forms (denoted u ⊥ r ), as well as the horizontal displacement measured at 60 • to this axis (denoted u ∠ r ) -see sketch in Fig. 7E for these definitions. (Quantifying the curvature in different directions is difficult because of the dimples' shape; we argue that measuring u ⊥ r and u ∠ r is both more reproducible and clearer.) Plotting these displacements as a function ofδ, see Fig. 7G, shows a sharp transition asδ decreases: forδ 10, we have u ⊥ r ≈ u ∠ r (see Fig. 7D), while forδ 6, u ⊥ r u ∠ r (Figs. 7A and B); there is also a notable intermediate range, i.e., 6 δ 10, in which u ⊥ r is slightly larger than u ∠ r ; this corresponds to the slightly curved state shown in Fig. 7C. We also note that this transition aroundδ ≈ 6 corresponds to the deformation and stress profiles within the two boundary layers of the interacting snappits no longer overlapping: from Fig. 5C and Figs. 6C and D we see that the effects of the boundary layer are localized within a region |r − L| 5 BL (|ξ| 5) of the join between element and skirt.

A hexagonal lattice
The cause of the bifurcation from rotational symmetry to the two-fold symmetry discussed above is the high level of strain within the sheet that is caused by the interaction of neighbouring elements. However, we find that, unlike the bilayer sheet, it is not simply the bare strain level that affects the point at which this bifurcation occurs. To demonstrate this, we consider removing a single snapping element from the centre of the sheet; this results in a hexagonal lattice as shown in Fig. 8E. The hexagonal lattice maintains the rotational symmetry of the triangular lattice already considered.
Figs. 8A-D demonstrate that the dimpled sheet based on a hexagonal lattice shows behaviour that is phenomenologically similar to that observed for a triangular lattice. In particular, the global shape retains rotational symmetry for large dimple spacingδ but at smaller values ofδ a bifurcation to a cylindrical shape occurs (compare Figs. 8A,B and C,D,for example). The profiles of the snapped dimpled sheets with d ∈ [30, 40] mm are also plotted in Fig. 8F. Comparing to the triangular lattice, the transition from the rotational symmetry to the two-fold symmetry is much sharper for the hexagonal lattice. We also measure the horizontal displacements along two different directions, u ⊥ r and u ∠ r , and plot them as functions ofδ in Fig. 8G. In this case, the bifurcation occurs withδ ≈ 3.2 (compared toδ ≈ 6 for dimples on a triangular lattice). Initially, this smaller separation at bifurcation might seem natural: there are only six dimples in the hexagonal lattice, rather than seven used in the triangular lattice -since the inversion of the dimples is the root of the strain and hence shape bifurcation, having fewer of them will thus delay the shape bifurcation. However, a more detailed analysis shows that even after accounting for this difference the hexagonal lattice delays the shape bifurcation, as we now show.

The effect of lattice spacing
The standard theory for the buckling of a bilayer sheet (Freund and Suresh, 2006) leads to the conclusion that the buckling transition occurs at a fixed strain level. It is therefore natural to compare the strain levels between the cases of triangular and hexagonal lattices. Since the strain level within the inverted dimpled sheet is far from uniform, we seek a measure that gives a value of the strain averaged over the whole sheet and use the total elastic energy, U e = 1 2 V σ ij ij dV = 1 2 E V 2 ij dV , of the sheet (with volume V ). We therefore introduce a measure of the volumeaveraged squared-strain, 2 ij defined as (15) Figure 9 shows how the anisotropy of the deformed shape, as measured by the difference between u ⊥ r and u ∠ r , varies as the volume-averaged strain, encoded by ψ e , varies. This comparison confirms that the shape bifurcation occurs significantly earlier (i.e. at lower strain) for the triangular lattice than for the hexagonal lattice. Crucially, this measure accounts properly for the fact that there are fewer snapping elements in this case so that the resulting effect can be attributed solely to the geometry of the lattice packing of snapping elements. Figure 9: (A) The normalized radial displacements (ūr) of the snapped dimpled sheets with both hexagonal and triangular patterns along perpendicular and inclined directions as a function of the strain energy density (ψe) obtained from FEM simulations. (B) The ratio between the normalized radial displacements along perpendicular and inclined directions (∆ū ⊥ r /ū ∠ r ) as a function of the strain energy density (ψe). Both quantities show a transition, corresponding to global shape bifurcation, at a critical strain energy density; this critical value is higher with a hexagonal lattice than for a triangular lattice.

Discussion and Conclusions
In this paper, we have considered the properties of bistable elements (snappits) embedded within a larger soft, shapeable sheet. Such embedded elements are able to induce a global deformation of the resulting dimpled sheet when they switch from their natural state to their other, inverted, state: in this way, 'activating' snappits allows them to induce deformation within the sheet.
We began by studying the deformation of a single, activated snappit in the framework of shallow shell theory. In particular, we showed that the lateral scale of the deformation induced by the inversion of a snappit scales with the Pogorelov length scale p = (tR) 1/2 . Moreover, by performing systematic FEM simulations, combined with the predictions of shallow shell theory, we found that the bistability of a single snappit is largely governed by the lateral size of the snappit relative to the Pogorelov scale, λ ∼ L/ p : the size of the outer skirt region has only a relatively weak influence on the threshold value of λ at the transition from monostable to bistable (see Fig. 4).
Our analysis of the deformation of a single activated snappit showed an azimuthal instability that breaks the rotational symmetry. Using a boundary layer analysis of the shallow shell equations, we showed that this azimuthal instability results from a large compressive hoop stress, though understanding the mode number of instability remains an open problem. We also showed that when several snappits are combined in a matrix, moderate compressive strains (achieved with moderately spaced snappits) lead to a global buckling mode that respects the rotational symmetry of the underlying lattice; this contrasts with the breaking of rotational symmetry in the single snappit case. Nevertheless, at larger values of the compressive strain (achieved with closely spaced snappits) the sheet deforms globally in a cylindrical shape, breaking the inherent rotational symmetry. Qualitatively, this transition occurs when the boundary layers surrounding each snappit no longer overlap and is similar to the shape bifurcation that is well known in a bilayer sheet with mismatch strain (Freund and Suresh, 2006). We also showed that the average level of strain at which the shape transition occurs can be increased by modifying the lattice on which snappits lie: on a hexagonal, rather than triangular, lattice, the strain at which buckling occurs is almost doubled. Controlling shape bifurcation in shape-shifting structures in this way could open the door to more robust morphing strategies.
CRediT authorship contribution statement

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.
geometrical since no external forces were imposed. Within FEM material constants of density ρ = 1.3 × 10 3 kg/m 3 , Young's modulus E = 6 MPa and Poisson's ratio ν = 0.4999 were taken for the majority of simulations, which were performed using the linearly elastic constitutive relation. To assess the importance of material nonlinearity effects, some simulations were also performed using the Neo-Hookean constitutive model with C 10 = 1 MPa and D 1 = 0.

Appendix B.1. Inversion of a single snappit
To model the axisymmetric deformations of a single snappit (Section II), we used 4-node hybrid bilinear axisymmetric quadrilateral elements (CAX4RH) for both the shell and the plate. Simulations were conducted in two steps (Static, General): first we applied pressure load to the top surface of the shell to invert the shell. The pressure was then removed in the second step, allowing the system to relax and find the closest equilibrium solution (an artificial global damping factor of η stab = 1 × 10 −8 was specified to stabilize the simulations). We consider two types of boundary conditions at the outer edge: (i) a hinged boundary, i.e. the vertical displacement and bending moment vanish (u z = M r = 0) at the outer edge (r = L + ∆L) and, (ii) a clamped boundary, i.e. both the vertical displacement and the radial rotation vanish (u z = ∂u z /∂r = 0) at the outer edge. For some combinations of geometrical parameters, the system returns to its initial, undeformed state after relaxation. In this case, the system is monostable (since no alternative equilibrium to the natural undeformed state has been found). If an alternative equilibrium is found in this way, the system is bistable. In the bistable scenario, the inverted shape is extracted from the bottom surface of the shape; this is consistent with what is measured experimentally.

Appendix B.2. Buckling of a single snappit and of a dimpled sheet
For the azuimuthal (non-axisymmetric) buckling of the snappit (in Section 2.5) and the global deformation of the dimpled sheet induced by the buckling of each dimples (in Section 3), we used 4-node doubly curved shell elements (S4R). The simulations were also conducted in three steps (here we chose the procedure type of Dynamic, Implicit): for a single snappit, we fixed the outer edge (r = L+∆L) and applied a displacement load at the centre (r = 0) in the first step, where the displacement is prescribed as twice the initial height of the shell (i.e. u z = W = 2(R−H)); we then removed the displacement load and released the fixed boundary condition at the outer edge in the second step, fixing the centre at the position as prescribed in the first step, allowing the system to relax; then a third step for relaxing the system is set to make sure the closest equilibrium solution is reached. For a dimpled sheet, the outer edge and central point were fixed, and a displacement load applied at the centre of each dimple (again, u z = 2(R − H)) in the first step; in the second step, we removed the displacement and released the fixed boundary condition at the outer edge to let the system relax and find the closest equilibrium solution. We then extracted the profile of the curved shape of the sheet in the snapped state, and measured the radius of curvature by fitting the extracted profiles.

Appendix C. Boundary conditions
To make the problem solved numerically in Section 2 clearer, we summarize the appropriate boundary conditions, expressed in terms of the unknowns of the problem, together with their physical meaning in Table C.1.