Dynamics of viscoelastic snap-through

We study the dynamics of snap-through when viscoelastic effects are present. To gain analytical insight we analyse a modified form of the Mises truss, a single-degree-of-freedom structure, which features an `inverted' shape that snaps to a `natural' shape. Motivated by the anomalously slow snap-through shown by spherical elastic caps, we consider a thought experiment in which the truss is first indented to an inverted state and allowed to relax while a specified displacement is maintained; the constraint of an imposed displacement is then removed. Focussing on the dynamics for the limit in which the timescale of viscous relaxation is much larger than the characteristic elastic timescale, we show that two types of snap-through are possible: the truss either immediately snaps back over the elastic timescale or it displays `pseudo-bistability', in which it undergoes a slow creeping motion before rapidly accelerating. In particular, we demonstrate that accurately determining when pseudo-bistability occurs requires the consideration of inertial effects immediately after the indentation force is removed. Our analysis also explains many basic features of pseudo-bistability that have been observed previously in experiments and numerical simulations; for example, we show that pseudo-bistability occurs in a narrow parameter range at the bifurcation between bistability and monostability, so that the dynamics is naturally susceptible to critical slowing down. We then study an analogous thought experiment performed on a continuous arch, showing that the qualitative features of the snap-through dynamics are well captured by the truss model. In addition, we analyse experimental and numerical data of viscoelastic snap-through times reported in the literature. Combining these approaches suggests that our conclusions may also extend to more complex viscoelastic structures used in morphing applications.


Elastic snap-through
Snap-through buckling is a striking instability in which an elastic object rapidly jumps from one state to another. Such instabilities are familiar from everyday life: umbrellas suddenly flip upwards on a windy day, while the leaves of the Venus flytrap store elastic energy slowly before abruptly snapping shut to catch prey unawares ( Forterre et al., 2005 ). Similarly, snap-through is harnessed to generate fast motions in technological applications ranging from fluidic actuators ( Gomez et al., 2017b;Overvelde et al., 2015;Rothemund et al., 2018 ), micro-scale switches ( Krylov et al., 2008;Ramachandran et al., 2016 ), Fig. 1. A jumping popper toy can be turned inside-out and released on a surface. It becomes unstable, and after a time delay ( ≈ 70 ms here) the popper rapidly snaps (in under 5 ms) back to its natural shape and leaps from the surface. responsive surfaces ( Holmes and Crosby, 2007 ) and artificial heart valves ( Gonçalves et al., 2003 ). In these applications, snapthrough has proved to be particularly useful among other elastic instabilities, such as wrinkling and crumpling, due to its ability to convert energy stored slowly into fast motions in a highly reproducible way.
Despite the ubiquity of snap-through in nature and engineering, its dynamics is not well understood, with classical work focussing on determining the onset of snap-through in simple elastic objects such as plates and shells ( Bazant and Cendolin, 1991;Patricìo et al., 1998 ). Because snap-through generically occurs when a system is initially in an equilibrium state that ceases to exist (a saddle-node/fold bifurcation), standard analytical techniques often cannot be used to study the dynamics. For example, it is generally not possible to perform a linear stability analysis to obtain an eigenvalue (natural frequency) that characterises the growth rate of the instability: beyond the fold point there ceases to be an equilibrium base state from which the system evolves. This is in contrast to the case when snap-through is caused by a bifurcation in which the equilibrium state becomes unstable without ceasing to exist ( Fargette et al., 2014;Pandey et al., 2014 ). The dynamics near a saddle-node bifurcation have been well studied in low dimensional systems, consisting of a few ordinary differential equations (ODEs), in various physical and biological settings ( Majumdar et al., 2013;Strogatz and Westervelt, 1989;Trickey and Virgin, 1998 ) including work on slow-fast systems -see Jones and Khibnik (2012) and chapter 2 of Berglund and Gentz (2006) (and references therein). However, it is much more difficult to extend this to an elastic continuum described by partial differential equations (PDEs). For this reason, previous work has mainly relied on experiments and numerical simulations (e.g. using commercially available finite element packages, or solutions of the governing PDEs using standard numerical methods) to quantitatively model the snap-through dynamics ( Arrieta et al., 2011;Brinkmeyer et al., 2012;Diaconu et al., 2009;Loukaides et al., 2014;Santer, 2010 ). Some progress has also been made using lumped mass-spring models ( Carrella et al., 2008 ), though there remains a general lack of analytical results in the literature; particularly lacking are closed-form expressions for the time taken to snap-through in terms of the physical system parameters. Analytical insight would be of interest both from the perspective of fundamental science and also for applications of snap-through, as such insight would provide a basis to control the dynamic response and guide more detailed simulations or experiments.
Moreover, some features of snap-through are not understood at a qualitative level, including delay phenomena: snapthrough often occurs much more slowly than would be expected for an elastic instability. This slowness is illustrated by children's 'jumping popper' toys, which resemble rubber spherical caps that can be turned inside-out. The inverted configuration remains stable while the cap is held at its edges, but leaving the popper on a surface causes it to snap back to its natural shape and leap upwards. As shown in Fig. 1 , the snap back is not immediate: a time delay is observed during which the popper moves very slowly, apparently close to equilibrium, before rapidly accelerating. The delay can be several tens of seconds in duration -much slower than the estimated elastic timescale, which is on the order of a millisecond ( Gomez, 2018 ).
To explain such discrepancies between estimates of the speed of snap-through and that actually observed, it is commonly assumed that some dissipation mechanism must be present. For example in the Venus flytrap, the estimated elastic timescale is orders of magnitude faster than the observed snap-through time, and air damping is not enough to account for the discrepancy. In this case the proposed mechanism is poroelasticity ( Forterre et al., 2005 ): the snapping leaves are saturated with water and may dissipate energy via internal fluid flow. Similarly, morphing devices that demonstrate delayed snap-through are commonly composed of silicone-based elastomers, which are known to exhibit viscoelastic behaviour ( Brinkmeyer et al., 2012 ). It is also easily demonstrated that holding a popper toy for longer in its inverted state causes a slower snap-back -an observation that is consistent with the importance of viscoelastic effects.
While attributing delayed snap-through to various dissipation mechanisms is natural, we recently demonstrated that anomalously slow dynamics are, in fact, possible in elastic systems with negligible dissipation ( Gomez et al., 2017a ). In such scenarios, the time delay arises from the remnant (or 'ghost') of the snap-through bifurcation, reminiscent of the 'critical slowing down' observed in other areas of physics such as phase transitions ( Chaikin and Lubensky, 1995 ) and electrostatic 'pull-in' instabilities : the saddle-node bifurcation continues to attract trajectories that are nearby in parameter space, producing a bottleneck whose duration increases without bound as the distance from the bifurcation decreases. Snap-through then appears to proceed much slower than what the elastic timescale alone would suggest. In the process, we were able to propose new analytical formulae for the snap-through time as a function of the material parameters. Nevertheless, a key feature of this slowing down is that the system needs to be very close to the snap-through transition: the amount of delay that is experimentally attainable may in practice be small. Moreover, in viscoelastic systems, it is not clear what role viscoelastic effects play in obtaining anomalously slow snap-through dynamics, as opposed to the purely elastic slowing down. While we have previously considered the influence of external viscous damping (e.g. due to air drag) ( Gomez, 2018 ), viscoelastic behaviour is fundamentally different because it modifies the stability characteristics of structures. Here, therefore, we seek to extend these studies to understand analytically how material viscosity affects the snap-through dynamics.

Viscoelasticity and pseudo-bistability
Unlike elastic solids, viscoelastic materials generally undergo stress relaxation when subject to a constant strain; this causes the effective stiffness of the structure to evolve in time. If a constant stress is imposed instead, the material may also exhibit a slow creeping motion ( Howell et al., 2009 ). Santer (2010) has demonstrated how these combined effects allow structures to exhibit 'temporary bistability' or 'pseudo-bistability' during snap-through. The idea of pseudo-bistability is that when a structure is held in a configuration that is near (but just beyond) a snap-through threshold, just as a popper toy may be held inside-out, the change in stiffness associated with stress relaxation may cause the structure to appear bistable (i.e. an elastic structure with the same instantaneous stiffness would be bistable). When the structure is released, the stiffness recovers during a creeping motion, until eventually this bistability is lost and rapid snap-through occurs. Similar to the phenomenon of creep buckling ( Hayman, 1978 ), the total snap-through time is then governed by the viscous timescale of the material and can be very large. This phenomenon may be useful in morphing devices that are required to cycle continuously between two distinct states; for example, dimples proposed for aircraft wings that buckle in response to the air flow to reduce skin friction ( Dearing et al., 2010;Terwagne et al., 2014 ), and ventricular assist devices which use snapthrough of a spherical cap under a cyclic pneumatic load to pump blood ( Gonçalves et al., 2003 ). In these applications, pseudo-bistability means that the actuation needed to move the structure between different states can be applied for a shorter duration, which may lead to a significant reduction in the energy consumed ( Santer, 2010 ).
Using finite element simulations, Santer (2010) has demonstrated pseudo-bistability in a single-degree-of-freedom trusslike structure, as well as spherical caps similar to the jumping popper toy of Fig. 1 . The phenomenon has been observed experimentally in spherical caps ( Madhukar et al., 2014 ) and truncated conical shells ( Urbach and Efrati, 2017 ), and generically appears to occur only in a narrow parameter range near the transition to bistability, i.e. the threshold at which snap-through no longer occurs. Brinkmeyer et al. (2012) performed a systematic study of the snap-through dynamics of viscoelastic spherical caps, using a combination of finite element simulations and experiments. Continuing this work, Brinkmeyer et al. (2013) studied the pseudo-bistable effect in viscoelastic arches. In these studies the phenomenon is found to have a number of common features, including (i) to obtain any time delay the structure needs to be held for a minimum period of time in an inverted state before release; and (ii) the resulting snap-through time depends sensitively on the parameters of the system and appears to diverge at the bistability transition. However, these basic features are not well understood quantitatively despite having important implications for applications of pseudo-bistability. The sensitivity of the snap-through time, for instance, means the system needs to be precisely tuned to obtain a desired response time. For this reason direct comparison between experiments and finite element simulations has revealed large quantitative errors ( Brinkmeyer et al., 2013;2012 ).
In addition, the numerical simulations referred to above are all based on two key assumptions regarding the viscoelastic response: (i) the material behaves elastically with an effective stiffness that evolves in time, and (ii) the response during recovery is the reverse of the response during relaxation, i.e. once the structure is released, the stiffness smoothly recovers to its initial, fully unrelaxed, value. These assumptions mean that modelling the dynamics is relatively simple compared to more general viscoelastic models, and the resulting equations are more easily implemented in commercially-available finite element packages. Furthermore, the different dynamical regimes can often be inferred by considering the elastic response in which the stiffness is fully unrelaxed and fully relaxed, as the instantaneous stiffness must be bounded between these two extremes ( Santer, 2010 ). However, the validity of these assumptions, and whether they can be justified from first principles, remains unclear.
An alternative approach is to start from the constitutive law of a viscoelastic solid, and derive the equations of motion that couple the stress to the deformation of the structure. While this approach is significantly more complicated, it eliminates the need to make any additional assumptions regarding the behaviour of the stiffness. This method has previously been used to obtain analytical expressions for the snap-through loads of simple viscoelastic structures ( Nachbar and Huang, 1967 ), as well as the conditions under which creep buckling occurs ( Hayman, 1978 ). More recently, Urbach and Efrati (2018) developed a general theoretical framework for modelling viscoelastic snap-through based on a metric description of the constitutive equations. While this approach yields insight into the phenomenon of pseudo-bistability, the dynamics are modelled quasi-statically by neglecting the system inertia, so that it is unclear precisely when pseudo-bistable behaviour is obtained. Elsewhere, due to the inherent complexity of viscoelastic effects, it is unknown what role inertia plays in the dynamics and why the snap-through time appears to diverge near the snap-through transition. Are we simply observing another instance of critical slowing down, similar to the purely elastic dynamics studied by Gomez et al. (2017a) ?

Summary and structure of this paper
In this paper, we aim to provide analytical understanding of the dynamics of viscoelastic snap-through, and in particular the features of pseudo-bistability. We consider a thought experiment in which a structure is indented to a specified displacement, and allowed to undergo stress relaxation before the indenter is abruptly removed. While we are motivated by continuous viscoelastic structures such as shells and arches, we first study a Mises truss for simplicity. This is a single-degreeof-freedom structure that exhibits bistability and snap-through, and enables us to make significant analytical progress. Focussing on the limit in which the timescale of viscous relaxation is much larger than the characteristic elastic timescale, we obtain three key results. (1) Inertial effects immediately after the indentation force is removed play an important role in determining when snap-through and pseudo-bistability occur.
(2) While the intuitive picture of pseudo-bistability as being caused by a temporary change in stiffness is correct, the assumption of reversibility made in previous numerical studies (i.e. that the stiffness smoothly reverses back to its fully unrelaxed value when the indenter is removed) leads to significantly different predictions of when snap-through occurs, compared to our first principles analysis. (3) Pseudo-bistability is a type of creeping motion governed by the viscous timescale, and so does not rely on critical slowing down to obtain slow dynamics, unlike purely elastic snap-through ( Gomez et al., 2017a ). Nevertheless, this creeping motion may be very slow indeed as the system may, in addition, be subject to critical slowing down in the pseudo-bistable regime. We then study a pre-buckled viscoelastic arch as an example of a more realistic structure that is used in morphing applications. Using direct numerical solutions, we show that the predictions of the truss model are qualitatively accurate for the arch system. This suggests that the analytical insight gained from the truss model may apply more broadly to the complex viscoelastic structures used in applications of snap-through.
The remainder of this paper is organised as follows. We begin in §2 by deriving the equations governing the motion of the Mises truss. We then discuss the equilibrium states and the stress relaxation during indentation. In §3 , we analyse the snap-through dynamics when the indenter is released, focussing on the limit in which the timescale of viscous relaxation is much longer than the characteristic elastic timescale. Using direct numerical solutions, we identify the different dynamical regimes and explain these asymptotically using the method of multiple scales. In §4 , we perform simulations of a viscoelastic arch, showing that its snap-through behaviour is well captured by the truss model. In §5 , we compare our predictions to experimental and numerical data of pseudo-bistable snap-through times reported in the literature. Finally, in §6 , we summarise our findings and conclude.

A simple model system: Mises truss
As a first step towards understanding the dynamics of viscoelastic snap-through, we follow Brinkmeyer et al. (2013) and consider a Mises truss (also referred to as a 'von Mises truss'). In its simplest (elastic) form, this features two central springs, assumed to be linearly elastic, that are pin-jointed at their ends and inclined at a non-zero angle to the horizontal in their natural state. To give the system inertia, we place a point mass where the springs meet, and restrict the mass to move only in the vertical direction; see Fig. 2 a.
The truss in its current form is bistable: as well as the undeformed or ' natural ' state, the truss may be in equilibrium in a reflected state, where the length of each spring is unchanged from its rest length. However, by connecting an additional spring of sufficient stiffness vertically to the point mass ( Fig. 2 b) ( Krylov et al., 2008;Panovko and Gubanova, 1987 ), the inverted state ceases to be a stable equilibrium: in an experiment in which the truss is held fixed in an inverted position using an indenter, the truss will immediately snap back to its natural state when the indenter is released, independently of how long it is held. This snap-through is reminiscent of a spherical cap ( Taffetani et al., 2018 ); in fact, we may consider the truss as a lumped model for a generic, continuous elastic structure that features an ' inverted ' state that snaps back to a ' natural ' state. The central springs represent the membrane (stretching) stiffness, since these springs can be viewed as corresponding to the midsurface of the structure. The vertical spring models the bending stiffness ( Krylov et al., 2008 ): this spring penalises rotating the truss about its pin-jointed ends, mimicking bending the structure about its edges as it is loaded to an inverted position.
We now suppose that the vertical spring is viscoelastic. To choose an appropriate viscoelastic model, we note that a typical snap-through experiment includes both displacement-control and force-control: during indentation we impose a given displacement, but releasing the structure corresponds to imposing zero indentation force. It is therefore insufficient to describe the viscoelastic response using a Kelvin-Voigt or Maxwell model, since these fail to accurately capture both stress relaxation (under displacement-control) and creep (under force-control) behaviour. Instead, we shall use the constitutive law of a standard linear solid (SLS), which is the simplest model that describes both of these effects ( Lakes, 1998 ). Physically, the SLS model is equivalent to placing a linear spring in parallel with a Maxwell element that features a second spring and While we could also incorporate viscoelasticity of the central springs in our formulation, this would introduce additional viscous timescales and hence make it much more difficult to make analytical progress; we will show that our simpler model successfully captures snap-through and pseudo-bistability without additional complexity.
When the truss is indented and held for a specified time, stress relaxation causes the effective stiffness of the SLS element to decrease, so that the behaviour upon release is no longer obvious: the truss may immediately snap back, or it may initially creep in an inverted state for a period of time. In particular, these regimes cannot be inferred by only considering the equilibrium states of the system. We note that Nachbar and Huang (1967) have analysed a similar truss using a Kelvin-Voigt model, and determined the onset of snap-through to an inverted state when a constant indentation force is suddenly applied. Here, we are interested in the dynamics of the snap-back when the indentation force is removed.

Governing equations
As shown in Fig. 2 c, in the natural state the central springs are assumed to be inclined at an angle α 0 > 0 to the horizontal, and the springs are at their rest length; the distance between the pin joints at each base is 2 w 0 . We assume that the central springs are linearly elastic with constant stiffness k . The rest length of the vertical SLS element is h 0 , the dashpot has viscosity η, and the upper springs have modulus E 1 (for the spring in parallel with the dashpot) and E 2 (for the spring in series).
Let x be the downward displacement of the point mass m from the natural state. To obtain an equation of motion for x , we calculate the various forces exerted on the point mass. We write α and l for the corresponding inclination angle and change in length of the central springs, respectively. For simplicity, we assume that the truss remains shallow in shape, i.e. α 0 1 and | α| 1. Neglecting terms of O (α 3 0 , α 3 ) , simple geometry gives that ( Gomez, 2018 ) l ≈ Because the central springs are linearly elastic, the vertical component of the total force exerted on the point mass by the central springs (directed downwards) is approximately 2 k α l .
The displacement x also leads to a strain in the upper SLS element of size e = x/h 0 . The corresponding stress σ satisfies the constitutive law of a standard linear solid ( Lakes, 1998 ) (with t denoting time) (1) Note that the limiting case of a Maxwell material is recovered by setting E 1 = 0 , while the constitutive law for a Kelvin-Voigt material is recovered in the limit E 2 → ∞ . The stress σ in the SLS element leads to a vertical force −Aσ exerted on the point mass (directed downwards), where A is the cross-sectional area of each element.
Combining the above, and also accounting for a downwards indentation force f , conservation of momentum gives Together with appropriate initial conditions specified below, the coupled ODEs (1) -(2) (together with the relation e = x/h 0 ) provide a closed system to determine the trajectory x ( t ) and stress σ ( t ).

Non-dimensionalisation
To make the problem dimensionless, it is natural to scale the displacement with the initial height of the truss in the small-angle approximation, i.e. x ∼ α 0 w 0 . We scale time with the characteristic timescale of stress relaxation, t ∼ η/ E 2 , obtained by balancing the final two terms in (1) . Balancing the remaining terms in Eqs.
(1) -(2) , we introduce the dimensionless variables Here we have chosen the stress scale so that the constitutive equation for an elastic solid is simply = X in dimensionless variables. Inserting these scalings into (1) , and eliminating the strain for the dimensionless displacement X , we obtain where we define The parameter β ∈ (0, 1) plays a key role in the stability of viscoelastic structures ( Urbach and Efrati, 2018 ), as it measures the degree of stress relaxation that occurs in response to a step increase in strain; we shall discuss this further below when considering the behaviour of the truss when the indenter is applied. The parameter (1 − β ) may also be interpreted as the ratio of the timescale of stress relaxation ( t ∼ η/ E 2 ) to the characteristic timescale over which creep , obtained by balancing the first two terms in Eq.
(1) . The related parameter E 2 /E 1 = β/ (1 − β ) is also referred to as the relaxation strength. The value of β is governed by the physical mechanisms causing viscoelastic behaviour, such as molecular processes (e.g. molecular rearrangement in polymers) or the effects of coupled field variables (e.g. fluid flow in poroelastic materials); for a detailed discussion see Lakes (1998) . Here we assume that β is a known material constant, which can be measured experimentally using relaxation tests ( Urbach and Efrati, 2017 ). We also note that in this nondimensionalisation, the limit β → 1 corresponds to a Maxwell material, since this is equivalent to setting E 1 = 0 in Eq. (1) .
This limiting case more closely resembles fluid-like behaviour in which the material has no preferred natural state and simply relaxes to the current configuration ( Urbach and Efrati, 2018 ). The opposite limit, β → 0, corresponds to a purely elastic material, in which stress relaxation does not occur and the solution of (3) is simply = X for all times. Note that the Kelvin-Voigt model (i.e. sending E 2 → ∞ in (1) ) cannot be obtained in this non-dimensionalisation, since we have scaled time by the relaxation timescale η/ E 2 . (This limiting case can only be obtained by first rescaling time by the creep timescale before sending β → 1.) In terms of dimensionless variables, the momentum Eq. (2) can be written as where we have introduced the dimensionless parameters Here, the Deborah number De measures the ratio of the timescale of stress relaxation to the characteristic timescale of the experiment ( Howell et al., 2009 ), which here is the timescale of elastic oscillations ( ∼ α −1 0 m/k ). (Hence, in this nondimensionalisation, the viscous timescale is T = O (1) while the elastic timescale is T = O ( De −1 ) .) We may interpret λ as the relative stiffness of the upper SLS element compared to the central springs. The cubic term on the right-hand side of (5) represents the dimensionless force due to the central springs. As expected, this vanishes in the undeformed state X = 0 , the reflected state X = 2 (when the central springs are also at their natural length), and the intermediate displacement X = 1 when the springs are aligned horizontally -in this state they are compressed but do not contribute any vertical force.

Steady solutions
When the system is in equilibrium with elastic constitutive law = X, the momentum Eq. (5) implies that the indentation force F must balance the total force exerted by the central springs and the SLS element, which we label F eq . In particular, the indentation force associated with a steady displacement X is When the indentation force is removed, any equilibria must satisfy F eq = 0 , which has roots For λ < 1/4, there are two real non-zero solutions, which coincide and disappear at a saddle-node bifurcation when λ = λ fold = 1 / 4 ; the corresponding displacement at this point is X = X fold = 3 / 2 . For λ > λ fold , the only real solution is the undeformed state, X = 0 . This behaviour is apparent in Fig. 3 a, which plots the force-displacement curve for different values . At zero force, the truss is bistable for λ < λ fold = 1 / 4 and monostable for λ > λ fold . (b) Response diagram for the steady roots of F eq (X; λ) = 0 as λ varies. At the critical value λ = λ fold = 1 / 4 , the stable non-zero root (upper solid curve) meets an unstable root (dotted curve) and disappears at a saddle-node bifurcation. of λ; we see that increasing λ (corresponding to a stiffer SLS element) acts to rotate the curve anti-clockwise about the origin, until eventually the turning point of the cubic lies above the line F eq = 0 . The corresponding behaviour of the roots to F eq = 0 is shown in Fig. 3 b. It can be shown that the roots for which F eq (X; λ) > 0 where = d / d X (solid branches on Fig. 3 b) are linearly stable, while the root for which F eq (X; λ) < 0 (dotted branch) is linearly unstable ( Panovko and Gubanova, 1987 ).

Indentation response
In a snap-through thought experiment, we imagine indenting the truss to an inverted state by imposing the constant displacement X = X ind ≥ 1 , for a time interval −T ind < T < 0 of duration T ind > 0 (for later convenience, we define T = 0 to be the time at which the indenter is released). To avoid introducing additional timescales into the problem, we suppose that the indentation is suddenly applied at T = −T ind , i.e. over a timescale much faster than the viscous timescale η/ E 2 . We can then approximate the behaviour for T < 0 as where H ( ·) is the Heaviside step function. Substituting this into the constitutive Eq. (3) , and solving for the stress in the upper SLS element, we obtain This solution is classical in the literature and represents the stress relaxation of a standard linear solid under a step increase in strain ( Lakes, 1998;Santer, 2010 ): the stress initially (i.e. at T = −T ind ) jumps instantaneously to a fully unrelaxed value = X ind / (1 − β ) when the indentation is applied, and then decays exponentially to the fully relaxed value = X ind associated with an elastic material.
Inspecting the momentum Eq. (5) , we see that the effect of this relaxation is to give an effective value of λ that changes in time. The corresponding indentation force can be written as where the effective value of λ is Note that λ eff decreases from λ eff (−T ind ) = λ/ (1 − β ) to λ eff (∞ ) = λ during relaxation. In terms of the force-displacement curve in Fig. 3 a, this corresponds to rotating the curve clockwise as stress relaxation occurs, so that the indentation force decreases in time. From this picture, we anticipate that there are different dynamical regimes when the indenter is released, depending on the values of λ and T ind . For λ < (1 − β ) λ fold , the turning point on the cubic lies below the line F eq = 0 in the fully unrelaxed state (since λ eff (−T ind ) < λ fold ), and moves further below this line as relaxation occurs. Hence the truss is bistable at the moment when the indenter is released, and we do not expect snap-through to occur if the indentation displacement is sufficiently close to the stable non-zero root of F eq = 0 . Similarly, for λ = λ eff (∞ ) > λ fold , the turning point lies above the line F eq = 0 when the structure is fully relaxed, and so the truss is always monostable; we expect snapthrough to occur for any value of X ind and T ind . For (1 − β ) λ fold < λ < λ fold , the turning point lies above the line F eq = 0 when the structure is fully unrelaxed (since λ eff (−T ind ) > λ fold ), but eventually decreases below this line as stress relaxation occurs. In particular, the truss is effectively bistable when the indenter is released (i.e. T = 0 ) provided that λ eff (0) < λ fold , which can be re-arranged to give We then expect snap-through to generally not occur if the inequality (6) is satisfied, and to occur otherwise. We will show that while this naïve argument correctly accounts for different dynamical regimes, it fails to quantitatively predict when snap-through occurs because of the effects of inertia. For later reference, we shall write F ind for the value of the indentation force F just before the indenter is released, i.e. at T = 0 − . From above, this is given by

Dynamics of release
At T = 0 , the indenter is suddenly released so that the indentation force We solve the momentum Eq. (5) for the corresponding stress and substitute this into the constitutive Eq. (3) . After rearranging, we obtain Due to the presence of inertia, X and ˙ X must be continuous across T = 0 (writing ˙ = d / d T ), giving the initial conditions The jump in acceleration here is necessary to balance the discontinuity in the applied indentation force.
Currently, we have five dimensionless parameters in the problem: the Deborah number De, relative stiffness λ, relaxation parameter β, indentation depth X ind , and indentation time T ind . Throughout this paper, we restrict our attention to indentation depths 1 ≤ X ind ≤ 2. As baseline values we use β = 1 / 2 (i.e. both springs in the SLS element have equal modulus, E 1 = E 2 ) and X ind = X fold = 3 / 2 , i.e. the displacement at the saddle-node bifurcation; in this case the initial conditions (9) are analogous to those studied by Gomez et al. (2017a) for purely elastic snap-through. We expect to recover similar behaviour here in the elastic limit β → 0, i.e. we expect the dynamics are governed by the elastic timescale and only slow down considerably near the saddle-node bifurcation at λ = λ fold . However, for values β > 0, it is not clear when the dynamics are instead governed by viscous relaxation. To gain insight, we focus on the limit De 1, which corresponds to a relaxation timescale that is much slower than the elastic timescale. This is the relevant regime for many structures composed of rubbery polymers, such as silicone-based elastomers typically used in morphing devices ( Brinkmeyer et al., 2013;2012;Urbach and Efrati, 2017 ), where the molecular rearrangement underlying viscoelastic behaviour occurs over slow timescales ( Lakes, 1998 ).

Numerical solution
Typical dimensionless trajectories X ( T ) in the limit De 1 are shown in Figs. 4 a-d. These are obtained by integrating the ODE (8) with initial conditions (9) numerically in matlab (routine ode45 , error tolerances 10 −10 here and throughout). Fig. 4 shows that the initial jump in acceleration causes oscillations to occur on the fast elastic timescale T = O ( De −1 ) ; these oscillations persist due to the absence of external damping in our model. As anticipated from the discussion in §2.4 , there are different regimes depending on the size of λ. For λ (1 − β ) λ fold ( = 1 / 8 with β = 1 / 2 ), the truss appears never to snap and instead approaches the stable non-zero root of F eq = 0 ( Fig. 4 a). For (1 − β ) λ fold λ λ fold , the truss snaps back to the natural state for small enough values of T ind , but remains in an inverted state indefinitely for larger T ind ( Fig. 4 b).
For λ ࣡ λ fold , the truss appears to snap for any value of T ind . However, the dynamics slow down considerably when 0 < λ − λ fold 1 and T ind is sufficiently large; see Fig. 4 c. In this regime the oscillations are rapidly damped out, and the trajectory features an initial plateau before abruptly accelerating towards the natural configuration (highlighted in the lower panel of Fig. 4 c), reminiscent of the dynamical bottleneck caused by a saddle-node ghost ( Gomez et al., 2017a ). For larger values of λ, the dependence on T ind decreases and this initial bottleneck phase is not observed ( Fig. 4 d).
These regimes are confirmed when we analyse the snap-through time, T snap (defined as the time at which the displacement first crosses the natural displacement, X = 0 ); the computed times are shown on the ( λ, T ind )-plane for the baseline parameter values in Fig. 5 a. The blank regions on the figure correspond to regions where snap-through does not occur (after integrating the equations up to T = 50 , which was found to be sufficient due to the limited amount of slowing down in Fig. 5 a). This shows that the critical value of T ind at which snap-through no longer occurs with (1 − β ) λ fold λ λ fold Fig. 4. Dimensionless trajectories obtained by numerical integration of (8) with initial conditions (9) (coloured curves). Here X ind = X fold , De = 10 , β = 1 / 2 and data is shown for (a) increases nonlinearly as λ increases, and appears to approach a finite value T ind ≈ 4 as λ → λ fold . For comparison, we have also plotted the naïve prediction (6) based on whether the truss is effectively bistable at the moment when the indenter is released (green dashed curve). This provides a good approximation for smaller values of λ, but increasingly over-predicts the critical value of T ind as λ increases, with the predicted value diverging as λ → λ fold . (For later comparison, the boundary predicted by the multiple-scale analysis in §3.2 is shown as a red dotted curve).
Another key feature of Fig. 5 a is that the snap-through time is very small throughout most of the parameter space. In fact, we will show that here the elastic oscillations cause the truss to immediately cross  Fig. 5 b, which shows that considerable slowing down can occur. In fact, the snap-through time appears to increase without bound as we take λ λ fold in this region. We will show that this is precisely the pseudo-bistable regime: here the displacement initially oscillates around an inverted state and does not immediately cross X = 0 . As with the trajectories in Fig. 4 c (lower panel), this inverted state also undergoes a slow creeping motion until the truss rapidly accelerates towards the natural state, so that T snap ࣡ O (1). This difference in timescales (i.e. a slow creep followed by a rapid snap-through event) is considered to be a distinguishing feature of pseudo-bistable behaviour ( Brinkmeyer et al., 2013;2012 ). Computed snap-through times for different values of X ind are shown in Figs. 6 a-b. These show that the boundary at which snap-through no longer occurs is qualitatively different depending on whether X ind < X fold or X ind ≥ X fold . For a shallower indentation X ind < X fold , the boundary appears to be shifted entirely to the left of the line λ = λ fold , and there is no longer a region where the dynamics slow down considerably ( Fig. 6 a). The truss also snaps at values λ < (1 − β ) λ fold when T ind is sufficiently small. In contrast, for deeper indentations X ind ≥ X fold , the boundary intercepts the line λ = λ fold and the The critical values λ = (1 − β ) λ fold and λ = λ fold are plotted as vertical black dotted lines. For later reference, also shown is the boundary predicted by Eq. (21) relevant for X ind < X fold (purple dashed curve), and the boundary predicted by Eq. (22) relevant for X ind ≥ X fold (red dotted curve). In each panel, the snap-through times have been computed on a 100 × 100 grid of equally spaced values. For ease of comparison the range of the colourbar is the same in both panels here and in Figs. 5 a-b. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.) size of the pseudo-bistable region may be significantly larger compared to the case X ind = X fold ( Fig. 6 b). For different values of the relaxation parameter β, we observe a qualitatively similar picture provided β ࣠ 1/2; see Appendix A .

Multiple-scale analysis
To understand the above observations, we now perform a detailed analysis of the dynamics in the limit De 1. The trajectories in Figs. 4 a-d indicate that the displacement undergoes fast oscillations (on an O ( De −1 ) timescale) around a value that varies on an O (1) timescale. This suggests that the dynamics can be understood asymptotically using the method of multiple scales ( Hinch, 1991 ). We introduce the fast timescale T defined by T = De −1 T . Treating T and T as independent, the chain rule implies that We seek an asymptotic expansion of the solution in the form

Leading order problem
We insert the expansion (11) into the ODE (8) , re-scaling time in terms of T . After Taylor expanding the F eq and F eq force terms about X 0 , and expanding the derivatives using (10) , we obtain to leading order in De −1 the homogeneous problem Integrating the above equation for X 0 with these conditions, and simplifying using the expression (7) for F ind , we obtain where A ( T ) is unknown and satisfies A (0+) = 0 .
To reveal the role that A ( T ) plays in the dynamics, we decompose the leading order solution into a "slow part" and a "fast part" respectively: We specify that X satisfies the "slow part" of the leading order Eq. (13) , i.e. F eq X ; Under the assumption that the "fast part" | X | 1 (see Appendix B ), we Taylor expand the F eq force term in (13) about X and retain only linear terms in X to obtain Provided that ω 2 = F eq (X ; λ 1 −β ) > 0 (as justified later in Appendix C ), the solutions are periodic; we denote the period by L = 2 π /ω, which will vary on the slow timescale as X varies. Integrating the equation from T = 0 to T = L then shows that Returning to the way we decomposed the solution in (14) , we see that X corresponds to the mean value of X 0 and varies on the slow timescale T . This evolution is captured by the variable A ( T ). The variable X describes the oscillations around this mean displacement that occur on the fast timescale T ; the property (17) guarantees that these oscillations do not influence the mean value if their amplitude is small. We now show that it is possible to obtain an evolution equation for A ( T ) without requiring detailed knowledge of X , using only the zero-mean property (17) . (While it is possible to obtain an analytical expression for X using the simplified Eq. (16) , we do not pursue this here, as this introduces a further unknown function of the slow timescale T -knowledge of X will be sufficient to determine when snap-through occurs and the associated snap-through time.)

First order problem
This represents a linear, inhomogeneous problem for X 1 . Setting the right-hand side to zero, we see that the homogeneous problem can be solved approximately whenever | X | 1 by taking X 1 = X 1 (T ) , as X 0 ≈ X (T ) in this case and so all T derivatives vanish. The Fredholm Alternative Theorem then implies that we determine A ( T ) from the solvability condition associated with the approximate homogeneous solution X 1 = X 1 (T ) ( Keener, 1988 ). To formulate this condition, we simply integrate the first order problem from T = 0 to T = L . We assume that for each fixed T , the solution X 1 is also a periodic function with period L ; this is reasonable, since X 1 is forced by the X 0 terms that have period L . It follows that all ∂ /∂ T terms vanish in the integration and we are left with Using the zero-mean property (17) , the first term can be evaluated as βλX / (1 − β ) . Eliminating A ( T ) for X using the relation (15) , we arrive at F eq X ; This equation is exactly the original ODE (8) in the limit De → ∞ , i.e. neglecting the terms associated with inertia. This is perhaps not surprising: when the zero-mean property (17) holds, the fast elastic oscillations 'cancel out' on the slow viscous timescale T and so do not affect the leading order dynamics. However, the above analysis does show that the correct initial condition is not the indentation displacement X (0+) = X ind , as might be expected. Instead, from (15) where the second equality follows from (7) . This correction arises from the initial transient around T = 0 in which inertia is always important; physically, the above equation states that the change in spring force in moving from X ind to X (0+) (when the SLS element is fully unrelaxed with effective stiffness λ/ (1 − β ) ) must balance the discontinuity in the indentation force. When viewed on the slow timescale, the mean value X then appears to change discontinuously from the indentation displacement X ind .
To check our multiple-scale analysis, we integrate the simplified ODE (18) numerically subject to the initial condition (19) . In Fig. 4 c solutions are superimposed (as black dotted curves) onto the trajectories obtained by integrating the full ODE (8) for the baseline parameter values. We see that the agreement is excellent, with the multiple-scale solution indeed capturing the average behaviour of the displacement during snap-through (see lower panel in Fig. 4 c). (The slight disagreement when the mean value changes rapidly on a timescale comparable to T is because the multiple-scale analysis is no longer applicable.) Fig. 4 c also shows that the initial value X (0+) may be much smaller than X ind depending on the indentation time T ind . We postpone a detailed analysis of X (0+) to section §3.2.3 below and Appendix C .

Snap-through dynamics
We have shown that while the amplitude of the oscillations is small compared to the mean displacement, the leading order behaviour is given by subject to the initial condition (19) . Since (20) represents a first-order autonomous ODE for X , the dynamics can be understood by considering the (X , d X / d T ) phase plane, for which the qualitative features are determined by the roots of F eq (X ; λ) = 0 and F eq (X ; λ 1 −β ) = 0 . For λ < λ fold , there are two distinct non-zero stationary points, which correspond to the stable and unstable roots of F eq (X ; λ) = 0 ; for λ > λ fold , the only stationary point is X = 0 . The roots of F eq (X ; λ 1 −β ) = 0 correspond to vertical asymptotes where | d X / d T | = ∞ . Denoting these roots by X ± , we find that These roots are real and distinct if and only if λ < 1 − β.
Different cases for the phase plane are illustrated in Fig. 7 . In Appendix C , we show that when λ < 1 − β (so that d X / d T may diverge), the initial value X (0+) must satisfy either X (0+) < X − or X (0+) > X + , i.e. X always starts outside the interval between the two vertical asymptotes. We have therefore lightly shaded the curve in this region in Fig. 7 . Fig. 7 shows how when λ > 1 − β, there are no vertical asymptotes and the trajectories smoothly approach the natural state; this explains the observation that pseudo-bistable behaviour does not occur for larger values of β, even though the snap-through time T snap ࣡ O (1) (discussed more in Appendix A ). We see that pseudo-bistability is only possible if λ < 1 − β: the asymptotes on the phase plane correspond to rapid snap-through events, which occur after a slow creeping motion provided X (0+) is sufficiently large (top right panel in Fig. 7 ). Since we are primarily interested in this pseudo-bistable regime, we restrict our attention to the case λ < 1 − β in the following analysis; in fact, we will restrict to β ≤ 1/2 so that we are safely in this regime in considering values of λ up to and slightly beyond λ fold 1 .
Focussing on λ < 1 − β, we deduce from Fig. 7 that there are three possibilities depending on where the solution starts in the phase plane: 1 As indicated by Fig. A.1 b in Appendix A , the truss may show different behaviour when λ is slightly below the critical value 1 − β. Further analysis shows that as the left asymptote X − on the (X , d X / d T ) phase plane approaches X − 1 (i.e. as λ 1 − β), the precise definition of T snap becomes important; for example, the mean value X (0+) may start to the left of both asymptotes so that the truss immediately snaps to near the natural configuration, but if X (0+) is close to 1 the truss does not cross X = 0 (our definition of T snap ) on the elastic timescale. By restricting to β ≤ 1/2, we are able to bypass these technicalities for values of λ in the interval of interest for pseudo-bistable behaviour.
• If X (0+) > X + and λ < λ fold , then the mean value starts to the right of both vertical asymptotes and approaches the stable non-zero root of F eq (X ; λ) = 0 . Because X + > 1 , the truss remains in an inverted position and there is no snapthrough 2 .
• If X (0+) > X + and λ > λ fold , then the mean value starts to the right of both vertical asymptotes. It then decreases to the vertical asymptote at X + where rapid snap-through occurs and inertial effects again become significant -because this approach to the asymptote occurs on the slow timescale T , the total time taken to snap-through is at least O (1). This regime corresponds to pseudo-bistable behaviour, in which the snap-through time is governed by the timescale of viscous relaxation.
• If X (0+) < X − , then the mean value starts to the left of the vertical asymptotes and smoothly decays to zero. Because X − < 1 , the truss immediately snaps back to near its natural configuration during the initial transient in which inertia is important. The amplitude of the elastic oscillations will therefore be large compared to the mean value in this case, and our assumption | X | 1 is no longer valid; nevertheless, we expect the truss to pass The final task is to determine when the initial value satisfies X (0+) > X + . This is not obvious because (19) implies that X (0+) is the root of a cubic polynomial, for which there may be multiple real solutions. The relevant solution can be found by analysing the phase portrait of Eq. (13) : setting A (T ) = 0 , this equation governs the elastic behaviour of the truss at very early times, and hence determines which root of (19) the solution approaches. When viewed on the slow timescale, this root corresponds to the relevant value of X (0+) . The full analysis is provided in Appendix C . The key result is that if 1 ≤ X ind < X fold , then X (0+) > X + if and only if F ind < 0; otherwise we have X (0+) < X − . Physically, this states that the indentation force needs to be adhesive for the truss to remain in an inverted state. This is intuitive: if the truss has to be 'pulled' upwards to the imposed indentation depth, it should move further downwards (increasing X ) when the indenter is released. Using the expression (7) for F ind , the condition F ind < 0 can be expressed as In the alternative case X fold ≤ X ind ≤ 2, the condition F ind < 0 turns out to no longer be relevant. Instead, in Appendix C we show that where Physically, this condition arises from bounding the amplitude of the elastic oscillations at very early times, so that these do not 'push' the truss sufficiently far from the inverted state and cause an immediate snap-back.
Combining this with the phase-plane discussion above, the predicted dynamical regimes are shown schematically in Fig. 8 . This explains how the qualitative features of the dynamics are very different in the two cases X ind < X fold and X ind ≥ X fold . When X ind < X fold , it may be shown that the boundary predicted by (21) reaches a vertical asymptote on the ( λ, T ind )-plane when λ < λ fold . For values T ind below the boundary we have X (0+) < X − , and the above discussion implies that the truss immediately snaps with T snap = O ( De −1 ) (shaded blue in Fig. 8 ). Above the boundary, X (0+) > X + and, because λ < λ fold here, the truss does not snap-through. Pseudo-bistable behaviour therefore cannot be obtained when X ind < X fold . Conversely, when X fold ≤ X ind ≤ 2, the boundary predicted by (22) reaches a vertical asymptote when λ > λ fold . Hence, there is a region where X (0+) > X + and λ > λ fold , in which pseudo-bistable behaviour occurs (shaded red in Fig. 8 ).
We deduce that T snap ࣡ O (1) precisely when To check the validity of the picture presented in Fig. 8 , we have superimposed the boundaries predicted by (21) -(22) (purple dashed curves, red dotted curves respectively) onto the numerical snap-through times in Figs. 5 -6 (and Fig. A.1 of Appendix A ). We observe that the agreement with the numerics is excellent when λ < 1 − β, despite the fact that the assumption | X | 1 made in the multiple-scale analysis is not formally valid throughout the range of values shown.  ( Brinkmeyer et al., 2013;2012;Madhukar et al., 2014;Santer, 2010;Urbach and Efrati, 2017;2018 ). We see 2 An additional case is possible in which the unstable root of F eq (X ; λ) = 0 is larger than the vertical asymptote at X + , so that two stationary points lie to the right of both asymptotes on the phase plane (unlike the top left panel of Fig. 7 , in which the unstable root lies in the lightly shaded region between the asymptotes). However, in focussing on values β ≤ 1/2, this regime is found to occur only for values of λ < λ fold extremely close to λ fold and so can generally be ignored.
In each panel, arrows indicate the direction of motion. The stable/unstable roots of F eq (X ; λ) = 0 are shown as filled/unfilled circles respectively. The vertical asymptotes at X ± are plotted as black dotted lines. The initial value X (0+) is determined by the transient around T = 0 in which inertia is important -this is the solution of the cubic (19) . Note that pseudo-bistability occurs for λ > λ fold , λ < 1 − β (top right). that pseudo-bistability occurs only in a narrow parameter range, near the threshold at which snap-through no longer occurs (i.e. λ = λ fold = 1 / 4 here) and the width of the pseudo-bistable region grows as the amount of stress relaxation increases (increasing T ind ). Pseudo-bistable behaviour is not obtained if T ind is too small, nor if the indentation depth X ind is below a critical value. In addition, the phase plane in Fig. 7 explains how the truss initially creeps in an inverted state before abruptly accelerating, leading to the difference in timescales that is characteristic of pseudo-bistable behaviour ( Brinkmeyer et al., 2013;2012 ). We emphasise that the analytical understanding of these features presented here is, to the best of our knowledge, new.
The importance of inertial effects immediately after the indenter is released has not been appreciated previously. Because this causes the displacement of the truss to change rapidly from the indentation displacement, the effective stiffness will also change rapidly from its value just before the indenter is released. This is in direct contrast to the viscoelastic models used by Santer (2010) and Brinkmeyer et al. (2013Brinkmeyer et al. ( , 2012 , which assume that (i) the stiffness reverses back to its fully unrelaxed value when the indenter is released, and (ii) there is no rapid change in the stiffness caused by the discontinuity in the applied indentation force. In Appendix D we show that when we make assumptions (i)-(ii) in the framework of the truss model, we obtain radically different predictions of when snap-through and pseudo-bistability occur. In the type of snap-through experiment considered here (a structure is allowed to relax in a specified displacement before being abruptly released), assumptions (i)-(ii) cannot therefore be derived from first principles starting from the constitutive law for a standard linear solid. Instead, we believe it is necessary to couple the stress within the structure to its displacement and account for inertial effects when the indenter is removed.

Snap-through time in the pseudo-bistable regime
Another key feature of the dynamics is that the snap-through time increases considerably as λ λ fold in the pseudobistable regime, becoming much larger than O (1) ( Fig. 5 b). This slowing down does not require X ind ≈ X fold (since it can be observed when X ind = 1 . 7 , Fig. 6 b). The phase plane in Fig. 7 (top right panel) suggests that this slowing down is due to a saddle-node ghost: when λ > λ fold the non-zero stationary point no longer exists, but as λ λ fold the trajectory passes increasingly close to the line d X / d T = 0 at X ≈ X fold . Because the velocity becomes very small but non-zero, this will lead to a slow passage through a bottleneck.
To analyse this slowing down in detail, we set where 0 < 1 and | χ| −1 / 2 ; here we anticipate an 1/2 scaling for the change in displacement during the bottleneck phase, which is the generic scaling for overdamped dynamics near a saddle-node bifurcation ( Strogatz, 2014 ). (We have also introduced a minus sign since we expect the displacement to decrease during snap-through.) Expanding the force terms appearing in the simplified ODE (20) , we obtain d χ which (up to numerical pre-factors) is the normal form for overdamped dynamics near a saddle-node bifurcation ( Strogatz, 2014 ); here the neglected terms are small compared to at least one retained term provided | χ| −1 / 2 . The solution is where The snap-through time is dominated by the time spent passing through the bottleneck, T b , which can be determined by finding the time at which χ first reaches O ( −1 / 2 ) ; after this point, we no longer have X ≈ X fold and so the truss is moving rapidly. Using the expansion tan x ∼ (π / 2 − x ) −1 as x → π /2, we have that It follows that there are different distinguished limits, depending on the size and sign of χ (0+) ; these are discussed in Appendix E . For example, for the baseline case X ind = X fold , we insert the expansions (23) into the initial condition (19) .
Upon neglecting terms quadratic in and e −T ind (e.g. with β = 1 / 2 , Fig. 5 b implies we have T ind ࣡ 4 in the pseudo-bistable regime, so that e −T ind 0 . 02 1 ), we obtain The above expression for the bottleneck duration then gives the leading order estimate The distinguished limits correspond to −1 / 2 e −T ind 1 and −1 / 2 e −T ind 1 , and we obtain   (27) approximates the numerically-computed snap-through times reasonably well, with the data indeed obeying an inverse square-root scaling law associated with an overdamped saddle-node ghost ( T snap ∝ −1 / 2 ) only when T ind is sufficiently large compared to log −1 / 2 . However, the O (1) error in (27) becomes significant if T ind ࣠ 5 or 10 −4 is only moderately small. While it is possible to obtain the O (1) correction analytically by integrating the ODE (20) directly, we do not compute this here. For deeper indentations X ind > X fold , we observe a similar picture, though the snap-through times are generally larger so that the O (1) correction is less significant; see Fig. 9

b.
A factor of β/ (1 − β ) = E 2 /E 1 (corresponding to the relaxation strength of the material) consistently appears in the expression (27) for the snap-through time. This means that the dimensional snapping time scales as η/ E 1 rather than the timescale of stress relaxation, η/ E 2 . To understand the origin of this timescale, we note that η/ E 1 corresponds to the creep timescale for a Kelvin-Voigt material, i.e. if we neglect the E 2 spring in the SLS element. Physically, this arises because the solution is close to the equilibrium at X = X fold as it passes through the bottleneck. In equilibrium, the E 2 spring (which is in series with the dashpot) will be relaxed at its natural length. Hence, sufficiently close to these solutions we do not expect this spring to play an important role compared to the E 1 spring, so the SLS element acts analogously to a Kelvin-Voigt element. Alternatively, we note that close to equilibrium the stress is nearly equal to its elastic value, = X. If we write = X + ξ for some small perturbation ξ and substitute this into Eq. (3) , at leading order we obtain [ β/ (1 − β )] ˙ X ∼ ξ (assuming ˙ ξ ξ due to the slow bottleneck timescale), which gives the same factor of β/ (1 − β ) ; this shows that the timescale η/ E 1 is the natural timescale over which the system undergoes bottleneck behaviour near an equilibrium 3 .

A continuous model system: Viscoelastic arch
Our analysis of the Mises truss indicates that in the limit of large Deborah number, only two types of snap-through are possible: the system either immediately snaps on the elastic timescale, or it is pseudo-bistable and first undergoes a slow creeping motion governed by the viscous timescale. This creeping motion may be very slow indeed, as we found that the snap-through time is subject to critical slowing down in this regime. As discussed in §2 , we may also consider the truss as a lumped model for a continuous viscoelastic structure (such as those used in morphing applications), when we identify λ with the analogous parameter that determines whether the continuous structure is bistable or monostable; for example the Föppl-von-Kármán number for complete spherical shells ( Knoche, 2014 ) and the analogous parameter for spherical caps ( Taffetani et al., 2018 ). Hence, we expect that our conclusions for the truss model may hold more generically and explain the features of pseudo-bistability observed in more complex systems.
To test this hypothesis, we now analyse the snap-through dynamics of a viscoelastic arch. This allows us to model a continuous structure without the additional complications that come with more complex structures such as spherical caps ( Gomez et al., 2016 ), where there is also the possibility of non-axisymmetric deformations during snap-through ( Seffen and Vidoli, 2016 ). Nevertheless, arches illustrate many features of snap-through that are present in more complex systems ( Harvey and Virgin, 2015 ), and approximate the behaviour of curved panels commonly used in engineering applications ( Wiebe, 2012 ). We consider a flat strip that is subjected to an end-shortening L in the horizontal direction, with one end clamped at an angle α to the horizontal while the other end is clamped horizontally. As shown schematically in Fig. 10 a, this The lower branch (solid blue curve) corresponds to the inverted shape, while the upper branch (solid red curve) corresponds to the natural shape. (The dotted branch is an unstable mode that is not observed experimentally.) (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.) causes the strip to buckle into one of two possible arch shapes -an 'inverted' shape (directed downwards) and a 'natural' shape (directed upwards). However, as a simple experiment illustrates (e.g. using an ordinary strip of plastic), the inverted shape needs a sufficiently large end-shortening L to be stable: for smaller values of L , the arch snaps from the inverted shape to the natural shape. Snap-through can also be initiated at a fixed end-shortening by increasing the clamp angle α. This allows us to study an analogous experiment to the truss system: the arch is indented to an inverted position and held fixed, allowing stresses to relax for a specified duration, and then the indenter is abruptly removed.
Previously, we have studied the snap-through dynamics of this system in the case of a purely elastic material, both for shallow arches using beam theory ( Gomez et al., 2017a ) and for deeper arches using the dynamic elastica equations ( Gomez, 2018 ). This analysis indicates that the bifurcation associated with snap-through is a saddle-node bifurcation, and so is qualitatively similar to the bifurcation observed in the Mises truss and spherical caps ( Taffetani et al., 2018 ). If instead the ends of the arch are held symmetrically (i.e. at equal angle α), the bifurcation changes type to a subcritical pitchfork. This alternative setup has been studied by Brinkmeyer et al. (2013) , using a combination of finite element simulations and experiments on arches composed of the rubbery polymer Sylgard 182. In their viscoelastic model, Brinkmeyer et al. (2013) use a Prony series expansion for the Young's modulus, which assumes that the modulus can be written as a sum of exponentially decaying modes; the coefficients in the sum and the relaxation timescales are fitted to experimental data for the relaxation of Sylgard 182. Here, we instead start from the constitutive law of a standard linear solid, and solve the equations of motion coupling the stress in the arch to its displacement. In this way, we will show that the behaviour of the arch is well captured by the truss model: the different dynamical regimes map directly onto those obtained for the truss, and the snap-through time in the pseudo-bistable regime obeys the same inverse square-root scaling law due to the ghost of the saddle-node bifurcation.

Theoretical formulation
The properties of the strip are its density ρ s , thickness h and natural length L . We model the arch shape using beam theory; this requires a small thickness h L , which guarantees that the strip remains in the limit of small strains (for a further discussion see Gomez, 2018 ), as well as a shallow arch shape (i.e. α 1 and L L ). Assuming the arch deforms only in the plane perpendicular to its width, we denote the transverse displacement by w ( x, t ), where x is the horizontal coordinate measured from the left end ( Fig. 10 a) and t is time. We model the indentation force as a transverse point force applied at the arch midpoint with magnitude f . Under the above assumptions, the transverse displacement satisfies the dynamic beam equation ( Howell et al., 2009 ) where P c ( t ) is the (unknown) compressive force applied to the arch and m is the bending moment (each per unit width).
Here we are also including viscous damping due to the environment, which is assumed to be linear in the velocity with constant coefficient Y (per unit width). In this framework the SLS constitutive law is applied using the same equation as for the truss, i.e. Eq. (1) in §2.1 , where now we interpret e as the axial strain field and σ as the axial stress within the strip. After relating these quantities to the displacement w and bending moment m in the small-slope approximation using standard relations in beam theory ( Wang and Chen, 2009 ), and integrating to solve for m , the bending term in (28) can be evaluated as where we assume initial data at t = t 0 . As t → ∞ , we recover the usual Euler-Bernoulli law for an elastic solid with Young's modulus E 1 , which here is written as where B = E 1 h 3 / 12 is the fully relaxed bending stiffness (per unit width). (This can be derived from (29) by using a Watsonlemma type argument to evaluate the integral in the final term to leading order as t → ∞ ; see Hinch, 1991 , for example.) With the Euler-Bernoulli law (30) , Eq. (28)  The boundary conditions at the clamped ends are (subscripts denoting differentiation) Neglecting the effects of extensibility, which is valid for shallow arches while h L ( Pandey et al., 2014 ), the imposed endshortening may be approximated as Eqs. (28) -(29) , the boundary conditions (31) and constraint (32) , together with appropriate initial conditions, then fully specify the problem.

Non-dimensionalisation
To make the problem dimensionless, it is convenient to scale the horizontal coordinate with the natural length L of the strip. The end-shortening constraint (32) provides a natural vertical lengthscale w ∼ ( L L ) 1/2 ; using the Euler-Bernoulli law (30) , this is associated with a typical bending moment m ∼ B ( L L ) 1/2 / L 2 , which motivates introducing Here we again scale time with the timescale of stress relaxation, t ∼ η/ E 2 . In terms of these variables, the beam Eq. (28) be- where we have introduced The Deborah number here measures the ratio of the viscous timescale η/ E 2 to the timescale of undamped elastic oscillations, t * ∼ ρ s hL 4 /B , so is analogous to the Deborah number defined in §2.2 for the truss. The other parameters correspond to the dimensionless values of the damping coefficient, compressive force and indentation force, respectively. Inserting the above rescalings into the bending term (29) , we obtain where β is defined as in §2.2 and T 0 = t 0 / (η/E 2 ) . The clamped boundary conditions (31) become where we have introduced the normalised clamp angle Finally, the imposed end-shortening (32) becomes

Steady bifurcation behaviour
In the absence of any indentation force, F = 0 , equilibrium solutions obey the steady beam equation together with the clamped boundary conditions (35) and end-shortening constraint (36) . The equilibrium behaviour of the arch is therefore entirely characterised by the geometric parameter μ (since τ is unknown and is determined as part of the solution). It is possible to solve the equilibrium problem analytically ( Gomez et al., 2017a ), which indicates that for 0 < μ < μ fold ≈ 1.7818, both inverted and natural equilibrium shapes exist and are linearly stable. The critical value μ = μ fold corresponds to a saddle-node bifurcation where the inverted shape intersects a higher, unstable mode and disappears; for later reference we write W fold ( X ) for the equilibrium shape at the bifurcation point, which has midpoint height W fold (1 / 2) ≈ −0 . 3476 . For values μ > μ fold , only the natural shape exists and is stable. The bifurcation diagram is shown in Fig. 10 b. The parameter μ therefore plays an analogous role to the stiffness parameter λ in the truss model (compare Fig. 10 b to Fig. 3

Snap-through dynamics
In a numerical snap-through experiment, we suppose that the arch is initially fully relaxed in the natural shape, labelled W nat ( X ), for each value of μ. For −T ind < T < 0 we then rapidly indent the arch to an inverted position with midpoint displacement W mid < 0. Typically we specify W mid = W fold (1 / 2) , i.e. the midpoint displacement at the saddle-node bifurcation, which is the analogue of the baseline value used for the truss system. With the imposed midpoint displacement, we integrate (33) -(36) with initial conditions For the release problem, i.e. T > 0, we instead impose zero indentation force F = 0 during the integration. Because W is continuous across T = 0 , we have that M is also continuous, giving the initial conditions We solve the dynamic equations using the method of lines ( Morton and Mayers, 2005 ), i.e. we discretise using finite differences in space to obtain a system of ODEs in time; details of the numerical methods are provided in Appendix F . The ODEs are integrated numerically in matlab using the ode45 routine (relative and absolute error tolerances 10 −8 , maximum time step 10 −4 ). We have verified second-order accuracy in the convergence of our scheme, and also that the numerical drift in the end-shortening constraint remains on the order of the integration tolerances. Because we are interested in performing a large number of simulations as μ and T ind are varied, we discretise using N = 50 grid points for all simulations reported in this section; we have checked that this yields quantitatively similar results compared to using a larger number of grid points, e.g. N = 100 , but requires much less computing time. We also specify a dimensionless damping coefficient υ = 0 . 5 , which provides sufficient numerical damping (needed due to the fast motions during the indentation stage) while ensuring the arch motions remain underdamped. Typical trajectories of the arch midpoint, W (1/2, T ), during the release stage are shown in Figs. 11 a-d. Here we have specified De = 10 and β = 0 . 1 ; this value of β is an approximate value measured experimentally by Urbach and Efrati (2017) for silicone rubber shells. (While pseudo-bistable behaviour is obtained with β = 0 . 5 for the truss, we will show that smaller values are required for the arch system.) Fig. 11 shows that the arch exhibits similar dynamical regimes to the truss model, with μ playing the role of the stiffness parameter λ: for μ (1 − β ) μ fold ( ≈ 1.6 with β = 0 . 1 ) the arch never snaps ( Fig. 11 a); for (1 − β ) μ fold μ μ fold the arch snaps if T ind is sufficiently small ( Fig. 11 b); and for μ ࣡ μ fold the arch appears to snap for any indentation time ( Figs. 11 c-d). Pseudo-bistable behaviour is obtained when 0 < μ − μ fold 1 and T ind is sufficiently large ( Fig. 11 c). Note that the oscillations are evidently underdamped, but do not persist long (compared to the truss, Fig. 4 ) due to the presence of external damping.
The different regimes just discussed are also evident when we analyse the snap-through times on the ( μ, T ind )-plane; the baseline case W mid = W fold (1 / 2) ≈ −0 . 3476 and β = 0 . 1 is shown in Fig. 12 a. (Because the displacement of the natural shape depends on μ, we define T snap to be the time when the arch midpoint first crosses W = 0 .) We observe very similar features to the analogous plot for the truss system (compare Fig. 12 a to Figs. 5 a-b in  §3 ), in that (i) the critical value of T ind at which snap-through no longer occurs depends nonlinearly on μ, and approaches a finite value as μ μ fold ; (ii) pseudo-bistability occurs only in a narrow region where 0 < μ − μ fold 1 and T ind is larger than the critical value obtained at μ = μ fold ; and (iii) slowing down occurs in the pseudo-bistable regime as μ μ fold .
The dependence on the indentation displacement W mid and relaxation parameter β is also similar to that in the truss model. For shallower indentation depths W mid > W fold (1/2) (corresponding to X ind < X fold for the truss), the boundary at which snap-through no longer occurs is shifted to the left of the line μ = μ fold and pseudo-bistability is never obtained; see Fig. 12 b. For a deeper indentation W mid < W fold (1/2) (corresponding to X ind > X fold for the truss), the boundary may shift to the right so that the pseudo-bistable region is enlarged ( Fig. 12 c). When β is increased, this picture breaks down as the snap-through time is O (1) throughout a large portion of the ( μ, T ind )-plane (see Fig. 12 d). Similar to the truss ( Fig. A.1 b in Appendix A ), an analysis of the trajectories in this regime confirms that this behaviour is different to pseudo-bistability, being instead a consequence of the extremely slow creep timescale associated with larger values of β. We note that because this behaviour is obtained with β = 0 . 5 here, as opposed to larger values for the truss, the value of β does not map directly between the two systems.
We have also simulated the snap-through times in the pseudo-bistable regime, setting μ = μ fold + μ where 0 < μ 1; data for two different indentation depths are shown in Figs. 13 a-b. These confirm that the expected inverse square-root scaling T snap ∝ μ −1 / 2 is obtained as μ → 0 when T ind is sufficiently large. While the data deviate significantly from the scaling for moderately small values μ 10 −3 , this is similar to the behaviour we have seen in the truss model, in which the error of the asymptotic prediction becomes significant if is not too small (compare to Figs. 9 a-b). Moreover, as with the truss model, this error is less significant for a deeper indentation depth ( Fig. 13 b), since the snap-through time is generally larger in this case compared to the baseline value, W mid = W fold (1 / 2) ( Fig. 13 a).
In summary, the truss provides an excellent lumped model of the continuous arch system, providing qualitatively correct predictions of the different dynamical regimes, the features of pseudo-bistability and the scaling laws for the snap-through time.

Data comparison
To further examine whether the inverse square-root scaling law for the snap-through time in the pseudo-bistable regime holds more generically, in this final section we examine numerical and experimental data for the snap-through times of viscoelastic shells and arches reported in the literature. As was the case in our numerical experiments, each structure is held in an inverted state for a duration t ind , before being instantaneously released. The snap-through time is measured to be the time taken between release and when the structure rapidly accelerates towards its natural state. We focus on results in which the indentation time t ind is fixed, while the analogue of the bifurcation parameter λ or μ is varied between each snapping experiment. In all cases examined the structure exhibits pseudo-bistability, undergoing a slow creeping mo- tion followed by a rapid snap-through, so that the snap-through times are easily measured. Where data is only available graphically, we have extracted the values using the WebPlotDigitizer (arohatgi.info/WebPlotDigitizer).
A summary of the conditions for each data set is provided in Table 1 . Here we have separated the data so that only a single parameter is varying within each data set (corresponding to a particular row in the table), and we have provided the relevant parameter values. These are the shell/arch thickness h , relevant horizontal lengthscale l (defined to be the base diameter of the shell/natural length of the arch), fully relaxed Young's modulus E , Poisson ratio ν, material density ρ s , viscous timescale [ t ] vis , elastic timescale t * , and the indentation time t ind . (Where a parameter varies within a data set, the range of values is provided, and we use the average to compute t * .) The simulation data reported by Brinkmeyer et al. (2013Brinkmeyer et al. ( , 2012) assume a Prony series expansion for the Young's modulus, so there is no single viscous timescale; we estimate [ t ] vis by the timescale that appears in the dominant term of the Prony series (i.e. the term with the largest coefficient). To estimate the elastic timescale, we balance inertial forces and bending forces for a shell/arch to obtain ( Gomez et al., 2017a;Ventsel and Krauthammer, 2001 )   where l is the horizontal lengthscale defined above, and B is the bending stiffness (in particular B = Eh 3 / [12(1 − ν 2 )] for a shell, and B = Eh 3 / 12 for an arch). Table 1 shows that for all data sets, the viscous timescale [ t ] vis is an order of magnitude larger than the elastic timescale t * , so these systems are effectively in the large Deborah number limit. The indentation times t ind are also much larger than [ t ] vis , so that the dimensionless indentation times T ind = t ind / [ t] vis are large. Hence, the pseudo-bistability observed in these systems occurs in an analogous parameter range to that in our truss and arch models.
For each data set it is found that as the analogue of λ varies, the snap-through time increases rapidly and appears to diverge near a critical value. Beyond this transition no snap-through occurs. This transition therefore appears to be the saddle-node bifurcation at which the inverted arch/shell becomes bistable 4 . We use the critical value that is reported to compute the normalised distance to the bifurcation, which we denote by eff ; for example, if the thickness h is varied and h c is the critical value when snap-through no longer occurs, then we define eff = | h − h c | / | h c | and similarly when other parameters are varied. The dimensional snap-through times are plotted as a function of eff on logarithmic axes in Fig. 14 a. Here we observe the characteristic signs of critical slowing down: as eff → 0 the snap-through time increases systematically, varying by over two orders of magnitude within a very narrow range of eff .
The key observation is that the data of Brinkmeyer et al. (2013) (yellow symbols) and Urbach and Efrati (2017) (orange symbols) are approximately consistent with an inverse square-root scaling law, i.e. t snap ∝ −1 / 2 eff as eff → 0. To be more quantitative, we fit each data set to a power law of the form t snap ∝ γ eff using least squares. We restrict the fit to values eff ≤ 10 −2 , since this is the range where the inverse square-root scaling is approximately observed for the truss and arch systems (see Figs. 9 and 13 ). (One data set of Brinkmeyer et al. (2012) , blue squares, has no values of eff ≤ 10 −2 so here we instead fit the six points closest to bifurcation.) The best-fit exponents are provided in Table 1 . This confirms that for the data reported by Brinkmeyer et al. (2013) and Urbach and Efrati (2017) the fitted values are dispersed around γ = −0 . 5 , lying in the range γ ∈ (−0 . 599 , −0 . 316) . We note that such dispersion is expected due to the small parameter range over which the data exhibit a power law: the fitted value of γ is sensitive to the precise range of values of eff used for fitting.
In addition, this dispersion may be due to sensitivity to the precise value of the bifurcation point used to calculate eff : a small error in this value (e.g. a rounding error in the reported value) introduces shifts in the values of eff , which can cause large variations when plotted on logarithmic axes. To eliminate this second type of sensitivity we plot t −2 snap as a function of eff on linear axes, where a linear relationship indicates that the inverse square-root scaling is obeyed. This plot is shown in Fig. 14 b, focussing on values eff 10 −2 where we expect to observe the inverse square-root scaling; to give a clearer plot here we have re-scaled time so that all data sets approximately pass through (10 −2 , 1) . Fig. 14 b shows that the data of Brinkmeyer et al. (2013) and Urbach and Efrati (2017) become approximately linear as eff → 0: for each data set the dotted line is the best-fit (least squares) line over the six data points that are closest to the bifurcation point (for the experimental data of Brinkmeyer et al. (2013) , yellow downward-pointing triangles, three data points lie in the range eff ≤ 10 −2 so we fit only these). Nevertheless, in the absence of more data points in this parameter range, it is not possible to state conclusively that the inverse square-root scaling law holds for these systems.
In contrast, the snap-through times reported by Brinkmeyer et al. (2012) (blue symbols) consistently do not follow the inverse square-root scaling: the best-fit exponents γ −1 and a linear relationship is not observed in Fig. 14 b (for clarity we do not plot the best-fit lines on Fig. 14 b for these data). Nevertheless, the qualitative features of the pseudo-bistable regime, including the sensitivity of the snap-through time to changes in eff , is well captured by the truss model.

Discussion and conclusions
In this paper, we have analysed the dynamics of snap-through when viscoelastic effects are present. We re-emphasise that such effects are fundamentally different to external damping such as viscous drag: viscoelasticity modifies the bistability characteristics and hence can change not only how, but also when snap-through occurs. Moreover, in an experiment where we first indent a structure to a particular configuration, the resulting dynamics depend on the history of stress relaxation during the indentation phase. Previous approaches have dealt with this complexity by modelling the structure as being elastic with an effective stiffness that evolves according to a Prony series ansatz ( Brinkmeyer et al., 2013;2012;Santer, 2010 ). Here, we have presented an alternative approach that derives the equations of motion from first principles using the constitutive law of a standard linear solid. This enables us to capture both stress relaxation and creep phenomena without making any additional assumptions regarding the dynamic behaviour.
To gain analytical insight we first studied a modified form of the Mises truss, a simple and commonly used model system that exhibits bistability and snap-through ( Brinkmeyer et al., 2013;Krylov et al., 2008;Panovko and Gubanova, 1987 ). By introducing an additional vertical spring and a point mass in our formulation, the truss becomes a more realistic lumped model for more complex structures such as spherical shells and arches. Using a small-angle approximation, we were able to reduce the number of dimensionless parameters in our problem to five. These are the Deborah number De, measuring the relative sizes of the timescales of viscous relaxation to the intrinsic elastic oscillations; the relaxation parameter β, which measures the ability of the structure to relax its stress; the relative stiffness λ, which acts as a bifurcation parameter and determines the bistability characteristics of the truss; and the details of the indentation stage are specified by the indentation displacement X ind and duration T ind . Regarding the truss as a lumped model, we then expect that analogous parameters will control the dynamics in more complex viscoelastic structures; for example, λ may be compared to the Föppl-von-Kármán number for spherical shells.
We focussed on the dynamics when De is large. Using direct numerical solutions, we showed that the onset of snapthrough cannot be inferred by whether or not the truss is effectively bistable at the moment the indenter is released. Instead, we turned to a detailed asymptotic analysis of the snap-through dynamics using the method of multiple scales. This analysis showed that the leading-order dynamics generally obey the equations of motion when we neglect the terms associated with inertia, as expected. However, immediately after the indenter is released, inertial effects become important, as the displacement moves rapidly to a new value to balance the jump in the applied force. It is this purely elastic behaviour at early times that determines whether the truss then creeps in an inverted state or immediately jumps back to near its natural configuration. In this way, we were able to build up a complete picture of the different dynamical regimes and determine precisely when pseudo-bistable behaviour is obtained ( Fig. 8 ). Our analysis describes many features of pseudo-bistability that have been reported previously in experiments and numerical simulations, such as why a minimum indentation depth X ind and duration T ind are needed to obtain any creep behaviour. We then analysed an analogous indentation experiment performed on a pre-buckled arch, a simple and common prototype of a continuous viscoelastic structure. By solving the dynamic equations numerically, we were able to confirm that our conclusions for the truss model are qualitatively accurate for the arch, with the normalised clamp angle μ playing the role of the stiffness parameter λ. The features of pseudo-bistability that we predict are also readily observed in a commercially available popper toy: this needs to be turned sufficiently far inside-out, and held for a few seconds, in order to not immediately jump upwards when placed on a surface.
In the pseudo-bistable regime, the truss undergoes a creeping motion until a rapid snap-back occurs. In our leadingorder description of the dynamics, this snapping event is associated with an infinite velocity, implying that inertial effects must become important again. This is very similar to the analysis of creep buckling, in which an infinite velocity is used as a criterion to determine the onset of instability ( Hayman, 1978 ). Many studies on creep buckling consider only force-control situations, such as a force of constant magnitude suddenly applied to a structure. This means that the Kelvin-Voigt model is often sufficient to study the dynamics (see Nachbar and Huang, 1967 , for example). The snap-back considered here is considerably more complicated: due to the initial indentation stage, both the history of stress relaxation and inertial effects at early times must be accounted for. Pseudo-bistable snap-through then requires a combination of both stress relaxation (during indentation) and creep effects. For this reason, we must use a standard linear solid constitutive law to capture both of these effects, rather than a Kelvin-Voigt or Maxwell model. Pseudo-bistability is fundamentally different to the slow snap-through studied by Gomez et al. (2017a) , who considered a purely elastic system in which the effects of viscoelasticity and external damping were negligible. In that case, snapthrough occurred on the elastic timescale, and only became slow because of the phenomenon of critical slowing down: this effectively introduces a dimensionless pre-factor that multiplies the timescale of snapping, which can grow much larger than O (1) very close to the snapping transition. If external damping instead dominates inertial forces, a similar slowing down occurs though with the elastic timescale replaced by the damping timescale; see Gomez (2018) . Pseudo-bistable snapthrough is therefore only possible in systems with internal damping, i.e. material viscoelasticity, to provide the required stress relaxation and creep.
As well as being characterised by overdamped dynamics, a key feature of the pseudo-bistable regime is that it occurs in a narrow parameter range at the transition between bistability and monostability. This corresponds to the saddle-node bifurcation at λ = λ fold (for the truss) and μ = μ fold (for the arch). As a direct consequence, the snap-through dynamics are susceptible to critical slowing down. Provided T ind is sufficiently large, we showed that for both systems, the snapthrough time t snap inherits an inverse square-root scaling law, i.e. we have t snap ∝ (η/E 2 ) −1 / 2 where η/ E 2 is the viscous timescale and is the normalised distance to the bifurcation in parameter space (for the truss we in fact showed that t snap ∝ [ η/E 1 ] −1 / 2 , i.e. the relevant viscous timescale is η/ E 1 not η/ E 2 , though as discussed at the end of §3.2.4 both timescales are of the same order because pseudo-bistability requires that E 2 /E 1 = β/ (1 − β ) = O (1) ). The −1 / 2 scaling arises because the dynamics here are overdamped so that the time derivative in the normal form (24) is first order, in contrast to underdamped systems in which the dynamics are second order in time and the bottleneck duration scales as ∼ −1 / 4 ( Gomez et al., 2017a ).
In analysing data from experiments and numerical simulations reported in the literature, we found that this slowing down explains the sensitivity of the snap-through time observed in the pseudo-bistable regime. While there is some evidence of the inverse square-root scaling law it is not conclusive ( Fig. 14 ). In future work, it would be interesting to perform further experiments and determine whether the inverse square-root scaling indeed holds more generically. We also note that this sensitivity is responsible for large quantitative errors between finite element simulations and experiments despite excellent qualitative agreement ( Brinkmeyer et al., 2013 ), and in morphing applications would mean that parameters need to be precisely tuned to obtain the desired response. In such applications, our analytical expression could help to resolve the issue, as it provides a simple power law with which to control and calibrate the dynamic response: once the coefficient in the power law is determined (e.g. by fitting experimental data), it is possible to make further predictions without the need for detailed simulations.
Finally, while the details of many practical morphing applications are likely to be more complicated than the relatively simple examples considered here, we believe that our analytical insight will be useful in guiding designers to the appropriate parameter regimes. Our analysis also reveals important features of the dynamics that cannot be neglected, which we expect to also be important in more complex structures: as well as the importance of inertia, we have shown that the effective stiffness does not smoothly reverse back to its fully unrelaxed value when the indenter is removed, as has been assumed in previous numerical studies. This means that any viscoelastic model used in finite element simulations should couple the stress to the deformation without making additional assumptions.
In this appendix we discuss the influence of the parameter β (defined in Eq. (4) ) on the snap-through dynamics of the truss. The computed snap-through times for different values of β are shown in Figs. A.1 a-b for the case X ind = X fold . When β < 1/2, we observe a similar picture to that for the baseline value β = 1 / 2 considered in the main text, though the size of the pseudo-bistable region shrinks considerably in the purely elastic limit β → 0. For example to obtain any pseudo-bistable behaviour when β = 0 . 3 , it is necessary to take values (λ − λ fold ) 10 −3 and T ind ࣡ 5 (see the inset of Fig. A.1 a), compared to values (λ − λ fold ) 6 × 10 −3 and T ind ࣡ 4 for the baseline case β = 1 / 2 ( Fig. 5 b).
When β > 1/2, however, we observe very different behaviour: the region of parameter space where T snap ࣡ O (1) is much larger, extending to values λ < λ fold ( Fig. A.1 b). Crucially, this slowing down is qualitatively different to pseudo-bistability, as may be observed from the typical trajectories in this regime -see Fig. A.2 . These show that for sufficiently large T ind , the truss simply relaxes back to the natural shape without rapidly accelerating (yellow dashed curve, purple dashed-dotted curve in Fig. A.2 ). In particular, a slow creeping motion followed by a rapid snap-through event, indicative of pseudo-bistable be- haviour ( Brinkmeyer et al., 2013;2012 ), is not observed. This is a consequence of the creep timescale becoming unbounded in the Maxwell limit β → 1 (recall the discussion in §2.2 ), with the truss displaying more fluid-like behaviour.
Appendix B. Assumption of a small "fast part", | X | 1 In this appendix we discuss the "fast part" X of the leading order solution in the multiple-scales analysis of §3.2 ,defined in Eq. (14) . In terms of X , the leading order Eq. (13) becomes where we have expanded the F eq force term in (13) about X 0 = X (higher-order terms vanish because F eq is a cubic polynomial). The initial conditions (12) imply that Due to the nonlinear terms in (B.1) , it is difficult to make analytical progress, though we note that by multiplying by ∂ X /∂ T and integrating twice, it is possible to obtain an implicit equation for X (up to quadrature). However, each integration introduces an unknown function of the slow timescale T (in addition to the unknown function A ( T ) in Eq. (13) ), which require solvability conditions to be determined. We therefore make the simplifying assumption that | X | 1 , so that (B.1) is as given in the main text.

Appendix C. Determining X (0+)
In this appendix we consider the solution X (0+) of the cubic Eq. (19) , which we reproduce here: This provides the initial value of the "slow part" of the leading order solution in the multiple-scales analysis, and so determines whether the truss initially creeps in an inverted position or immediately snaps to near its natural state. We focus on the case λ < 1 − β, so that the force-displacement curve F eq (X ; λ 1 −β ) has distinct real turning points at X ± . We now discuss the cases λ < (1 − β ) λ fold and (1 − β ) λ fold < λ < 1 − β separately.
In this case the turning point X + on the force-displacement curve F eq (X ; λ 1 −β ) lies below the horizontal axis, i.e. F eq (X + ; λ 1 −β ) < 0 . Furthermore, restricting to β ≤ 1/2 and X ind ≤ 2, it is possible to show that i.e. the line of height βλX ind 1 −β 1 − e −T ind lies below the turning point at X − . There are therefore three distinct real roots of the cubic Eq. (19) . This is illustrated in the left panel of Fig. C.1 , which highlights the roots as red circles. However, it is not immediately clear which root is the relevant solution for X (0+) .
To determine the relevant root, we return to Eq. (13) , i.e. the leading order problem in the multiple-scale analysis. Recall that in our multiple-scale analysis, we first re-scaled the ODE in terms of the fast elastic timescale T . Setting A (T ) = 0 , Eq. (13) then governs the dynamics of the truss at very early times, before viscous relaxation of the SLS element becomes important. Defining V 0 = ∂ X 0 /∂ T , this equation can be written as the first-order system The significance here comes from the fact that the critical points of this system are precisely the solutions of the cubic (19) , and hence correspond to the possible values of X (0+) . The relevant value is then determined by which critical point the solution oscillates around on the phase plane, when we follow the trajectory emerging from (X 0 , V 0 ) = (X ind , 0) -these oscillations correspond to the "fast part" of the solution, X , as opposed to the slowly varying mean. For later reference, we also note that (C.1) is Hamiltonian with conserved energy 1 2

(C.2)
A typical phase plane of (C.1) in the case λ < (1 − β ) λ fold is shown in the right panel of Fig. C.1 . Here the three real roots of Eq. (19) give rise to three critical points. By considering the Jacobian of (C.1) , we find that the two roots where F eq (X ; λ 1 −β ) > 0 correspond to centres, while the intermediate root where F eq (X ; λ 1 −β ) < 0 is associated with a saddle point. Fig. C.1 also shows that there are two homoclinic orbits that emerge from the saddle point (black solid curves) that act as separatrices: all trajectories that oscillate around the left centre are enclosed in the homoclinic orbit to the left of the saddle point, while all trajectories that oscillate around the right centre are enclosed in the right orbit. We therefore expect that X (0+) corresponds to one of the centres rather than the saddle point, depending on whether the solution starts to the left or the right of the saddle point on the phase plane. Moreover, because F eq (X ; λ 1 −β ) > 0 at the centres, and F eq (X ; λ 1 −β ) < 0 when X ∈ (X − , X + ) (as seen from the force-displacement curve), we must have X (0+) < X − or X (0+) > X + . (In addition, because we expect the solution to oscillate around these centres as T increases and A = 0, this verifies our earlier assumption that F eq (X ; λ 1 −β ) > 0 made immediately after Eq. (16) .) Suppose that 1 ≤ X ind < X fold . From the force-displacement curve in Fig. C.1 , we see that X ind lies to the right of the saddle point if and only if F eq X ind ; is monotonically decreasing on the diagram). Using the expression (7) for F ind , this is precisely the statement When this is satisfied the relevant solution is the right centre with X (0+) > X + ; otherwise X (0+) < X − .
Suppose instead that X ind ≥ X fold . In this case, X ind always lies to the right of the saddle point (the force-displacement curve illustrates how the intermediate root is smaller than X fold when λ < (1 − β ) λ fold ). However, it is also possible that X ind is large enough to fall outside the homoclinic orbit to the right of the saddle point. If this occurs, the phase plane in Fig. C.1 shows how the amplitude of the oscillations becomes very large, with the trajectory enclosing both centre points. In particular, the trajectory crosses the horizontal axis again near the origin, and so the truss will immediately snap in this regime. If this occurs, we consider the relevant solution to be the left centre, i.e. X (0+) < X − . Now imagine that λ and X ind are fixed, while the indentation time T ind is varied. As T ind increases, the value of βλX ind 1 −β 1 − e −T ind increases, and so both the right centre and the homoclinic orbit are shifted further to the right on the phase plane. There will be a critical value of T ind for which the initial point ( X ind , 0) lies exactly where the homoclinic orbit crosses the horizontal axis (highlighted as a yellow square on Fig. C.1 ). Only when T ind is larger than this value does the initial point fall inside the separatrix and we have X (0+) > X + . This change in behaviour is an instance of a homoclinic bifurcation ( Strogatz, 2014 ), which has been observed in other dynamic snap-through ( Nachbar and Huang, 1967 ) and pull-in instabilities ( Krylov, 2007 ).
To determine this critical value, we note that there is a discontinuous change in where the trajectory starting from ( X ind , 0) later crosses the horizontal axis. Setting ∂ X 0 /∂ T = 0 in (C.2) , the value of X 0 at this point satisfies In this regime, Eq. (19) may have three distinct real roots (labelled as case 1) or a single real root (case 2). satisfies (C.3) , otherwise X (0+) < X − ; while if X ind ≥ X fold , we have X (0+) > X + if and only if T ind satisfies (C.4) , otherwise X (0+) < X − .

Appendix D. Assumption of an evolving stiffness
In the approach of Santer (2010) , later adopted by Brinkmeyer et al. (2013Brinkmeyer et al. ( , 2012 , the viscoelastic response is modelled as an elastic structure with an effective stiffness that changes in time. In particular, the evolution is specified by a Prony series, which assumes that the stiffness can be expressed as a sum of decaying exponential functions. For example, if E ( t ) is the Young's modulus at time t , and E 0 = E(0) is the initial modulus, this can be written as ( Brinkmeyer et al., 2012 ) E(t ) = E 0 1 − where N ≥ 1 is an integer; the coefficients K j and timescales τ j are specified parameters that can be fitted to experimental data from relaxation tests. In the case of a step increase in strain applied at t = 0 , this model is physically equivalent to a superposition of SLS elements, with each term corresponding to the stress relaxation of a particular element ( Brinkmeyer et al., 2012;Kim et al., 2010 ); the value E 0 corresponds to the fully unrelaxed modulus, i.e. just after the strain is applied. As t → ∞ , the modulus E ( t ) decays to the fully relaxed value To apply this model when the indenter is released, Santer (2010) and Brinkmeyer et al. (2013Brinkmeyer et al. ( , 2012 assume that the evolution during recovery is the reverse of the behaviour during indentation, and that there is no jump in the value of the stiffness. The stiffness is also allowed to fully relax before the indenter is released. With t = 0 now denoting the point when the indenter is released, this corresponds to setting K j e −t/τ j , for t > 0. Thus E(t = 0+) = E ∞ , and E ( t ) decays to the fully unrelaxed value E 0 as t → ∞ . If instead the stiffness is only allowed to relax for a time duration t ind before the indenter is released, as in our truss model, this modifies to since when t = 0+ this corresponds to (D.1) evaluated at t = t ind .
To understand this reversibility assumption within the framework of our truss model, we recall from §2.4 that in response to an indentation displacement X ind suddenly applied at T = −T ind , the dimensionless stress is Since X is the strain in the SLS element, the effective stiffness is therefore This corresponds to the Prony series (D.1) when we identify for some pressure scale [ σ ]. Hence, when the indenter is released, the above assumption (D.2) would be equivalent to specifying for T > 0 (Here we instead identify t/τ 1 = T .) Setting F = 0 in the momentum Eq. (5) , and substituting the effective stiffness above, the trajectory X ( T ) would then obey De −2 d 2 X d T 2 = −F eq X (T ) ; λ rel eff (T ) , where the effective value of λ is which increases from λ rel eff (0+) = λ[1 + β 1 −β e −T ind ] to λ rel eff (∞ ) = λ/ (1 − β ) during release. The initial conditions are X (0+) = X ind , ˙ X (0+) = 0 .

(D.4)
To understand how the solution of (D.3) differs to the model we consider in the main text when De 1, we neglect the inertia term as a first approximation to obtain F eq X (T ) ; λ rel eff (T ) = 0 , Recall that indentation corresponds to rotating the force-displacement curve in Fig. 3 clockwise as stress relaxation occurs. The assumption that the response during recovery is the reverse of relaxation then suggests that the graph simply rotates back anticlockwise as soon as the indenter is released. For each value of λ, the displacement is found as the value of X at which the force-displacement curve crosses the horizontal axis. Hence, the different dynamical regimes follow from those discussed during indentation in §2.4 . For λ < (1 − β ) λ fold (so λ rel eff (∞ ) < λ fold ), the truss is always bistable during indentation, and so remains bistable during recovery -we generally expect the truss not to snap-through. For λ > λ fold , the truss is monostable when the indenter is released (since λ rel eff (0+) > λ fold ), and so remains monostable during recovery -to satisfy F eq = 0 the truss must immediately jump to X = 0 and the snap-through time is governed by inertia. For (1 − β ) λ fold < λ < λ fold , the truss is temporarily bistable when the indenter is released if and only if T ind satisfies the condition (6) , which is equivalent to λ rel eff (0+) < λ fold ; in this case, rapid snap-through occurs as soon as the anticlockwise rotation is enough to put the turning point on the force-displacement curve above the line F eq = 0 . This regime corresponds to pseudo-bistable behaviour in which the snap-through time is O (1). If (6) is not satisfied, the truss is initially monostable and so immediately snaps.
Comparing the above picture with Fig. 8 of the main text, we conclude that the assumption of reversibility leads to very different behaviour to that derived without this assumption: both the regions where snap-through and pseudo-bistability occur differ significantly between the two models. This conclusion also holds when we account for the effects of inertia and directly solve the ODE (D.3) with initial conditions (D.4) numerically; see Fig. D.1 , which plots the computed snap-through times on the ( λ, T ind )-plane for the baseline case X ind = X fold and β = 1 / 2 . This shows that the region where snap-through occurs is approximately in agreement with the above analysis (contrast this to Fig. 8 ), though the boundary is shifted slightly to the left of the line λ = (1 − β ) λ fold due to the de-stabilising effects of inertia. displacement W ( X, T ). To avoid this, it is helpful to eliminate the bending moment in favour of the variable I(X, T ) = e T (1 − β ) ∂ 2 M ∂X 2 + ∂ 4 W ∂X 4 .
The bending term (34)  In our numerical scheme, we approximate the spatial derivatives appearing in (F.1) -(F.2) and (35) -(36) using centered differences with second-order accuracy, and we apply the trapezium rule to approximate the integral in the end-shortening (36) . During the release stage when the indentation force is zero, we apply the scheme at all interior grid points in the numerical mesh, using ghost points to evaluate the derivatives near the mesh boundaries without losing accuracy. During the indentation stage, we instead apply the scheme at all interior points away from the arch midpoint; the equation at the midpoint is replaced with where W mid < 0 and κ > 0 are prescribed constants. This form imposes a displacement that smoothly decreases from the initial value W nat (1/2), associated with the natural shape, and approaches the value W mid corresponding to an inverted position. The (T + T ind ) 2 term in the exponential ensures that the arch is initially at rest (see below). The parameter κ governs the rate of indentation. While taking κ → ∞ is analogous to the discontinuous indentation considered for the truss, we choose κ = 10 3 to avoid large truncation errors in the numerical scheme at very early times; see for example Morton and Mayers (2005) . In this formulation, the jump conditions at the indenter (due to the δ-function in (F.1) ) are automatically satisfied, and the unknown indentation force F does not explicitly enter the discretised equations. For both indentation and release stages, we avoid solving a DAE system by differentiating the end-shortening constraint twice in time, from which an explicit expression for the compressive force τ 2 can be determined (more generally see Ruhoff et al., 1996 ); the solution then satisfies the end-shortening constraint provided the initial data are compatible, which is guaranteed by starting fully relaxed and at rest in the natural shape when the indenter is first applied.