Happy Catastrophe: Recent Progress in Analysis and Exploitation of Elastic Instability

A synthesis of recent progress is presented on a topic that lies at the heart of both structural engineering and nonlinear science. The emphasis is on thin elastic structures that lose stability subcritically --- without a nearby stable post-buckled state --- a canonical example being a uniformly axially-loaded cylindrical shell. Such structures are hard to design and certify because imperfections or shocks trigger buckling at loads well below the threshold of linear stability. A resurgence of interest in structural instability phenomena suggests practical stability assessments require stochastic approaches and imperfection maps. This article surveys a different philosophy; the buckling process and ultimate post-buckled state are well described by the perfect problem. The significance of the Maxwell load is emphasised, where energy of the unbuckled and fully-developed buckle patterns are equal, as is the energetic preference of localised states, stable and unstable versions of which connect in a snaking load-deflection path. The state of the art is presented on analytical, numerical and experimental methods. Pseudo-arclength continuation (path-following) of a finite-element approximation computes families of complex localised states. Numerical implementation of a mountain-pass energy method then predicts the energy barrier through which the buckling process occurs. Recent developments also indicate how such procedures can be replicated experimentally; unstable states being accessed by careful control of constraints, and stability margins assessed by shock sensitivity experiments. Finally, the fact that subcritical instabilities can be robust, not being undone by reversal of the loading path, opens up potential for technological exploitation. Several examples at different length scales are discussed; a cable-stayed prestressed column, two examples of adaptive structures inspired by morphing aeroelastic surfaces, and a model for a functional auxetic material.


INTRODUCTION
Bernard Budiansky famously used to say "everybody loves a buckling problem" [1,2] and this resonates today in a number of significant ways. First, it reflects that buckling is a process of instability, and instabilities have always held a macabre fascination for children and adults alike. From the collapse of a pile of bricks to the deliberate demolition of buildings or the catastrophic failures of large urban areas in the aftermath of a natural disaster, instabilities hold a central place in human experience and consciousness. But there is a newer, more modern, contextthat structural instabilities can also be harnessed for the greater good. See, for example, the recent paper by Reis [3] who coined the phrase buckliphilia for the exploitation of instability in counterpoint to buckliphobia, the safety-conscious avoidance of collapse.
The purpose here is to review modern developments in the theory and analysis of buckling instabilities, both in the work of the present authors and by others. In the process, we draw particular attention to new techniques of analysis-often applied to classical thorny buckliphobic problems-and highlight potential areas of buckliphilic exploitation. We place particular emphasis on the interplay between analytical, numerical and experimental techniques, showing how we pick our way through a plethora of unstable post-buckling equilibrium states, to focus on practically relevant solutions.
With its origins in singularity, or catastrophe, theory [4,5], the non-linear analysis of structural buckling can be cast in the framework of static bifurcation theory. Broadly speaking, the instability of a trivial, unbuckled state, upon varying an external load parameter, falls into one of three categories: supercritical, subcritical or transcritical; see Figure 1. The former is sometimes called a safe bifurcation because, as shown in panel (A), the post-buckled path emerges smoothly out of the unbuckled equilibrium path, and hence stability can be safely tracked under slow variation of the applied load. In contrast, subcritical bifurcations, as seen in panel (B), have been termed dangerous [6], because the structure would irreversibly jump to a post-buckled state (not shown) that is a long way from the trivial one. Such jumps in elastic structures tend to give rise to energy loss, accompanied by a significant "bang, " and often lead to permanent non-elastic deformation or collapse. Transcritical instabilities (C) lie somewhere in between and are often associated with a loss of reflection symmetry. There is an extensive literature on understanding and classifying buckling instabilities, see for example Thompson and Hunt [7] and references therein. This paper shall concern instabilities that lead to large irreversible jumps.
Subcritical bifurcations, exemplified by the classic responses of thin elastic shells, are known to carry distinctive features, such as the likelihood of extreme imperfection-sensitivity and wide experimental scatter, and certainly merit a general overview. The canonical example is that of the axially-loaded cylindrical shell, of interest to rocket designers, aircraft and storage tank manufacturers, as well as in the construction of buffers to absorb mechanical energy. Here, instability under realistic conditions occurs significantly below the critical load of the system as determined from linear stability analysis of the perfect problem, absent from imperfections.
One approach to deal with such imperfection-sensitivity is through stochastic methods. Eliashakoff, Arbocz, and others pioneered developments, such as the international databank of imperfections (see [8,9]) and references therein. While such methods appeal at one level, from a modeling point of view there is also significant sensitivity to the precision of the chosen numerical method, and useful analyses typically require many Monte Carlo realizations. Also, from a practical perspective, estimating a safety margin of a particular specimen would necessitate comprehensive imaging and analysis of all its imperfections. Unfortunately for modern, lightweight composite structures, such imperfections often occur beneath the surface layer and are hard to characterize in practice. There has therefore long been a search for a lower-bound criterion below which a violently subcritical structure, such as a cylindrical shell cannot buckle. One such phenomenological idea is that of a "reduced stiffness" approach [10], but such simple design formulae still require a full understanding of the non-linear elastic equilibrium states of the shell.
Methods based on sensitivity to perturbations have received a recent resurgence of interest inspired by theories of critical transitions in fluid dynamics, see for example [11,12]. The focus in these works is sensitivity not to imperfection but to dynamic shocks; how large a dynamic perturbation in the form of a small tap or localized impact would be needed to induce buckling. It is not clear yet whether such ideas represent a realistic prospect for a practical non-destructive test for a particular specimen. Nevertheless, we shall argue in section 5 below that small localized shocks can allow engineers to explore the critical mountain pass unstable equilibrium that provides the route to buckling.
Fundamentally, this paper takes a deterministic rather than stochastic point of view. Starting with the perturbation methods introduced by the Dutch engineer [13], there is a rich tradition of using pseudo arc-length continuation to track unstable postbuckling solutions emerging from classical bifurcation points [14][15][16][17][18]. Supplemented by energy landscape considerations providing information on stability, these methods can be interpreted fundamentally as implicit analytical approaches to studying equilibrium solutions and their buckling paths. On the other hand, stochastic formulations, primarily explicit, after the introduction of quasi-realistic imperfection types and shapes and often employing Monte Carlo techniques, can be used to track expected dynamical behavior for particular specimens (see e.g., [19,20]). Clearly, both approaches carry inherent advantages and disadvantages, and should be regarded as complementary.
A key idea is that the perfect problem, devoid of imperfections or shocks, can give theoretical and practical insight into how structures buckle subcritically. We shall emphasize the significance of the Maxwell load, the level at the fundamental and periodic buckle patterns have the same energy (see e.g., [21]) and section 4 below. Nearby spatially localized equilibria are energetically preferable. But to find such states, we need to overcome an energy barrier, in the form of snaking or concertina pattern of unstable states connected by sequences of folds. Numerically, these patterns can be captured using pseudoarclength continuation, going back to Riks [22]. But how can we embed such a methodology in modern finite-element analyses? How can one access such unstable paths in an experiment? What is the best approach to understanding energy barriers? These are the questions this paper seeks to answer. We shall also keep in mind the perspective [3] that buckling instabilities, rather than to be avoided at all costs, can, in principle, be beneficial happy catastrophes.
The rest of this paper is outlined as follows. Section 2 gives a brief overview of non-linear post-buckling analysis of subcritical problems, starting from the pivotal work of Koiter, and including some general comments on analytical perturbation methods. A motivating simple pin-jointed "knee" model is presented as well as the classical problem of the axially loaded cylindrical shell. Section 3 surveys recent progress in computational path-following methods applied directly to a finite element representation to compute stable and unstable paths, with illustrations for a simple snap-through structure as well as the more complex cylindrical shell. Section 4 then considers computational energy-based methods that are able to identify Maxwell loads and mountain-pass solutions, again with reference to the cylindrical shell. Section 5 considers emerging experimental ideas to implement the numerical methods from the previous sections, via carefully controlled laboratory procedures. Section 6 surveys three examples, at different lengthscales and from distinct engineering domains, that attempt to exploit subcritical buckling instabilities: prestressed stayed columns, adaptive aeroelastic structures, and a structural model for auxetic materials. Finally, section 7 draws conclusions and suggests avenues for future work.

LARGE AMPLITUDE POST-BUCKLING ANALYSIS
Before the advent of modern computer-based methods, nonlinear post-buckling of elastic structures was largely dealt with by systematic asymptotic analysis; i.e., perturbation procedures based on Taylor expansions about the critical bifurcation point [7,23]. The pioneer of this field was Warner T. Koiter (1914Koiter ( -1997, beginning with his PhD thesis, completed during the Second World War [13]. Local perturbation analysis can be highly instructive in highlighting fundamental properties like underlying symmetries and symmetry-breaking, yet its range of validity typically is limited. This observation is by no means new. Koiter, for example, in discussing the buckling of a spherical shell under external pressure in 1969 [24], states: "In the problem of the spherical shell under external pressure the systematic perturbation procedure is only valid in a range of load factors within a fraction of the order h/R of the critical load, and for deflections of the order (h/R) 1/2 times the shell thickness. It follows that the systematic perturbation procedure at the critical point has little, if any, practical significance for the present problem." Here h represents the shell thickness and R the radius. He went on to suggest the following: "A far more powerful method of achieving a second approximation to the post-buckling behavior was also developed already in our earlier work [13, section 38]. It consists of an evaluation of the quartic terms in the energy expression not at the critical load itself, but at the actual value of the load factor under consideration." The reference is of course to his thesis [13], which did not appear in English until 1967 and was never published in the open literature. Koiter was a deeply humble man.
In more recent years Koiter's ideas have been supplemented by other asymptotic techniques, such as expansions at the so-called Maxwell load [25] and pseudo arc-length continuation which will be explained in section 3 below. Throughout this article we will apply our ideas to the recurring infinite-degree-of-freedom example of the axially-compressed cylindrical shell. But first, let us introduce a simple single degree-of-freedom example.

A Simple Motivating Example
Consider the simple mechanical model of Figure 2. A linear spring k is placed in-line with a "knee" comprising two finite-length rigid links hinged with a rotational spring, and compressed by an axial force P as shown. The rotational spring can take various characteristics-elastic, rigid-plastic or elastoplastic, as shown in the insets. It is also assumed that the arms can pass through one another without restriction.
The system has two degrees-of-freedom, with associated generalized coordinates and θ that respectively describe in-line displacement and rotation of the rigid link elements. It has three distinct possibilities for equilibrium. First, we have the simplest state in which the rotational spring does no work and the in-line spring simply squashes to give a fundamental equilibrium path describing the pre-buckling state: Second, under the condition that the knee rotation θ continues to grow in the positive sense, the potential energy function for the post-buckling response is either when the rotational spring is elastic with rotational stiffness K, or if it has passed its limiting plastic moment M p . In each case, the first two terms are the strain energies in the rotational and inline springs, respectively; the final subtracted term is the work done by the dead load P. When the rotational spring is in the plastic state, its energy contribution can be seen as quasi-strain energy, i.e., the work done in moving the joint through positive rotation, without necessarily being able to release this work if the rotation is reversed. Responses are then readily obtained in closed form from the two equilibrium equations ∂V/∂θ = 0 and ∂V/∂ = 0. The fundamental and post-buckling paths are plotted in Figure 3. for the three possibilities of Figure 2. There are several points worthy of note. For the elastic system of Figure 2A the initially-stable pure-squash fundamental equilibrium state reaches a supercritical bifurcation point B where it becomes unstable, whereupon one branch of the stable post-buckling path is then followed. Post-buckling analysis of this type of behavior responds well to the perturbation method [7].
The rigid-plastic system on the other hand has no bifurcation point, and this highlights one of the key issues to be addressed in this paper. Over much of the range, the fundamental and postbuckling paths, although being relatively well-separated in the P−θ plot, lie close to one another in the P− diagram. Moreover, although the paths approach each other asymptotically, they never meet; the bifurcation exists only at infinite load. Yet when the load is high, buckling could certainly be triggered by small imperfections or fluctuations. So the real problem becomes to determine a practical range of loads over which the system can be assumed to be relatively safe from instability. Another point to note in this case is that although the buckled path represents pure plastic behavior with zero stiffness, there remains an effective stiffness on this path beyond its minimum load point, once the arms have passed through vertical. This stiffness is due to geometric effects (deflection perpendicular to the load remaining Frontiers in Applied Mathematics and Statistics | www.frontiersin.org constrained). A number of circumstances are known where such geometrically non-linear, infinite-buckling-load problems arise in practice, for instance the buckling of railway tracks in heatwave conditions [26], and kink-banding of layered materials [27]. A major theme here is to review techniques for obtaining ball-park estimates for effective buckling loads and post-buckling information for such problems.
The bifurcation point is restored for the elasto-plastic system as shown in the bottom row of panels in Figure 3, where the behavior follows the purely elastic response of the top panels Frontiers in Applied Mathematics and Statistics | www.frontiersin.org until the rotational moment reaches its plastic limit, whereupon it switches to the rigid-plastic response of the middle row. The absence of a first bifurcation point is avoided, but an effective secondary bifurcation point S also needs to be negotiated. As with many interactive buckling problems there is the danger that the continuation of the elastic path remains as a possible equilibrium solution. This simple example emphasizes that modern pathfollowing computational methods need to be constantly vigilant in checking for further instability points and bifurcations as they track along an equilibrium path.

The Axially-Compressed Cylindrical Shell
Much has been written on the classical problem of the buckling of an axially loaded cylindrical shell (see e.g., [2,28]) and references therein. It carries all the hallmarks of a classical subcritical buckling problem, in particular its notorious sensitivity to imperfections. Its apparent simplicity, yet underlying complexity, has ignited significant academic dueling over the course of the last century (cf. different "resolutions" of the "paradox" by Zhu et al. [29] and Elishakoff [30]), and caused design engineers many a headache.
It might seem strange in a forward-looking review paper to focus on such a problem, but with the modern impetus toward ever stronger and lighter structures and new materials, understanding the cylinder response remains a fundamental issue of continuing research. We will therefore use axiallycompressed cylindrical shell buckling as the exemplar problem on which to illustrate the methodology reviewed in this paper.
When the fundamental deformation mode of a long, slender structure loses stability, it can either transition into a periodic buckling mode, spread equally over the domain, or a localized mode that is concentrated only over a portion of the domain. Such different kinds of buckle patterns were illustrated in the beautiful experimental and computational work of Yamaki [31]. It is well-known that localized post-buckling modes have a proclivity to develop in systems governed by subcritical bifurcations [32]. Moreover, Hunt and Lucena [33] showed that axially localized post-buckling modes exist for the axially compressed cylinder. This circumferential ring of an axially localized, diamond-shaped waveform can then undergo homoclinic snaking in the compressive loading parameter, leading to sequential ring formation along the length of the cylinder [34].
The work by Horák et al. [35] showed that an unstable equilibrium-localized axially and circumferentially-in the form of a single dimple is also possible (see Figure 4A). This state corresponds to the mountain-pass solution separating the stable pre-buckling and restabilized post-buckled states. The initial dimple, which can be found using the mountain pass algorithm of section 4, can undergo circumferential snaking, culminating in one ring of diamond-shaped buckles (see Figure 4B). This type of circumferentially-driven snaking is different from the axially-driven snaking described by Hunt et al. [34], where circumferentially-complete rings of buckles grow axially in a sequential manner.
We shall continually return to the squashed cylinder problem throughout this study.

Post-buckling Analysis-A Post-Koiter Reflection
Another canonical shell buckling problem that has received a resurgence of interest due to the recent work of Hutchinson [36] is that of the spherical shell under uniform external pressure. That problem too exhibits violently subcritical bifurcation and shock-sensitivity. In an extended paper published in 1969 [24] writes "An important result of Beaty's analysis [37] was that the numerical factor of the quartic term is much larger than the coefficient of the cubic term, indicating that the quartic term becomes already important for very small deflections in terms of the shell thickness, and that it is dominant over the cubic term for larger deflections. A similar evaluation of the quartic term in the energy expression at the critical load factor and for rotationally symmetric deformations was made independently by Walker [38], who also evaluated the next higher-order term, namely the quintic term, with an even larger numerical coefficient." He thus puts the poor performance of the perturbation method down to ever increasing influence of higher and higher-order terms-quartics larger than cubics, quintics more than quartics, and so on. This significant observation seems odd from the viewpoint of perturbation theory; using von Kármán-Donnell equations for the cylindrical shell [39] for example, a discrete formulation comprising doubly-periodic shapes generates energy terms only up to quartic level [7].
The need for higher-order terms can be explained through the process of elimination of passive coordinates, as espoused in the book by Thompson and Hunt [7]. Consider a conservative elastic structural system whose stable equilibrium configurations are described by minima of the energy function W(q i , ), where {q i } describes a set of n incremental generalized coordinates measured from a monotonically-increasing (fundamental) equilibrium path in loading parameter . Suppose that n − m of the q i are deemed passive (not actively involved in the buckling process). These passive terms are represented parametrically in terms of the remaining m active coordinates and loading parameter thus, q α = q α (q j , ), where now 1 ≤ j ≤ m and m + 1 ≤ α ≤ n. A new energy function W(q j , ), equal in value to the W-function, but written in terms of just the active coordinates and loading parameter, is introduced: Differentiation using the chain rule then gives derivatives of W in terms of W. Specifically, if W is diagonalized such that W ij = 0 for i = j, subscripts denoting partial differentiation with respect to the appropriate generalized coordinate, then derivatives up to cubic level pass over unchanged, but at quartic level we see contamination from lower-order derivatives. In particular, for a significant quartic term (see [7] for more details). Similar contamination from lower derivatives likewise appears at quintic level and above, leading to a lack of convergence as described in the Koiter quote above. The derivative W αα appearing in the denominator of (3) is the so-called stability coefficient for the passive coordinate q α , and would have equated to zero had the coordinate been active and directly involved in the buckling process, If critical loads tend to bunch together on the fundamental path, as occurs for both the axially-compressed cylinder and pressurized sphere discussed above, then contamination from higher-modes close to the critical point of interest can clearly be extreme.
Modal analysis in the form of spectral or pseudo-spectral numerical methods made a resurgence in the 1990s and 2000s, allowing numerical continuation (path-following) methods to scale to models with hundreds of degrees-of-freedom (see e.g., [21,40]). Nowadays, with modern computers being able to cope easily with millions of degrees-of-freedom, because of its geometric versatility the finite element method is the preferred technique for solving complex problems in engineering mechanics. Results provided in the next section are therefore presented in a finite-element setting.

NUMERICAL PATH-FOLLOWING FOR SUBCRITICAL INSTABILITIES
In applied mathematics, methods for multi-parameter analysis, branch-switching and bifurcation tracking are well-established theoretically using the language of catastrophe (singularity) theory [41] and differentiable dynamical systems (see [42]), including in infinite dimensions [43]. Using the concept of pseudo-arclength continuation due to Keller [44] these methods have been implemented numerically and incorporated into a variety of numerical continuation software packages, such as AUTO [45]. Typically such formulations apply to systems governed by ordinary or partial differential equations. In structural mechanics, specialized arc-length techniques were developed for non-linear formulations by Riks [14] and Crisfield [15]. Classically, those studies tended to be restricted to a single parameter-the applied load. However, the formulation can easily be extended to the general setting, as described here.
Our formulation considers a discretized model of a quasistatically evolving, conservative and elastic structure, where the internal forces, f (u), and tangential stiffness, K T (u), are uniquely defined from the current displacements, u, by means of the first and second variations of the total potential energy. Here u may represent all the degrees-of-freedom of a simple system, or a large-scale reduction of an infinite-dimensional problem. We first present a general framework, then illustrate the results for a simple toggle frame, before discussing implications for the cylindrical shell problem.

The General Setting
Equilibrium is defined as a balance between internal and external forces acting on the structure. In a displacement-based finite element setting, this balance is written in terms of n discrete displacement degrees-of-freedom u, and a scalar loading parameter λ: The vectors p(λ) and f (u) are the external (non-follower) load and internal force, respectively. In the case of linear and proportional loading we have p(λ) ≡ λp ,λ 1 = λp, wherep is a constant reference loading vector (dead loading).
The system (4) of n equations in (n + 1) unknowns-n displacement degrees-of-freedom and one loading parameter-is solved for a solution point, To turn this into a well-posed system of equations, one needs to add an additional scalar constraint, the most natural of which is the arclength constraint where n u and n λ take different forms depending on the nature of the arclength constraint. By linearizing about the current equilibrium state, x, and applying Newton's method for the iterative correction, δx, we can find a set of solution points describing a continuous equilibrium curve. Note that the partial derivative of the residual with respect to the displacement vector, F ,u = f ,u (u), is equal to the tangential stiffness matrix K T (u). More generally, Equation (4) can adapted to incorporate any number of additional parameters: where is a vector containing p control variables. Typically, 1 corresponds to parameters that influence the internal forces (e.g., material properties, geometric dimensions, temperature and moisture fields) and 2 relates to externally applied mechanical loads (e.g., forces, moments, tractions). The number n of equilibrium equations in Equation (7), correspond directly to the n displacement degrees-of-freedom. Because the structural response is parameterized by p additional parameters, a p-dimensional solution manifold in R (n+p) will be computed-the so-called equilibrium hypersurface [46]. By defining additional auxiliary equations, g, specific solution subsets on this p-dimensional manifold are recovered. In the general setting we therefore find solutions to the augmented system When r auxiliary equations are defined, the solution to Equation (8) is (p − r)-dimensional. Hence, r = p − 1 auxiliary equations are required to define a one-dimensional equilibrium curve in R (n+p) . Posing the problem in this general manner allows the structural response to be viewed not only as a function of a varying load but also as a function of other parameters that define the structure. By treating these additional parameters as "forcing" variables in an arc-length solver, their effect on the structural response is readily obtained.
This general treatment naturally lends itself to the tracing of loci of singular points in parameter space. To constrain the system of n equilibrium equations to such a locus, we simultaneously enforce a criticality condition, for example, i.e., at least one eigenvector φ of the tangential stiffness matrix K T spans the nullspace. In the most general form, a vector of q auxiliary variables v may be added to the auxiliary equations g. Hence, Equation (9) describes n equilibrium equations and r auxiliary equations in (n + p + q) unknowns leading to a (p + q − r)dimensional solution. To determine a one-dimensional curve of singular points, we thus require r = p + q − 1 auxiliary equations to constrain the system. Following the above example, when the n-dimensional null vector at the critical state is introduced as the auxiliary variable v, a singular curve in p = 2 parameters is appropriately constrained by the associated r = n + 1 auxiliary equations where the scalar equation restricts the magnitude of the eigenvector. When evaluating one-dimensional curves (r = p + q − 1), one additional constraining equation is needed to uniquely solve the system of for a solution point y = (u, , v) on the curve described by G(y). Hence, where N is a scalar equation that plays the role of a multidimensional arc-length constraint along a specific direction of the subset curve. Note that the system of equations for classical load-displacement equilibrium paths can be recovered by setting p = 1 and q = r = 0. A solution to Equation (10) can be obtained through a consistent linearization coupled with Newton's method, FIGURE 5 | Schematic diagram of a toggle frame under transverse load P causing the frame to snap through with displacement w. The height H of the toggle frame is a free parameter that can be varied in a numerical continuation solver to investigate the system's snapping behavior.
where the superscript denotes the jth equilibrium iteration and the subscript the kth load increment. The iterative correction cycle is typically started by a predictive forward Euler step. The above framework is quite general and can be adapted to find many different kinds of curves on an equilibrium surface; see [17] or [18] for further details. The key is to define pertinent auxiliary equations that constrain the equilibrium equation to the locus of points required. Examples include: 1. Classic equilibrium paths in load-displacement space (a loading parameter is varied). 2. Parametric paths in parameter-displacement space (a geometric, constitutive or secondary loading parameter is varied). 3. Pinpointing singular points (bifurcation and limit points) on either of the two paths mentioned above. 4. Bifurcated branches emanating from a bifurcation point. 5. Singular paths that describe a locus of bifurcation and/or limit points in load-parameter-displacement space. 6. Branch-connecting paths that connect points on distinct equilibrium curves, e.g., a fundamental and a bifurcated path.

Illustrative Example-Snap-Through of a Toggle Frame
As an example, consider the snap-through behavior of the centrally loaded toggle frame with clamped ends, shown in Figure 5. We start this idealized model with pre-defined geometry, material properties and loading to illustrate how the general algorithms can be used for a comprehensive investigation of structural stability and design parameter sensitivity. We shall evaluate the frame's fundamental load-displacement behavior, including pinpointing all relevant singular points. Additional non-singular and singular curves are then traced by starting from a chosen solution on the fundamental path to explore the surrounding design space. The toggle frame initially deforms symmetrically on the fundamental equilibrium path. This deformation mode becomes unstable at a symmetry-breaking bifurcation just before the maximum limit point on the curve. Because the connected nonsymmetric path branching from the bifurcation point is unstable, the toggle frame snaps dynamically into the inverted stable shape. In Figure 6, blue segments denote stable equilibria, red segments denote unstable equilibria and black dots denote critical points. Figure 6A restricts path-following to the classical displacement-load space. To illustrate generalized pathfollowing capabilities, Figure 6B extends the analysis to changes in the height H of the frame. Figure 6B shows an isometric view in displacement-load-height space of the fundamental and bifurcation paths, plus two additional parametric paths. For these parametric paths, the applied load is held constant at P = 37.4 and P = 64.8, respectively, and the relationship is traced between the height H and central displacement w.
By imposing a singularity condition in the generalized pathfollowing algorithm, the locus of limit and bifurcation points can be traced, illustrating how changes in the height of the frame affect the load-displacement solution of these singular points. There are multiple benefits of tracing such fold lines. First, they can be used to identify interesting points, such as the coincidence of limit and bifurcation points-the hilltop-branching points at H = 0.567 and H = 0.581-or points where bifurcation and limit points cease to exist-the cusp catastrophes at H = 0.506 and H = 0.346. These points are clearly marked in the orthographic projections of Figure 6C (w vs. H) and Figure 6D (P vs. H). Second, fold lines can be used in design studies to determine the sensitivity of singular points with respect to design parameters, without having to perform computationally expensive Monte Carlo studies. Finally, fold lines can be used for optimization purposes. For example, the displacement at the first instability can be maximized by reducing the height of the toggle frame to coincide with the hilltop-branching point at H = 0.567 (see Figure 6C).

Application to the Cylindrical Shell
Consider a thin-walled isotropic cylindrical shell of thickness t = 0.247 mm, radius R = 100 mm and length L = 160.9 mm loaded in uniform axial compression via displacement control. The cylinder is linear elastic and isotropic with Young's modulus E = 5.56 GPa and Poisson's ratio ν = 0.3, chosen to model Yamaki's longest cylinder (Batdorf parameter Z = L 2 √ 1 − ν 2 /Rt = 1, 000) [31]. To represent a typical experimental setup as closely as possible, the cylinder is rigidly clamped at both ends with axial compression/displacement u imposed at one end of the cylinder and the other end completely constrained.
The cylinder is modeled using isoparametric, geometrically non-linear finite elements based on a total Lagrangian formulation. The finite elements used are so-called "degenerated shell elements" [47] based on first-order shear deformation theory assumptions [48]. The cylinder is discretized into 97 axial and 241 circumferential nodes, that are assembled into 25-noded spectral finite elements using the element formulation of Payette and Reddy [49] and the large rotation parametrization described by Bathe [50]. To reduce computational effort and complexity, only a quarter of the cylinder is modeled. The circumferential domain is described by s/R ∈ [−π, π] and the axial domain by x/L ∈ [−0.5, 0.5] such that we model the quarter-segment s/R ∈ [0, π] and x/L ∈ [0, 0.5] with reflective symmetry conditions applied along the lines of symmetry. In all figures that follow, blue segments denote stable equilibria, red segments unstable equilibria, and black dots critical points.  Figure 7A shows the equilibrium path starting from a single dimple, superimposed on the pre-buckling curve (fundamental path) in terms of normalized axial compression (uR/Lt) vs. normalized load (P/P cl ). The classical buckling load is given by P cl = 2πEt 2 / 3(1 − ν 2 ). The stable pre-buckling curve runs diagonally in blue with the unstable single dimple solution running almost coincidentally alongside it. The unstable equilibrium branch of the latter starts at a limit point close to the first critical point on the pre-buckling path. This limit point is denoted by 0 in Figure 7A with the corresponding normalized radial (out-of-plane) displacement (w/t) shown in Figure 7B. The deformation mode clearly shows a localized dimple in the center of the domain.
Path-following in the direction of decreasing displacement leads to a snaking sequence. The reason behind snaking has been established in a number of related contexts as the behavior of homoclinic orbits in the unfolding of a heteroclinic connection between flat and periodic states (see [21]) and references therein for an application to structural mechanics. The phenomenon is closely related to the notion of front pinning around a Maxwell point, first described by Pomeau [51] in the context of fluid dynamics. For a recent review of snaking see [52].
Starting from limit point 0 in Figure 7A, the single dimple becomes more pronounced with decreasing end-shortening and stabilizes at a limit point (uR/Lt = 0.293, P/P cl = 0.479). This critical point corresponds to the smallest possible compression to allow a single dimple as an equilibrium solution. Tracing the equilibrium path further, a series of destabilizations and restabilizations add further buckles to the left and right of the original single dimple. Proceeding along the snaking path, the single dimple thus grows in a sequence of 1, 3, 5, 7, and 9 waves until an entire ring around the cylinder exists. The mode shapes corresponding to limit points I-V in Figure 7A are shown in Figure 7B and depict the series of increasing odd-numbered buckles (1, 3, 5, 7, and 9) spreading around the cylinder circumference. This snaking sequence of odd buckles connects to another equilibrium path that preserves an additional symmetry group at pitchfork bifurcation point PB. This additional path is described next.
The equilibrium path in Figure 7C is an additional snaking sequence starting from two sets of two dimples located to the left and right of the original dimple (see Figure 7D). In Figure 7C the snaking path of the single dimple (from Figure 7A) is shown in gray for reference and the new equilibrium path starting with two sets of two dimples is shown in red/blue. The four-dimple snaking path also originates at a limit point (O) close to the first critical point of the pre-buckling path (see inset A of Figure 7C). With decreasing end-shortening, the snaking sequence grows from 4 to 8, and finally to 10 buckles. The mode shapes corresponding to various points on the red/blue path of Figure 7C are shown in Figure 7D. The two equilibrium paths (gray and red/blue) are seen to connect at a pitchfork bifurcation (point D in inset B of Figure 7C). In the immediate vicinity of this connection, the four-dimple snaking path regains stability at a limit point. The ten-buckle waveform is stable from uR/Lt = 0.268 until it destabilizes at a pitchfork bifurcation (point E) at uR/Lt = 0.626. Beyond this bifurcation, an additional snaking sequence occurs, leading to the full Yoshimura post-buckling pattern. Additional rings of buckles all initiate from a single localization and then spread circumferentially (see [53] for more details). An additional snaking sequence starting from two dimples and representing growth of an even number of waves (2, 4, 6, 8, and 10) also exists. The even snaking sequence mirrors the behavior of the odd snaking sequence in its pattern formation and in the connection to another equilibrium path at a pitchfork bifurcation. In systems featuring spatial localization, snaking of both even and odd number of localizations is typical [52] and these solutions are often intertwined. This behavior is also confirmed for the axially compressed cylinder and is shown in Figure 8. Figure 8A shows the equilibrium path of the even snaking sequence in red/blue superimposed on the snaking solution of odd buckles in gray (from Figure 7A). The even snaking sequence is also broken away from the pre-buckling equilibrium path and starts with the formation of two adjacent inward buckles (see point ∅ Figure 8B). These two buckles then multiply throughout the snaking sequence, with the equilibrium paths of the even and odd snaking sequences intertwined. The different mode shapes corresponding to limit points i-v in Figure 8A are plotted in Figure 8B and show the series of increasing even buckles (2, 4, 6, 8, and 10) growing around the cylinder circumference.
The snaking solution of even buckles also ends at a pitchfork bifurcation (point PB in Figure 8A) where it connects to another segment of the equilibrium path. This connecting equilibrium path is shown in Figure 8C with the even-buckle path from Figure 8A superimposed in gray. The two segments (red/blue and gray) of the equilibrium path connect at a pitchfork bifurcation (point d in inset B of Figure 8C). The mode shape of point d in Figure 8D confirms the expected ten-buckle waveform. In conclusion, both the odd-and even-buckle curves lead to an axially localized post-buckling state of a single ring of ten diamonds.
In closing this section we remark that the snaking results for the present paper were obtained with a mesh 4× denser than those in reference [53]. While the overall behavior of the snaking sequence and the nature of pattern formation is unchanged, the refined results presented here update and eliminate the secondorder snaking features originally observed in reference [53].

ENERGY BASED METHODS-MAXWELL AND MOUNTAIN-PASS CRITERIA
While continuation methods are an integral part of unraveling the often complex behavior associated with subcritical instabilities, they do not tell us which state (energetically) the system would prefer, nor quantify the sensitivity of a locally metastable equilibrium state, so-called shock sensitivity. We can therefore supplement our continuation approach with energy-based methods to explore these questions more directly.

Maxwell Load vs. Maxwell Displacement
To address the problem of an infinite pre-buckling critical load (as in Figure 3B), one option is the Maxwell equal-energy criterion. This criterion originates from the concept of the Gibbs free energy in thermodynamics, to address state transitions triggered by statistical fluctuations or disturbances [54, p. 53].
Here, stability is governed only by the global minimum of total potential energy. Under some form of parametric variation, such as load or applied displacement, the Maxwell criterion provides the first circumstance for which the energy of a post-buckled state first falls below the energy minimum of the pre-buckled state. The reasoning is that, upon increase of the parameter, this would be the first time the system could be shaken out of its trivial state and transition to the post-buckled global minimum. Although the Maxwell criterion cannot be considered a hard and formal point of instability, it may, nevertheless, serve as a useful and robust lower-bound estimate for instability in systems where small disturbances and imperfections have pronounced effects. A long structural system loaded axially typically prefers a localized to a distributed post-buckled response [32]. If the localized buckle subsequently restabilizes, additional cells of buckling will often develop in adjacent positions to the first via a sequence of localized instabilities, according to the snaking mechanism described for the cylinder above. Typically, the load tends to fluctuate or snake between upper and lower limits as the sequence progresses. The Maxwell load, defined (as above) as the lowest load for which the post-buckled energy matches its prebuckled counterpart, lies between the two limits, and effectively acts as an organizing center about which the post-buckled load oscillates as the snaking progresses.
This snaking sequence with localizations developing over the length of a structure has now been recognized in a number of different circumstances (see for example [55][56][57][58][59]). However, in section 3.3 we describe an alternative snaking scenario, in which localized buckles trigger not axially but orthogonal to the direction of the applied load, around the circumference of a buckling cylindrical shell. In this sequence, a Maxwell displacement rather than a Maxwell load acts as the organizing center, with the system fluctuating between two limits of endshortening as the load continues to fall. Two examples of such snaking behavior are seen in Figures 7A,C.

Mountain Pass Algorithms
First introduced by Ambrosetti and Rabinowitz [60], the Mountain Pass Lemma is a fundamental mathematical tool for proving the existence of stationary points of non-linear functionals. Excluding some technical details, the key ingredients of the theorem are: 1. a suitable (energy) functional W(x) 2. a stationary point e 1 , which is a local minimum 3. a second point e 2 , for which W(e 1 ) > W(e 2 ) We note that a suitable function is normally available in the form of total potential energy, with local minima appearing on a stable fundamental equilibrium path [7].
The theorem states that, over the set of all continuous paths connecting e 1 and e 2 , i.e., one can find the infimum of the maxima of the energy functional W(u) along any path γ ∈ Ŵ. This infimum is the mountain pass solution and is a saddle point The physical significance of the mountain pass is that it represents the connecting point in solution space with the smallest energy hump, required to escape the local minimum at e 1 and transition to a lower energy state at e 2 . Therefore, the Maxwell load/displacement (depending on the loading regime), at which W(e 1 ) = W(e 2 ), marks the onset of the ability to jump to such a lower energy state, and therefore the onset of "shock sensitivity" in the system.
The application of the Theorem provides a computable energy hump to assess shock sensitivity; the mountain pass state x c itself is significant, since at this point the system has just one negative eigenvalue for which the system is unstable. This eigenvector marks a direction in solution space m tangent to the mountain path γ , and suggests a mode shape that if applied to the system would most easily induce transition from e 1 to e 2 . This eigenvalue at the mountain pass point therefore indicates the imperfection or probing modes that a subcritical system could be most sensitive to.
The literature gives a variety of algorithms for finding mountain pass solutions e.g., the nudged elastic band method [61], the dimer method [62], and conjugate peak refinement [63]. Here we briefly describe the latter, as it is used later to illustrate shock sensitivity of the axially-compressed cylinder. Conjugate peak refinement is an iterative scheme performing alternating line search maximization and minimization steps to find the mountain pass solution x c . The approach generates a sequence of piecewise-linear approximations to a path γ ⋆ which passes through x c . For the k th iteration, we denote this approximation γ (k) characterized by a set of points Ŵ (k) : = . We start the process by defining γ (0) as the straight line connecting e 1 and e 2 , so that Ŵ (0) = {e 1 , e 2 }. Then starting with k = 1, each iteration comprises three steps: 1. Line Search Maximization: We maximize the functional W(x) along the piecewise linear path γ (k−1) , to obtain line maximization pointx (k) lying between points x (k−1) j and x (k−1) j+1 .

Line Search Minimization:
We then find the scalar α such that W(x (k) + αs) attains a minimum, where the search direction s is chosen to be: At any iteration of the statex (k) + αs is the best approximation of the mountain pass solution. Convergence can be tested by ensuring (1) the gradient of the energy at this point zero (2) determining the lowest two eigenvalues of the energy's Hessian, since only the smallest must be negative. We now demonstrate the mountain pass procedure geometrically with a generic, two degree-of-freedom energy landscape given by a modified Müller-Brown potential [64]. This particular potential (Figure 9) has application in computational chemistry [63]. It is chosen since it has no symmetry, and is characterized by two local minima at e 1 and e 2 , with a non-trivial mountain pass connecting them. The approach provides a good approximation to the saddle point in just two iterations. In the first iteration, we see the algorithm starts by approximating the mountain path with a straight line between the e 1 and e 2 . A maximum is located along this line; see the • in the far left panel of Figure 9. A minimum in the conjugate direction (12) is then found, as indicated by the symbol in the next left-most panel. Thus, the first iteration provides a reasonable approximation to the saddle. The path γ (0) is updated to γ (1) , characterized by three points connected by the pair of straight lines. For iteration 2 the procedure continues in the same way, first a maximization step over the path to produce the second • in the third panel of Figure 9, followed by a minimization problem in the conjugate direction to the path through the maximum point. For this example, the newly found minimum on this path (the second in the right-most panel of the figure) satisfies the tolerance condition ∇W(x c ) e < ǫ = 1 × 10 −5 , and the algorithm terminates.

Application to the Cylindrical Shell
Classically, stability of equilibrium is governed entirely by the local Hessian of the total potential energy; wells with respect to all degrees-of-freedom denote stability, whereas saddles or maxima denotes instability. This framework fails in the case of infinite critical loads presented in Figure 3, where two equilibrium paths are separated by a small energy barrier but never strictly intersect. Thus, although the idealized structure never transitions out of the trivial equilibrium state, that state becomes metastable, where small disturbances could trigger instability.
As is seen in the insert of Figure 7A, the cylinder too features an equilibrium path that asymptotically converges to, but never actually intersects, the pre-buckled path. Here, we compare the energy levels of the pre-buckling path and the circumferentially periodic equilibrium path of a single ring of diamonds (path ending in point E in Figure 7D). As load is applied in a rigid manner (controlled end-shortening), the stability threshold corresponds to a Maxwell displacement and this computes to be M u = u M /u cl = 0.486 (u M R/Lt = 0.294). It is interesting to note that this value correlates well with the limit point (uR/Lt = 0.294) in Figure 7A.
The Maxwell displacement could serve as a lower-bound estimate for the cylinder's first instability load, by marking the onset of "shock sensitivity" [65]-i.e., that of metastability of the pre-buckling path. This snaking sequence marking the development of a single dimple into a buckle pattern that is periodic circumferentially but localized axially, is in marked contrast to the snaking at the lower Maxwell load of such rings developing to a fully periodic pattern identified in earlier work [40,66].
Rather than apply a computationally expensive infinite degree-of-freedom mountain-pass algorithm as in Horák et al. [35], another way of looking for the mountain-pass solution is to test the cylinder's resilience to the singledimple localization; how much displacement is necessary to trigger a dynamic escape from stable pre-buckling to a post-buckled state via the mountain-pass saddle? We could envisage perturbing the cylinder from the side using a hypothetical infinitesimally-thin, infinitely-stiff, probe or poker; experimental implementation of this idea is explored in section 5 below.
To implement such an analysis numerically, consider applying such a poker at right angles to the cylinder midsurface, half-way along its length. Such a process involves two fundamental parameters, applied end-compression u and lateral probing force F. We consider applying such a probe repeatedly as the axial compression is quasi-statically increased. The results are presented in Figure 10. At low levels of axial compression, we find a non-linear softening/stiffening relationship of strictly positive stiffness between the probe force F and the ensuing dimple displacement w; see path 1 in Figure 10A. Here w denotes radial displacement relative to the radial (Poisson) dilation that naturally occurs in the pre-buckling state.
For increased levels of end-shortening, the equilibrium manifold traces S-shaped curves; as the dimple develops, lateral resistance reduces, until limit points are traversed leading to regions of negative stiffness (paths 2-3 in Figure 10A). For even greater end-shortening, the probe force reduces significantly, dipping below the zero load axis (e.g., F = 0 on path 4). At this point an unstable saddle state is encountered, corresponding to the single-dimple mountainpass solution. Also shown in Figure 10A is a black fold line connecting maximum and minimum limit points, and thereby describing a boundary that separates the domain into stable and unstable regions. Figure 10B expands this landscape into three dimensions, providing an interesting stability landscape that qualitatively matches the experimental results of Virot et al. [67] on a different cylinder. The area between the stable pre-buckled and unstable single-dimple solutions under the F vs. w curve represents the energy barrier, and thereby the "shock sensitivity" of the prebuckling state. The size of this energy barrier can be understood qualitatively by plotting the fold line connecting maximum and minimum points; see Figure 10B. This curve slopes down toward the buckling point on the pre-buckling path (point CL). Indeed, the fold line intersects the buckling point on the pre-buckling path, confirming that resilience of the pre-buckling state to small perturbations (i.e., the linear stability) indeed vanishes at that point.

EXPERIMENTAL METHODS FOR EXPLORING INSTABILITIES
While numerical methods for the analysis of non-linear structures are well-developed, experimental methods tailored to such structures, in particular shell buckling, have received comparatively little attention (see e.g., [20,31,[67][68][69][70][71][72]). The trend in modern engineering is to test experimentally for single parameter values, and then use computational models, virtual testing, and "digital twins" wherever possible to extend the envelope. However, for fundamentally imperfection-sensitive buckling problems, such an approach will not explore the stability landscape reliably. Hence there is a fundamental barrier to researchers hoping to exploit non-linear structures concepts for industrial applications. This barrier has arguably led to over-conservative designs in safety-critical industries like the aerospace sector, which requires stringent testing before new components are allowed to fly. Organizations, such as NASA have therefore put renewed emphasis on experimental testing, and in particular, feeding high-fidelity imperfection measurements of specimens into models [72].

Experimental Path-Following
Conventional test methods fail to capture all but the simplest non-linear behavior, and consequently researchers lack reliable methods to validate their ideas experimentally. The main reason traditional test methods fail is the difficulty in measuring unstable parts of the response. Any structure whose equilibrium curve features limit points can snap under force-or displacementcontrolled test methods, as illustrated in Figure 11A. Such snapthrough thus gives rise to regions of the equilibrium curve that might seem inaccessible experimentally.
Numerical analysis succeeds where experiments fail because in a numerical setting the force and displacement at a control  point can be controlled independently and simultaneously. This freedom allows the solver to set combined limits on force and displacement, and prevents jumping to other solutions when an equilibrium becomes unstable (see Figure 11B for an illustration of the arc-length method). Consequently, displacement limit points can be traversed, unstable paths followed, and the full non-linear response of a structure described. In an experimental setting, the force and displacement of a control point are linked by the elasticity of the structure; meaning that one can control the displacement of a control point and generate a reaction force, or vice versa, but not both. This coupling makes it impossible to apply numerical techniques to experiments without some additional control.
There are several interesting published approaches to work around this problem. Wiebe and Virgin [73] use a hammer to trigger snap-through of a shallow arch. By analyzing the transient dynamic trajectories of the structure during the snap, locations of unstable static equilibria are deduced. By intentionally allowing snap-through, this method can locate unstable equilibria without needing to actively control them. Virot et al. [67] use a poker to laterally probe a cylinder under increasing axial load. By tracking changes in the probe force-displacement curves, they can estimate the load at which the cylinder becomes globally unstable, before the instability load is reached. The concept of probing is especially relevant to the work presented here: an experimental path-following method which utilizes probes to stabilize and control unstable equilibria.

Application to a Shallow Arch
Consider the centrally-loaded shallow arch studied by Neville et al. [74] shown in Figure 12A, with dimensions L = 205 mm, h = 20 mm, t = 1.5 mm, and depth = 5 mm. For symmetric deformations, the structure exhibits the complex non-linear behavior shown in Figure 12C. Looking at the first few "petals" (starting at one of the two fundamental equilibria and following the equilibrium path toward the other), it is clear that the response comprises many successive displacement limit points. A displacement-controlled experiment would only obtain the first two segments of the equilibrium curve (the solid lines in Figure 12D), snapping from one to the other at limit points L 1 and L 2 .
At u a = 5 there are several equilibria available; each with distinct values of F a . Each equilibrium is also associated with a unique deformation shape (Figure 12E), where the unstable equilibria correspond to more complex shapes. Controlling the structural shape allows us to stabilize unstable equilibria, and indirectly influence the reaction force on a displacementcontrolled control point. Having decoupled the force and displacement experimentally, numerical approaches become viable. Two extra probes provide control over the deformation shape as shown in Figure 12B. Such probes produce reaction forces that change the structure. However, when the probe reaction force equals zero, then it is an equilibrium of the unperturbed structure. Crucially, reaction force is also zero at unstable equilibria-in this case the probes provide the infinitesimal restoring force required to prevent the instability from growing. By applying a large probe perturbation at a stable equilibrium, the probes can be used to search for other equilibria.
By moving the probes and actuation point in concert, a simple form of path-following can be performed [75]. The actuation point steps forward, then the probes search for equilibrium (F p = 0). Small perturbations can be used to avoid large deviations from the equilibrium curve. If the actuation point steps past a limit point, the probes will not be able to find FIGURE 14 | Test for assessing the resilience of the cylinder to perturbations, proposed by Thompson [65] and Thompson and Sieber [11]. equilibrium and the actuation point direction is then reversed. This approach allows the equilibrium path of the shallow arch to be followed around a displacement limit point, as shown in Figure 13B. The FEA prediction is also shown in Figure 13A, for comparison. Deformation shapes of the arch at several points in the experiment are also shown in Figure 13C. Shapes 1 and 2 are the two fundamental equilibria (also shown in Figure 12C). Shape 3 corresponds to a stable equilibrium, and resembles shape 1 in Figure 12E. Shape 4 corresponds to an unstable equilibrium, and consequently is more complex. It resembles shape 2 in Figure 12E, which corresponds to the same segment of the equilibrium curve.
Experimental results are naturally affected by phenomena and imperfections not included in theoretical models. The shallow arch example, for instance, is sensitive to changes in geometry and probe location, as well as displaying complex behavior in response to the two input parameters (u a and u p ). Virtual testing is a technique that can address these issues, and aid in experimental design and interpretation of results. A successful example of such a virtual testing environment coupled to the commercial FE solver ABAQUS is presented in Groh et al. [75]. A finite-element model is used to simulate a "perfect" experiment-i.e., one in which the equipment and test specimen behave exactly as intended. The model includes limitations of the experimental setup-e.g., sensor noise, limited number of control points, etc.-and provides the same inputs and outputs that are available in the real experiment. This virtual testing environment provides a useful middle ground between numerical solutions and experiments, and serves as a "sandbox" or digital twin to explore the effects of different test configurations and imperfections.

Experimental Mountain-Pass Methods-"A Game of Poker"
Inspired by the theoretical work of Horák et al. [35], Thompson et al. [11,65] observed that the mountain-pass state for an axially-compressed cylinder looks remarkably similar to the small dimple induced by probing the cylinder laterally with a finger. The authors discuss a thought experiment where the cylinder's resilience to perturbations (i.e., linear stability) is tested by probing the cylinder radially inwards using a poker (see Figure 14).
In fact, the idea of poking axially-compressed cylinders from the side to assess resilience to buckling has a long history predating any mountain-pass considerations. By tapping axiallyloaded cylinders with a finger, Eßlinger and Geier [77] found stabilized single-dimple states. Hühne et al. [71] realized that the single dimple can act as an imperfection, that excites the characteristic observed buckling behavior. They set out to out to determine a robust lower bound to buckling load that could replace the empirically determined knockdown factors proposed in NASA's SP-8007 [78] guideline. Over the last decade, this methodology has led to a battery of tests on composite cylinders [79], and extensions to probabilistic perturbation [80] and multiple perturbation loads [81].
The poker force vs. displacement response of the cylinder for different levels of axial compression was shown previously using FE simulations in Figure 10A. Once these S-shaped curves pass through the zero-poker-force axis, the single dimple exists as an unstable saddle solution and a dynamic escape via the mountain pass is possible. The likelihood of such an escape can be quantified by the mountain-pass energy barrier, represented by the area enclosed by the associated force/displacement curve and force axis. In Figure 10A it is readily observed that the greater the applied axial compression, the smaller the energy barrier provided by the single-dimple mountain pass. Thus, by repeating the poking procedure for multiple levels of axial compression and computing the energy barrier up to the mountain-pass point, a non-destructive testing method can be established that provides a safety cushion before buckling is likely to occur.
Such a probing experiment was successfully implemented by Virot et al. [67]. By controlling the displacement of the lateral poker, the unstable mountain-pass point was determined when the reaction force on the poker vanished; this state replicating that of an unprobed cylinder. Furthermore, by aggregating the poker force vs. poker displacement curves for multiple levels of compression, a stability landscape emerged that qualitatively matched (Figure 10B). An experimental buckling load of the cylinder could also be determined by recording the level of compression for which it lost its ability to resist the poking displacement, i.e., the reaction force fell beneath a specific tolerance.
Even though the idea is simple to implement, in practice the system can bifurcate by pivoting around the point load. To offset  this symmetry-breaking effect, as highlighted by Thompson and Sieber [11], a second control probe is often required. Second, the choice of probing location needs to be carefully chosen, as different locations can lead to differing buckling load predictions [82]. Finally, there are practical shortcomings to obtaining the mountain-pass state for arbitrary systems. For the cylinder, the mountain-pass state happens to be a relatively simple single-dimple localization, but in general, it could be more intricate, and may be difficult to impose by one or even multiple pokers. Combining poking experiments with experimental pathfollowing algorithms may well prove a fruitful avenue for future research.

EXPLOITING BUCKLING INSTABILITIES
As stated in the Introduction and indeed reflected in the title of this paper, instability need not solely be considered as something to be avoided or designed against; it is also possible to utilize instability in a positive manner [3]. We highlight three areas where such benefits can be found in different engineering domains and length-scales.

Prestressed Stayed Columns
Prestressed stayed columns are important elements of many modern large-scale structures; see Figure 15.
Such columns tend to be slender and have intermediately placed cross-arms and associated pretensioned cables, thereby reducing the buckling effective length L e . The length L e provides a measure of the critical buckling eigenmode wavelength and the Euler strut buckling load is proportional to 1/L 2 e . The system of cable stays and cross-arms in the prestressed stayed columns provide intermediate restraints that reduce L e significantly and hence provide a commensurate increase in critical buckling load and ultimate capacity. Depending on the overall geometry, this change in critical buckling load can also be associated with quantitative and qualitative changes in the triggered buckling mode within the non-linear range. The behavior has been discussed at length in previous work, with the focus falling on qualitative critical [83] and post-buckling behavior determined by the pretensioning force [84], physical experiments [85], triggering of modal interactions and associated symmetrybreaking [59,86] and tuning behavior for different cases [87,88].
Some of these works use conventional finite element modeling, where post-buckling shapes are initiated by introducing imperfections that are affine to linear buckling modes. A drawback is that the full picture of modal interaction only becomes available under a combination of symmetric and anti-symmetric imperfections. Other modes may also be drawn in, for example, should it be thin enough, localized buckling in the main tubular column, and the numerical methodology discussed in section 3.3, can be useful. Nevertheless, there has been a sequence of increasingly-sophisticated low-dimensional models, to capture mode interaction [59,89]. These models are based on the Rayleigh-Ritz method, with a finite number of linear eigenmodes of the unstayed column being used to discretize the deformation response, with cross-arms acting as beams, and cable stays as tension-only members.
The particular complication in stayed columns is produced by the cable stay, where there is the possibility of a sudden loss in elastic stiffness caused by cables slackening. The outcome is similar to that described in Figures 2D, 3C. Using numerical continuation it has been possible to track post-buckling paths of perfect systems within the package AUTO-07P [45]. Figure 16 shows a realistically-proportioned symmetric stayed column system of length L with three cross-arms of lengths a e and a m , the exact details of which are presented in Yu and Wadee [59]. It shows that the mechanical response is likely to destabilize when distinct modes are triggered or mode interaction dominates; the kinks in the post-buckling paths are signatures of the portions of the cable stays going slack, instantly losing axial stiffness, and causing a sudden unstable jump in the response.
The column subsequently restabilizes once it finds a configuration that restores equilibrium. Both the numerical continuation procedure for the analytical model, and the Riks algorithm used in ABAQUS, can capture this behavior. One advantage of the former is that it tends to crystallize the detailed mechanical response into a few distinctive characteristics; the main column buckling modes are discretized into a Rayleigh-Ritz type model, and the non-linear results provide straightforward output of the contributions of the linear buckling modes to the post-buckling profile, Q 1 and Q 2 being amplitudes of the first two main column buckling modes. This analysis allows the interpretation of the effects of symmetry-breaking, and the potential to trigger higher pure or interactive modes in the post-buckling range.
All the consequences for the post-critical strength, stiffness and potential to jump between different equilibrium states owing to the cable stay behavior, can be determined directly. This information can then be used to determine parametric spaces where practical geometric quantities, such as stay diameter, layout of the stayed column system and initial prestressing forces, can generate qualitatively different, yet predictable, responses [59,87], as shown in Figure 16.
The simplest configuration with a single-cross arm can also be considered as a single cell within a larger lattice material. The performance of metal lattices, for example with a criss-cross structure as in Queheillalt and Wadley [90], may be enhanced by including internally woven tension ties, to make them similar to woven composites [91]. The behavior of the sandwich panel of Figure 17, depends on non-linear interactions between the individual cells, and is observed to have a similar structure to the stayed column. By adjusting the pretension in the ties at the production stage, alongside the geometry of the cells and the overall configuration, it is possible to engineer the post-buckling stiffness to a desired level. If, for instance, the post-buckling stiffness of an overall panel can be tuned to be practically zero [92], but with each buckling cell having a significant critical load, then the structure would be highly effective in absorbing energy. Moreover, since stress propagation depends on structural stiffness, if the post-buckling stiffness on the panel were zero, it would potentially minimize stress transfer to any attached structure. This type of application, using internal buckling of a lattice material for dynamic isolation purposes, has potential for applications where lightweight reusable elements are required for impact and blast protection or seismic isolation [93].

Adaptive Structures
So-called adaptive structures are able to change shape and/or material properties in response to varying external stimuli [94]. The application of adaptivity in engineering has the potential for significant improvements in performance, by making structures more efficient over a broad range of operating conditions [95]. A fascinating natural example of shape-adaptivity is the Venus flytrap, whose rapid transition from open to closed to capture its prey occurs as a consequence of snap-buckling instability [96].
Hence, when operating conditions require large displacements and/or multiple stable configurations, structural instabilities can be viewed as potentially advantageous rather than a source of failure.
Consider for simplicity the buckling response of a simplysupported Euler strut, illustrated in Figure 18. In particular, consider the stable post-buckled configuration in the inset of Figure 19A. For given combinations of compression and symmetry-breaking defects, such a post-buckled structure can behave like an arch, exhibiting dynamic "snap-through" behavior between the two stable states when subjected to an external transverse load, F. For a beam compressed just beyond the first buckling load, these configurations connect via an S-shaped equilibrium path in the space of a centrally-applied force, F plotted against the midpoint respective deflection δ, as shown graphically in Figure 19A. The structure initially deflects in a stable manner before reaching a maximum limit point, where it snaps dynamically through a region of instability (2) onto the equilibrium branch for its second stable configuration (3).
Compression levels required to produce any meaningful shape-change from δ to −δ are typically sufficient to cause snapthrough at a symmetry-breaking bifurcation (see asymmetric red shape in Figure 19A). Figure 19B shows how the equilibrium path in the F vs. δ plane connects the two stable branches of the broken pitchfork. Figure 19B also shows how the limit points of the F vs. δ equilibrium path change as a function of the compressive displacement, u, via the black fold line.
Specifically, the fold line tracks the two limit points with respect to changes in the compressive displacement, u, thereby illustrating the border between stable and unstable equilibria. By reducing u, the two limit points of the equilibrium path in Figure 19A collide in a cusp singularity (see Figures 19C,D). This cusp singularity therefore determines the critical value of compression, u, at the onset of dynamic snap-through behavior [97]. Indeed, depending on the value of compression, u, three distinct types of post-buckling behavior can be observed when the transverse load, F, is applied: 1. For values of compression, u, greater than the limit point on the broken-away pitchfork branch, the structure snaps from its first stable shape to its second configuration, traversing the region of instability delimited by the fold line. A selfequilibrated second configuration exists (stable even when F is removed). The structure is said to be bistable. 2. Reducing the compression, u, into the region between the limit point on the broken-away pitchfork branch and the cusp singularity, allows the beam to traverse a region of instability when F is applied, thereby still exhibiting snapthrough behavior. However, at these values of compression, it does not have a second stable configuration for F = 0. When the external force is removed, the strut then snaps back to its primary state. In other words it shows "superelastic" monostability. 3. By decreasing the level of compression, u, even further the structure deforms non-linearly, displaying stiffness adaptation but without snap-through. It is then elastically monostable or simply stable.
The control of geometrical parameters, material properties and/or boundary conditions can be used to tailor the equilibrium manifolds and adapt the multistability of the system to specific working and environmental conditions. We now consider a practical example. A passive adaptive air inlet can regulate the opening aperture of a connected duct by interacting with the fluid flow around it. As shown in Figure 20, the inlet comprises a deformable, glass-fiber panel poised in an open state. The panel has been buckled into the region between the limit point of the brokenaway branch and the cusp (refer to the taxonomy above), making it monostable.
As the airflow streaming over the panel accelerates into the connected duct, the decreasing pressure field creates an upwards force on the panel causing it to snap shut at a critical airspeed.  If the airspeed is lowered beneath another threshold, the inlet automatically opens again. Unlike traditional shape-changing systems, the inlet does not rely on auxiliary devices for actuation. By increasing the amount of compression beyond the limit point on the broken-away path, the inlet can be transformed into a bistable structure that remains closed once the airflow is reduced. The greater the applied compression, the higher the airspeed required to actuate snap-through; the system's parameters can be tailored to meet specific operating requirements (see [98,99] for further details).
This device has potential for engineering applications where cooling and drag reduction create competing design drivers. Examples include air inlets on cars or cooling ducts on jet engines, which use fresh-air cooling for reliable engine operation although this cooling induces a drag penalty. For additional engineering examples that use the non-linear taxonomy described above (see e.g., [100][101][102][103][104][105][106][107]) and references therein.

Buckling-Induced Auxetic Materials
The fact that materials exist with negative Poisson's ratio is not only intriguing but also of practical significance. So-called auxetic materials typically exhibit high energy absorption and fracture resistance, and have a broad range of practical applications from blast curtains and shock absorbers to running trainers and the ability to control waves (see for example [108][109][110][111][112][113] and references therein). Auxetic materials are known to occur naturally (various honeycomb structures for instance or even crumpled paper [113] but with recent technological advances in 3D printing they can also be readily manufactured. This ease of manufacture opens exciting new opportunities to tailor such materials, and work has exploded in this area. To show something of the flavor of this exciting new field, we now briefly outline two examples of auxetic behavior studied recently by the present authors, induced by local instabilities within the material's structure. Experimental and numerical work led by Bertoldi et al. [111] has demonstrated how mechanical instabilities in periodic porous structures can lead to the dramatic reorganization of the material from the original configuration, giving rise to the auxetic possibilities. One widely studied structure involves the uni-axial loading of a square array of circular holes made of an elastomeric matrix, which can be readily manufactured using a 3D printer. The hexagonally-shaped sample shown in Figure 21 at the right is manufactured with an OBJET Connex 3D printer (Stratasys Ltd., USA), a machine that employs PolyJet Matrix Technology to dispense material-in this case Tango, a rubber-like elastomerfrom designated micro-scale inkjet printing nozzles [114,115]. When a compressive load is aligned with the square array, geometric reorganization is seen, as the elastic instability induces periodic deformation patterns of tessellating ∞-shaped voids. At the macroscale, this generates a non-linear auxetic response by the simulations shown in blue in Figure 21. Interestingly if our hexagonal specimen is rotated through 45 • , so that the compressive load acts on a diamond arrangement of circular holes, the material response is close to that of classical linear elastic (non-auxetic) material (red line).
Most examples of auxetic behavior in the literature are based on re-entrant structures [108,109,113,116], and we next briefly review a recent contribution [117] detailing a variant that allows phase transitions from auxetic to non-auxetic phases and vice versa, based on the unit cell shown in Figure 22. The cell comprises two back-to-back single degree-of-freedom arches with displacements measured by Q 1 and Q 2 , linked by the linear spring of stiffness k [118]. Four different aspects of the full nonlinear response, including the important plot of load against its corresponding displacement of panel B, are given in Figure 23.
The continuous smooth curve replicates the response of a single arch [118], apart from the fact that under displacement control it goes unstable at bifurcation point A. Here the stable post-buckling solution becomes non-homogeneous, with displacement in one of the arches outpacing the other. This asymmetry continues until the first component passes through the horizontal and begins to stiffen, whereupon displacement in the other starts to take over. Symmetry is again restored at bifurcation B.
The absence of homogeneity in the natural loading path gives the potential for considerable complexity of response once cells are combined, as in Figure 24. The responses are shown in Figure 25. With a relatively modest extension of the single cell model, a complex tangle of stable and unstable equilibrium paths is generated. Further details can be found in Hunt and Dodwell [117].

CONCLUSION
This theme issue focuses on the notion of stability in a variety of different contexts, both mathematical and practical. It could be argued that there is no more classical context in which one thinks of stability than structural engineering. It is fundamentally the job of the structural engineer to avoid buckling, failure or collapse. This paper takes a slightly different point of view on the topic. We focus on emerging ideas of elastic stability and post-instability behavior of structures that fail subcritically, via irreversible jumps in energy. Such problems are of current interest for at least three reasons.
First and foremost, because of their sensitivity to shocks and imperfections, there is difficulty in certifying such structures for safety. We have argued that despite over 70 years since Koiter's pioneering work, a robust methodology for analyzing the stability of such structures has yet to emerge. We have promoted here a promising line of attack, based on the Maxwell equal energy criterion and the concept of the mountain pass, as well as emerging ideas on how such ideas might be applied experimentally. However, there remains much to be done before such ideas can provide a practical assessment and design tool. It is also interesting to note how the method relies on understanding the structure of unstable, localized post-buckling paths, which form the energy barriers or basin boundaries of the problem. In that sense, there is a strong connection to other active areas of stability-related research; tipping points in natural systems (see e.g., [119]), and transition to turbulence in pipe flow and related fluid-dynamic problems (see e.g., [120]).
Second, the structural engineering domain is changing. Across numerous lengthscales, there is a quest for ever more lightweight structures. It could be said that the revolution in composites and other nano-structured materials has been threatening to revolutionize just about the whole of the built environment for almost 50 years. Yet, despite the huge investment within academia and industry, why are we not yet seeing carbon fiber motor cars come off the production line, wholly composite airplanes in our airports, or fiber-reinforced polymer buildings being constructed en masse? There are doubtless a range of reasons for this slow penetration of composite technologies, and as most disruptive technologies, the revolution may actually be just around the corner. Nevertheless, we would argue that one of the bottlenecks still to be overcome, is that we do not understand how such structures fail. Most lightweight structures are optimized for strength, but such optimization typically leads to subcritical failure modes (see e.g., [28]). But, for structures made from composite materials, to the uncertainty and sensitivity of classical steel and concrete structures that buckle subcritically, we have additional complexities of anisotropy, internal micro-structural lengthscales and buried failures, cracks or delaminations that are hard to characterize and inspect. We have also seen, through the example of an auxetic material, that the distinction between materials and structures is fundamentally being blurred. It would seem then, the seeds of the robust methodology we have been trying to sow in this article are more important now than ever, if the true potential of nano-structured materials is to be realized in the built environment. We hope that future researchers may be inspired to water these seeds.
Finally, there is the point of view that we have been also trying to promote in this article that instability is not necessarily a bad thing. We have highlighted three areas of possible engineering exploitation of non-reversible structural instability. More generally though, we are quite used to the notion of things that snap and pop into instability. These include the pressure required to depress the keys on a computer keyboard being controlled by dome buckling, to old-fashioned bi-metallic strips being used to control switches, as in a motor car indicator light. Crash barriers and crumple zones also exploit the idea that elastic deformation of a subcritical structure can lead to transfer of significant amounts of energy into permanent plastic deformation. Origami also provides an inspiration to engineers in how small energy barriers need to be overcome in order to fold (or unfold) a structure into a new shape (see e.g., [121][122][123]). Most interestingly, there is great potential to draw inspiration from biology. Irreversible transitions are the norm in biology, for example in cell division, cell polarity formation and most morphogenesis problems. Although such processes are often controlled by genes and other signaling proteins, there is an increasing body of work that looks at the biomechanics of such transitions. Indeed, at many different lengthscales, processes that are crucial to development, or to the maintenance of an organism or ecosystem require sudden, irreversible response to continuous variation of external conditions. It is easy to think of examples like the springing of a Venus fly trap, or the opening of seed pods that can easily be regarded as buckling events of the nature described in this paper.
Clearly, there remain many lessons that engineers and designers need to learn by taking inspiration from the natural world. Not least among such lessons, as we seek to build a more resilient world in the face of global change, must surely be that there need not necessarily be anything to fear from an instability. Not only are sudden irreversible instabilities not necessarily to be feared, they can in fact be designed to be exploited for the greater good. Happy catastrophes indeed!

DATA AVAILABILITY
All data generated for this study can be found in the manuscript and/or supplementary material.