Harnessing energy landscape exploration to control the buckling of cylindrical shells

The complexity and unpredictability of postbuckling responses in even simple thin shells have raised great challenges to emerging technologies exploiting buckling transitions. Here we comprehensively survey the buckling landscapes to show the full complexity of the stable buckling states and the transition mechanisms between each of them. This is achieved by combining a simple and versatile triangulated lattice model for modelling the shell morphologies with efficient high-dimensional free-energy minimisation and transition path finding algorithms. We show how the simple free energy landscapes of short, lightly compressed cylinders become vastly more complex at high compressive strains or aspect ratios. We then exploit these landscapes to introduce an effective method for targeted design - landscape biasing. This is used to inform thickness modifications enabling landscape redesign, and the development of structures which are highly resistant to lateral perturbations. Our methods are general, and can be extended to studying postbuckling responses of other geometries.

Here, we demonstrate how the buckled states and buckling transitions of complex systems can be comprehensively surveyed and controlled. Our key methodological contribution is to combine a simple and versatile triangulated lattice model for modelling the shell morphologies with efficient high-dimensional free-energy minimisation and transition path finding algorithms, in order to develop a powerful computational methodology for exploring the buckling landscapes. We apply this approach to the problem of cylindrical shell buckling, where the extreme landscape complexity arises from a combination of subcritcality, multiplicity and snaking in the postbuckled states [17]. Harnessing these tools we also explore the landscape biasing technique as an effective method to design and control buckling responses.
We begin by surveying the (meta) stable states -the free energy minima in the landscape. It has long been recognised that the cylindrical postbuckling states are strongly subcritical; coexisting with the unbuckled states in a loading interval which spans between the lower buckling load [23], and critical load (see [24] for a detailed review). This means that, as has been revealed historically, cylindrical shells are capable of failing at even 20% of their critical load [25]. Furthermore, the particular sensitivity to lateral loads [26], lead to NASA's development of empirical predictions for the practical load bearing capacity of imperfect cylinders [27]. Previously, buckled states have been elucidated by solving the von Kármán-Donell equations relating the stress to the radial displacement in an elastic cylindrical shell, but only by assuming the solutions exhibit axial periodicity (see for example [28,29]), reminiscent of the diamond pattern shown by Yoshimura to enable global, inextensible buckling [30]. Similarly, group theory has also enabled the study of high-symmetry solutions [31]. However, a plethora of postbuckling solutions exist, discussed recently in the context of spatial localisation of the elastic deformation leading to snaking (pinning) in the solution space [17].
The energy landscape approach is then used to survey how the minima transform into each other by the lowest energy routes. We connect all postbuckling morphologies via such transition pathways, and so explore the entire stability landscape of cylindrical buckling. In this we reveal a diverse variation in the landscape properties, ranging from  Figure 1: Summary of the stable buckling morphology classes and the minimum energy states, surveyed across a range of aspect ratios A 0 and end shortening ratios λ. a-c Visualisations of representative minima, shown in 3D and as radial displacement field contour plots. The axial and angular coordinates z and θ of the contour plots are shown in c(i), with the displacement d expressed as a fraction of R 0 . d Phase diagram indicating the global free energy minimum across a range of A 0 and λ. The control ratio, k stretch R 2 0 /k bend is fixed throughout at 2.5 × 10 5 . The global minimum is either unbuckled (grey region, simulation data shown as black circles), or multiply dimpled (pink region, simulation data shown as red squares). The singly dimpled state is never the global minimum, but the existence region is shown outlined in blue with unfilled data points. very simple funnel-shaped landscapes at low aspect ratios and end shortenings, to broad and highly complex glassy landscapes at long aspect ratios.
Previously, only the first transition capable of buckling the unbuckled cylinder has been investigated [32,33]. This has received significant interest as capturing the minimum energy pathway (MEP) enables the minimum energy barrier to be obtained, which provides an absolute lower bound to the energy required for a compressed cylinder to buckle. An explicit link has therefore been made between the ease of single dimple formation and the sensitivity of loaded cylinders to lateral loading [32]. As important for structural applications, it has been suggested that these theoretical minimum energy barriers can be accessed experimentally via a local probing technique for cylindrical [26,[34][35][36] and spherical shells [37,38]. However, it is only via comparison to the minimum energy pathways obtained here that we are able to verify this.
Finally, we introduce a new method to begin to exert control over the buckling landscape -landscape biasing. In this, we are able to stabilise or destabilise targeted features in the landscape, such as transition states and minima. This is achieved by making local modifications to the elastic spring constants in the triangular lattice model to simulate thickness modifications. Thus, the knowledge of the energy landscape proves highly complimentary to experimental processes aimed at exerting postbuckling control [39][40][41]. We demonstrate the principal of landscape biasing by first showing how biasing against the unbuckled-single dimple transition state produces a 20% increase in buckling resistance of the unbuckled cylinder for a 1% increase in mass. We then show how biasing for a multiply dimpled state simplifies the local landscape, tripling the targeted state stability at 0% mass change.

Free energy minima
The triangular lattice model, detailed in Methods, discretises the shell into a triangulated mesh of extensional and angular elastic springs. Respectively, these allow for the decomposition of the total free energy into a sum of stretching and bending terms. To begin with, the stretching and bending spring constants, k stretch and k bend , are uniform throughout the shells. Each shell is generated with a welldefined aspect ratio A 0 = L 0 /(2R 0 ), where L 0 and R 0 are the length and radius of the cylinder when all springs assume their equilibrium configurations. When axially compressed, the shortening ratio is defined λ = L/L 0 , where L is the length of the compressed cylinder. The top and bottom edges of the cylinder are simply supported: the coordinates of the mesh are fixed, but the planes attached to the ends can bend freely. We also choose k stretch , k bend , and R 0 to maintain a constant dimensionless elastic con-trol ratio k stretch R 2 0 /k bend =2.5 × 10 5 , and free energies E reported throughout are nondimensionalized such that the reduced free energy E r = E/k bend . The elastic control ratio is chosen to be representative of a physical system, corresponding to an aluminium drinks can. The choice of elastic control ratio and nondimensionalization are discussed further in Methods.
The construction of the free energy landscape begins by surveying the free energy minima. To access the many different buckled states for each fixed aspect ratio A 0 and shortening ratio λ, a basin hopping step is employed prior to energy minimisation [42], detailed in Methods. Three characteristic cylinder morphologies are observed: unbuckled, singly dimpled, and multiply dimpled, visualised in 3D and as radial displacement fields d in Fig. 1a-c respectively.
The multiply dimpled states in Fig. 1c form the largest set of minima, within which is contained the often-studied morphologies of high rotational symmetry -the Yoshimuralike diamond dimpling pattern, an example of which is shown in Fig. 1c(ii) [30]. However the largest multiply dimpled subset is the irregularly dimpled morphologies, a characteristic example of which is shown in Fig. 1c(iii). A random perturbation applied to these cylindrical shells is therefore most likely to result in an irregularly dimpled state, showing that cylindrical shell buckling responses are inherently hard to predict.
The phase diagram in Fig. 1d summarises the minimum survey. At high λ in the phase diagram, the global free energy configuration is the unbuckled state, indicated by black circles. Upon decreasing λ the multiply dimpled states become the global minima, indicated by red squares. The solid red line indicates the point at which the buckled and multiply dimpled states are isoenergetic, at which the shortening ratio therefore produces an axial load equal to the Maxwell load (a detailed discussion regarding the loading limits is given in [24] for example). Across all tested scenarios, the most stable multiply dimpled states are those exhibiting a high degree of rotational symmetry, most commonly those with Yoshimura-like diamond patterns. However, we also find examples where more exotic high-symmetry multiply dimpled states form the global free energy minima, such as the example shown in Fig. 1c(i) at A 0 = 10, λ = 0.999.
The singly dimpled state, shown in Fig. 1b is of significant interest due to it's frequent role in the first buckling transition (which we consider further in the proceeding section), and also it's characteristic role of being the unit excitation in the postbuckling landscape. However, across a broad range of aspect ratios and elastic control ratios, detailed further in Supplementary Note 1, we observed that the single dimple is never the global free energy minimum. When it is energetically unfavourable to form a dimple, the unbuckled state is lower in energy; when it is energetically favourable to form a dimple, the energy is always lowered further by subsequent dimpling. The single dimple is therefore only metastable. This metastability region is outlined in blue in Fig. 1d. The non-monotonic form of the low-λ boundary arises from the complex deformation profile surrounding the dimple. At high aspect ratios, this profile extends around the circumference of the cylinder, such that self-interaction effects contribute to the dimple stability (detailed further in Supplementary Note 1).  (1). c Three example transition pathways connecting the unbuckled state and 1row by 9-dimples (1×9) state illustrated in d. All pathways are shown to begin with the unbuckled -single dimple transition, which is magnified in a, and the number of dimples are labelled at each minimum in Path A.

Buckling transitions
In order to describe the minimum energy mechanisms by which the cylindrical buckling morphologies interconvert, we must obtain the minimum energy pathways (MEP). Between any two states in the free energy landscape, the MEP is defined as a path in which the gradient of the free energy is parallel to the path tangent vector. The MEP will also pass through at least one saddle point in the landscape, a local energy maximum along the pathway. The buckling morphology at this point is known as the transition state. Several methods exist for finding the MEP and transition states, see for example refs. [43][44][45][46][47]. The string methods we use here are detailed in Methods.
Computationally, the only transition which has been followed previously is the simplest unbuckled-singly dimpled pathway, where the dimple is centrally located on the cylinder [32,33]. Meanwhile, local probing of cylindrical shells has been suggested as an experimental technique which may allow the true dimpling transition state to be accessed [26,[34][35][36]. In Fig. 2a, we compare the reduced-energy profiles E r (s) of the MEP (black series) with the pathway generated by simulating the local probe technique (blue line) for an example cylinder with A 0 = 0.8, λ = 0.9986. Local probe simulation methodologies are detailed in Methods. In order to usefully compare the paths, the path distance coordinate s is the Euclidean distance between the triangulated mesh of a point along the pathway, with that of the initial (unbuckled) state where a i and a o i are the position vectors of node i in the buckled and unbuckled mesh respectively.
On comparison, we observe that the local probe technique does meet the MEP at the transition state (labelled '*' and shown in shown in Fig. 2b), but does not access the minimum energy pathway generally. At the point of crossing the barrier, the locally-probed system snaps to a dimpled-like configuration: a small probe displacement resulting in a large change to the surrounding morphology, and a concomitant jump in E r and s.
This comparison shows that the local probe technique is capable of measuring the minimum energy barrier to the first dimpling transition. This is consistently shown across all the test cases summarised in Supplementary Note 2. Previous studies were unable to prove that the local probe technique could access the minimum energy barrier, as it was assumed that the true MEP was not too curved [26]: namely the direction of motion along the transition always has a component in the direction of the applied force (i.e the path never curves against the applied force).
However, the methodology presented here allows for the pathway between any two states to be investigated, not only the 0-1 transition. We therefore extend the first pathway found in Fig. 2a to find complete pathways from the unbuckled state to the multiply-dimpled global minimum. Examples are shown in Fig. 2c, with Fig. 2d showing the pathway endpoints: the unbuckled state, and the (1×9) global minimum. Two key observations are made: multiple competing pathways exist between the end points, and each pathway is complex, featuring many intervening minima. Out of the large number of possible pathways, three examples are highlighted in Fig. 2c, labelled A, B, and C. Movies showing the conformational changes along each pathway are shown in Supplementary Movies 1, 2, and 3 respectively. Path A is distinguished from other paths: out of the set of barriers along path A, the maximum energy barrier is the smallest out of all possible pathways. In Path A, eight separate dimpling transitions occur. In the first seven, a single dimpling event occurs to build a train of dimples. The final transition sees two dimples forming simultaneously to complete the ring of nine dimples. In this final transition, the path distance decreases as all dimples become shallower on formation of the final two. However, the system is capable of undergoing dimpling transitions not linked to the growing dimple train, leading to example alternative pathways B and C.

Energy landscapes
By connecting any pair of minima with an MEP, we may thus explore the complete energy landscape for any fixed A 0 and λ. Here we examine the extent of the landscape complexity as a function of A 0 and λ (varying the elastic constants is presented in Supplementary Note 3). As will be shown, cylindrical shells exhibit a diverse range of landscape types. We will first compare the energy landscape of a lightly compressed short cylinder where the single dimple is stable (A 0 = 0.8,λ = 0.9986), with a heavily compressed short cylinder where the single dimple is unstable (A 0 = 0.8,λ = 0.9980). We then compare the short, lightly compressed cylinder, with a long, lightly compressed cylinder (A 0 = 3.0,λ = 0.9990), where the single dimple is stable in both cases.
As the network of minima connected by MEPs is in general highly complex, it is instructive to consider simplified network representations. In Fig. 3, the free energy landscapes are visualised as disconnectivity graphs (for a comprehensive discussion of the disconnectivity graph representation of energy landscapes, we refer the reader to refs. [42,48]). In this, the network of minima and pathways is reduced to a spanning tree showing only the energy of the minima (the end points of each branch) and the lowest energy barrier connecting any two minima, read by tracing the path between two branches and finding the highest energy point. For example, in Fig. 3a, unbuckled state and singly dimpled state are labelled '0' and '1' respectively. On tracing between the two branches, the highest energy point along the path, labelled '*' marks the largest transition state energy. In this case, this is the 0-1 transition state shown in Fig. 2b. However, as the 1D disconnectivity graph does not show which states are directly connected, in general the highest energy point between two states is simply the largest energy encountered in the possible multi-step transition pathway.
In Fig. 3a, the disconnectivity graph is presented for A 0 = 0.8, λ = 0.9986, and represents the full energy landscape which was partially described in Fig. 2. Under these subcritical conditions, the unbuckled, singly dimpled, and multiply dimpled states coexist. However, the buckling landscape is remarkably simple: qualitatively, the states are (approximately) uniformly distributed across the stable energy range. To quantify this and subsequent observations, we partition the minimum-energy range into 100 bins of equal width and total the number of minima within each bin; this histogram is shown in Supplementary Note 4. We then calculate the variance in bin populations as a measure of the distribution uniformity. Here, the small variance in the bin frequency, 0.38, describes a relatively uniform distribution of minima across the energy range. The uniformity of the landscape is further reflected in the range of energy barriers  The example minimum (i) is a characteristic state of the system, featuring clusters of dimples closely aligned around the central circumference. The (1×9) global energy minimum (highlighted in red) exists in a deep well, with the minimum energy barrier greater than the first transition by a factor of 7. Thus, if an unbuckled state is subject to perturbations with sufficient energy to overcome the first dimpling transition, although other states may be sampled along the way, the tendency is to quickly become trapped in the global energy minimum. The notable exception to this picture however is that a second deep branch also exists at the base of the disconnectivity graph. This represents a competing set of deep states which are likely to split the population between the lowest minimum (1×9), and second-lowest minimum (1×8), labelled (iii) in Fig. 3a.
Upon decreasing λ to 0.9980, although the system is still subcritical, the singly dimpled state looses stability. The disconnectivity graph for this landscape is shown in Fig.  3b. Here, the landscape is markedly different to the lesscompressed case shown in Fig. 3a: although the number of minima is O(10) in both cases, at λ = 0.9980 the majority of states are concentrated at the lower stable energy range, indicated by the greater variance in bin population, 2.09, detailed further in Supplementary Note 4. Additionally, the range of energy barriers is large, varying from 10 −3 to 10 1 , with many states featuring high energetic barriers. This latter point is most pronounced when considering the (2×11) multiply dimpled state, labelled (i), which has an energy barrier 1000× greater than the minimum energy barrier from the unbuckled state. A further contrast in this disconnectivity graph is that the global minimum (2×9) does not have a large energy barrier compared to other transitions. Thus, random perturbations made to the unbuckled state may result in the system becoming trapped in several states different from the global minimum. Two highlighted examples of these which are close in energy to the global minima are the (2×10) system, labelled (iii), and a defective (2×9) system with two adjacent dimple vacancies, labelled (iv).
Finally, we return to a subcritical shortening ratio where the unbuckled, singly dimpled, and multiply dimpled states coexist, but now extend the aspect ratio: A 0 = 3.0, λ = 0.999. The disconnectivity graph for this system is shown in Fig. 3c. Three prominent features of this landscape offer significant contrast to the short aspect-ratio landscapes: the number of minima has increased by a factor of 100 compared to the A 0 = 0.8 systems, the minimum distribution is highly non-uniform -the bin population variance is 93, and the landscape becomes rough over a range of energy scales.
Expanding on these observations, the increase in the number of minima is due to two effects. Firstly, at large aspect ratios, all minima observed are no longer characterised uniquely by a single well-defined energy and morphology, but exist as clusters in which the intra-cluster energy variability is approximately ∆E r < ×10 −3 . Thus, on the finest scale, the stability landscape is rough and glass-like. In the stability landscape shown in Fig. 3c, we have clustered minima which share the same number of dimples with interconversion barriers < 10 −3 , reducing the number of minima shown by a factor of 10. The second effect is due to dimple confinement introduced by the fixed ends. At A 0 = 0.8, the fixed ends tightly constrain the dimples to lie within either one or two rows, due to the characteristic dimple size being similar to L 0 . At the longer aspect ratio of A 0 = 3.0, the constraining strength of the fixed ends is diminished, yielding a larger number of possibilities of dimple arrangements.
The large phase space for dimple arrangements within certain energy ranges enables numerous minima to exhibit similar energies and similar barriers. This is most pronounced in the range 10.1 < E r < 10.3, dominated by irregular systems with between 7 and 11 dimples. A representative example is shown, labelled (ii). In this region, the number of dimples is large enough to produce a significant number of variations in arrangement, yet not so large that packing constraints become dominant. On average, the inter-cluster energy barrier is O(10 −2 ). A similar glassy region exists at larger energies, where irregularly dimpled systems feature between 3 and 6 dimples. A representative example here is shown, labelled (i). Thus, the stability landscape becomes rough on two energy scales: (1) ∆E r ≈ ×10 −3 associated with intra-cluster variability, and (2) ∆E r ≈ ×10 −2 associated with inter-cluster variability in the absence of packing constraints (when comparing clusters of similar numbers of dimples). The distributions of energy barriers associated with this roughness are shown in the Supplementary Note 4.
For larger dimple numbers than 11, efficient packing on the cylinder is required, leading to a severe reduction in the phase space of dimple arrangements. Thus, in the vicinity of the global minimum, the (2×6) regularly dimpled state highlighted in red, the local landscape becomes significantly less glassy. Nonetheless, the overall landscape roughness coupled with a large number of deep states means that a perturbed unbuckled cylinder may buckle to any number of states, explaining the difficulty in designing cylindrical postbuckling states.

Controlling the landscape
Despite the complexity of the buckling landscapes, we now demonstrate how to control the stability of target features, by introducing a process we term landscape biasing. This enables us to design buckling responses by locally thickening or thinning the cylinder, complimentary to experimental realisation; see for example [39][40][41]. We demonstrate two examples of landscape biasing, by first biasing against a target transition state, and then biasing for a target minimum. The examples shown here significantly increase the stability of the target structures to lateral perturbations. These biased structures are therefore highly suited to scenarios where sudden morphological changes would be detrimental to device performance, a key example being aeronautical applications [27]. For these examples, we apply this method to A 0 =0.8, λ=0.9986 system, for which the buckling landscape is shown in Fig. 3a.
To begin with, it is observed that the minimum energy barrier from the unbuckled state to the singly dimpled state is small compared to both the overall landscape energy range, and other deep states, generating the extreme imperfection sensitivity of cylinders to sub-critical buckling transitions. The energy profile for this transition, shown originally in Fig. 2a, is re-plotted in Fig. 4a (solid black line), in which the reduced energy is referenced to the energy of the unbuckled state, E o . We aim to increase the energy barrier of this transition, in order to make the unbuckled cylinder more robust against lateral perturbations, by biasing the landscape against the transition state.
The landscape biasing workflow is shown in Fig 4b-d, and detailed further in Supplementary Note 5. Firstly, as shown in Fig. 4b, we obtain the radial deformation field for the unbiased transition state (as well as that of the unbuckled state). Secondly, we compute the fractional change in local elastic potential energy E f when transforming from the unbuckled to the transition state. It is observed that the stored elastic potential energy is highly localised about the centre of the dimple deformation. We then reason that in order to increase the energy of this transition state (and hence the barrier to the transition), we must modify the cylinder to energetically penalise this localisation of the potential energy, effectively biasing the landscape against the transition Biasing amplitude / % 0 5 10 15 20 ( 1 x 8 ) (1x9) state. A choice exists in how to perform this modification, but for this example we choose to simulate a local thickening of the shell by modifying k stretch (∝ t) and k bend (∝ t 3 ), facilitating experimental realisation. A more sophisticated yet complex treatment would alter k stretch and k bend independently, according to the separate local stretching and bending energies respectively. A comparison of alternative geometric methods to modify cylindrical shell buckling are presented in [39]. In the local thickening treatment, detailed in Supplementary Note 5, we weight the thickening according to the local energy change. Due to the symmetry breaking of the transition, in order to suppress dimple formation anywhere around the circumference of the cylin-der, at each z we average the thickening profile over all θ. Finally, the thickening profile is rescaled in order to achieve a prescribed total mass increase, which is set as 1% for the results presented in Fig. 4. The final thickening profile is shown in Fig. 4d, which sees the a thickness increase localised around the centre of the cylinder.
On attempting to dimple this biased cylinder, the transition state is now forced off-centre, shown in Fig. 4e. The energy profile for this transition is shown as the solid red line in Fig. 4, showing that for a 1% increase in mass, a 20% increase in buckling resistance is achieved. This improvement is over twice that of a uniformly thickened cylinder, 9%, with the same mass increase, the transition profile for which is shown as the dotted black line. This landscape biasing against the transition is the antithesis to modal nudging [49], the recently formalised technique for slender structures in which minimal structural modifications are made in order to select a specific failure mode.
The second way to design the bucking landscape is to bias for a target structure. We observe the landscape shown in Fig. 3a to exhibit a deep global minimum (1×9) and the shallower (1×8) state. Here, we choose to stabilise the (1×8) state through minimum-targeted landscape biasing. It will be shown how a target minimum can be significantly stabilised, thus realising a postbuckled state which is highly resistant to lateral perturbation. Furthermore, this example will show that through biasing we can select which highsymmetry morphology forms the global minimum.
In Fig. 5a, we show the radial displacement field of the (1×8) state. As before, we evaluate the local stored elastic potential energy, then weight the local elastic constants to exact a local thickening, detailed further in Supplementary Note 5. As the (1×8) state is to be stabilised, in regions of high stored elastic energy we locally thin the structure to reduce the energetic cost of the specific buckling mode. We also weight the thickening so that there is no overall mass change, and prescribe a biasing amplitude -the maximum percentage change in thickness allowed. To obtain the local thickness change, we therefore scale the weighting field w shown in fig. 5a by the biasing amplitude.
By systematically increasing the biasing amplitude from 0% to 20%, we observe how the buckling landscape changes at the bottom of the funnel, shown in Fig. 5b. At 0% bias, we show a magnification of the low-energy portion of the disconnectivity graph shown in Fig. 3a, featuring the two deep wells decorated with multiple stable minima. The wells corresponding to the (1×8) state and (1×9) state are shown highlighted in red and blue respectively. In Fig. 5c, the percentage change in the (1×8) and (1×9) barriers are shown relative to their respective barriers at 0% bias.
On application of a 5% bias, the landscape changes significantly relative to the unbiased case: the landscape is simplified as the biasing destabilises many minima, the (1×9) state increases in energy, and the targeted (1×8) state decreases in energy to such an extent that it becomes the global minimum. Furthermore, the landscape simplification and (1×8) state stabilisation effects act cooperatively to increase the barrier out of the target (1×8) state by 207% relative to the unbiased (0%) landscape. At 10% bias, these effects are further magnified. At 20% bias, there is no further change in the lower landscape structure, but the stabilisation of the (1×8) state and destabilisation of the (1×9) state continues. This leads to an ultimate barrier increase 302% for the (1×8) state, and barrier decrease of 91% for the (1×9) state.

Conclusions and Outlook
In this work, a triangular lattice model is used to evaluate the free energy of postbuckled states of elastic thin shells. This is implemented in efficient energy-minimisation and path finding algorithms in order to fully describe the buckling landscapes. Here, we have demonstrated this for the complex problem of buckling of fixed-end cylindrical shells, subject to axial compressive strains. To begin with, we surveyed the free energy minima, observing unbuckled, singly dimpled, and multiply dimpled states whose stabilities were evaluated for different aspect ratios and compressive strains. We then systematically used the string method to connect pairs of minima within the same cylindrical system in order to find the minimum energy pathways and transition states between these states. This enabled a global description of the buckling landscape: in which a simple funnel-shaped landscape became complex and glassy when increasing the aspect ratio, or featured many deep states when increasing the compressive strain. We then finally introduced the landscape biasing method to control the stability of targeted features of the landscape, in order to design structures with improved resistance to lateral forces.
Overall, by being able to both survey the free energy landscape and design specific transition modes through landscape biasing, we may now design dynamic buckling responses for diverse applications, ranging from energy harvesting devices to complex morphable materials.
One important consideration we highlight for future work is that of the role of imperfections in buckling responses, a significant concern in real-world applications. The ability to generalise our model to consider shapes other than the perfect cylinder, as well as including diverse elastic modulations and boundary conditions, lead us to emphasise the applicability of this model to studying the impact of a large range of different geometric or elastic imperfections on the buckling landscape.

Discretisation and free energy
To evaluate the free energy of an arbitrary thin shell (or composite of thin shells), we discretise the surface into a triangulated mesh of nodes, defining a set of neighbouring nodes and a set of neighbouring planes, in a manner based on [50] although other similar methods have also been reported, for example [51]. The local form of this discretisation is shown in Fig. 6. In this, neighbouring nodes i and j are connected by an extensional spring of equilibrium bond length r 0 ij and elastic constant k stretch ij . Neighbouring planes α and β are connected by an angular spring of equilibrium angle θ 0 αβ and elastic constant k bend αβ . In general, as in our triangulation scheme, r 0 ij and θ 0 αβ are non-uniform across the lattice. The discretisation of the shell into a set of extensional and angular springs allows the total free energy  Figure 6: The thin shell discretisation scheme. The nodes are indicated with black circles, in which nodes i and j are separated by a distance r ij . The planes are indicated with coloured triangles, in which the dihedral angle between planes α and β, θ αβ , is shown as the angle between the respective normal vectorsn α andn β .
to be decomposed into a sum of stretching and bending energies such that generally, where r ij is the separation distance between nodes i and j; θ αβ is the dihedral angle between planes α and β, defined as the angle between the respective normal vectorsn α and n β . Throughout this work, we report the nondimensionalised free energy E r = E/k bend ref , where k bend ref is a reference dihedral elastic constant. For cylinders of uniform elasticity, we define k bend ij = k bend ref . Furthermore, the bond lengths are nondimensioanlised by expressing r ij and r 0 ij relative to a reference length scale R 0 , which we choose to be the cylinder radius.
The single parameter defining the cylinder's elastic behaviour then becomes the control ratio k stretch ref R 2 0 /k bend which unless otherwise stated we fix at 2.5 × 10 5 . Through comparison with continuum elastic theory [50], in terms of Young's modulus Y , plate thickness t, and Poisson ratio ν . Hence, the control ratio is given by 9 To demonstrate the physical significance of our prescribed control ratio of 2.5 × 10 5 , if we choose a Poisson ratio appropriate for aluminium, ν = 0.3, the resulting ratio R 0 /t=247 is similar to that of aluminium drinks cans (R 0 /t ≈ 300). The cylinder radius R 0 is fixed throughout, such that to change the uncompressed aspect ratio A 0 , only the length L 0 is varied. In order to accurately calculate the free energy while balancing computational cost, the number of nodes must be sufficient to capture the deformation profiles of single dimples, the length scale of which depends on A 0 and the control ratio. For the cylinders studied here, ≈ 10 4 nodes per cylinder are required (an illustrative resolution test is shown in Supplementary Note 6). Our triangulated lattice model is also validated against ABAQUS/Explicit commercial software [52], shown in Supplementary Note 7.

Minimisation and path finding
The L-BFGS algorithm [53,54] is employed to efficiently minimise the free energy with respect to the large number of degrees of freedom (O(10 4 −10 5 )). For this, the total free energy is required as well as the derivatives of E with respect to each degree of freedom (the x, y and z coordinates of each node). By setting selected derivatives to zero prior to minimisation, we can constrain specific node positions. Here, we fix the x and y coordinates of the nodes which cap each end of the cylinder to the uncompressed configuration, forbidding deformation or relative rotation of the ends. By choosing the z coordinates at which to fix these nodes, we can achieve the desired cylinder end shortening. An example minimisation convergence plot is shown in Supplementary Note 8.
To simulate local probe experiments, in addition to fixing the end caps we also fix the position of a single node in the centre of the cylinder (thus mimicking a point probe). This point is moved radially inwards by a small increment and the free energy minimised. This increment-minimisation procedure is repeated until the entire pathway from the unbuckled state to a second minimum has been obtained.
In the free energy minimum survey, we access the many different dimpled states by performing a basin hopping step prior to each minimisation [42]. To perform this step, we begin with the unbuckled cylinder, and make a random number of trial dimples to the initial node coordinates. Each trial dimple consists of a paraboloidic indentation radially into the cylinder, in which the indentation depth is allowed to vary up to R 0 /2. The minimum energy pathways (MEPs) between any two minima of equal end shortening are found using the string method [55], which we augment for use with highdimensional systems. To begin with, the end points are maximally aligned through rotation and reflection of the displacement fields. An initial string of 30 images is then formed which interpolates the coordinates of the two end points. One iteration of the algorithm consists of evolving each image in the downhill direction, then re-interpolating the images along the new string. The Euler and Runge-Kutta methods used in [55] are however highly inefficient for the high-dimensional energy landscape considered here. Instead, we use 300 L-BFGS steps to rapidly converge the string to the MEP. A simple linear re-interpolation scheme is used, with the image density concentrated at the highest energy points along the string. This process is iterated until the E r of the highest energy point along the string changes by less than 10 −6 from the previous iteration. If intermediate minima exist along the pathway, a separate string is evolved for each, such that each pathway connects two minima via a single transition state. The Euler method is employed in the final stage to fine-tune the pathway, such that convergence is achieved when the RMS distance between the strings is less than 10 −6 . The transition state is then fine-tuned using the climbing string method with Euler steps [44], finishing once the RMS gradient is reduced below 10 − 5. Repeating the string algorithm to connect multiple end points forms a network of connected minima.
In order to show the general validity of this model, we further apply it to analyze the energy landscapes of the buckling of spherical caps in Supplementary Note 9. Our model has similar accuracy as the finite element model implemented in ABAQUS, and successfully captures the stable axisymmetrically inverted configuration of the spherical cap [56]. The analysis is suitably rich that we reserve further discussion for another publication. Across the range of aspect ratios A 0 and shortening ratios λ, it was observed in Fig. 1d (main text) that the singly dimpled state was never the global free energy minimum. For cylinders of low aspect ratio, this observation is explained by the high localisation in both the axial and azimuthal directions, such that multiples dimples may form on the structure which have sufficient spacial separation to not interact. If a single dimple reduces the overall free energy of the system by exchanging stretching contributions for bending contributions, then a second spatially separated dimple reduces the free energy further by an equal amount.
Thus, for any system with physical or elastic properties that enable multiple non-interacting dimples to form, the single dimple may never be the global energy minimum.
We next consider the alternative case where multiple dimples interact strongly, an extreme example of which is illustrated in Supplementary Figure 1. As will be shown in Supplementary Figure 2 and Supplementary Figure 4, spatial localisation of a dimple is reduced on increasing A 0 and reducing the elastic control ratio, allowing for multiple dimples to interact over the shell. We therefore aim to probe the single dimple stability at extreme values of A 0 and elastic control ratios. The limit of this investigation is shown in Supplementary Figure 1, where A 0 =10, λ=0.990, and the elastic control ratio = 2.5×10 3 . For reference, a material with Poisson ratio = 0.3 and elastic control ratio = 2.5×10 3 yields a shell with radius to thickness ratio R 0 /t=25. The singly dimpled state shown in Supplementary Figure 1a exhibits extreme deformation, with a maximum radial displacement of 75% of R 0 , and a reduced free energy of E r =34.8. The free energy is still reduced on forming a second dimple however, shown in Supplementary  Figure 1b, where E r =31.9. Overall, over a large range of physically relevant conditions, the singly dimpled state is never the global free energy minimum.

Self-interaction
The low-λ single dimple stability limit, shown in Fig. 1d (main text), exhibits a non-monotonic variation with A 0 when A 0 ≥5. In Supplementary Figure 2, we plot the radial displacement of the single dimple about the centre of the cylinder at the low-λ stability limit. It is observed that for A 0 =2, the dimple is highly localised in the azimuthal direction, and does not extend around the full circumference. This is shown in Supplementary Figure 2(inset) as the radial deformation tending to zero far from the dimple centre. This angular localisation decreases however as A 0 is increased, so that for A 0 ≥5, the radial displacement far from the dimple centre is of the order of 1% of R 0 . Thus, self-interaction effects modify the single dimple stability for A 0 ≥5, leading to the non-monotonic low-λ stability limit on increasing A 0 .

Supplementary note 2: Further analyses of the local probe technique
Accuracy of the unbuckled-single dimple local probe In Tables 1 and 2, we compare the transition state energy E TS from the MEP, with that obtained via the local probe technique for the unbuckled-single dimpled transition. In every case tested, the local probe achieves highly accurate estimates of the MEP transition state energy with a small percentage difference between the two (typically of the order of 10 −3 %). This is consistent with the typical accuracy of the local probe: successive increments change the energy by 10 −5 − 10 −4 close to the transition state.

Efficacy of multiple local probes
We have shown that the unbuckled-single dimple minimum energy barrier can be accurately accessed through the local probing technique for A 0 = 0.8, λ = 0.9986. In Supplementary Figure 3, we investigate whether local probing can access more complex pathways. This is tested via attempting to reproduce the sets of transitions which occur in the multi-step unbuckled-(1×9) pathways. This is attempted for the selected paths A, B, and C. We focus on Path A initially, reproduce in Supplementary Figure 3a.
To begin the local probe setup, we must estimate where the dimples will form without prior knowledge of the MEP. It is reasonable to assume the (1×9) state is formed by nine separate dimpling transitions. which are equally spaced around the circumference of the cylinder centre. These locations are labelled in Supplementary Figure 3b. A number of choices exist with how to proceed with locally probing at nine locations, regarding in particular: whether the probes are applied individually or in unison, the depth each probe should be tested to, and whether probes should be removed or remain in place after a dimpling event. Here, we choose to test the pathway as simply as possible, in which: (1) A single probe is initiated at radius R = R O ; (2) this probe is incremented radially inwards, up to the point where the energy of the system begins to reduce; (3) the system is then relaxed with the probe removed to form a new (dimpled) equilibrium configuration; (4) The next single probe is applied, and steps (1-3) repeated. By selecting the order of the probes, we may access different dimpling pathways. To compare as closely as possible to the MEPs here, we observe the dimpling order in paths A, B, and C, and perform the local probes in this order. The dimpling order is shown for paths A, B, and C, in Supplementary Figure 3a,c, and d respectively. The parenthesised pair of dimples in these orders indicate both dimples form in a single transition, but in general this would not be known a priori.
For path A in Supplementary Figure 3a, each of the eight separate transitions are labelled Ai-Aviii, and shown individually magnified. It can be seen that transition Ai is well reproduced in both energy and location by the local probe. As the series of transitions proceeds up to Avii however, although the local probe is able to accurately capture the transition state energy, the transition state location is obtained with decreasing accuracy. This can be rationalised by observing the movie of path A in Supplementary Movie 1. It can be seen that the dimples shift slightly in their location as the transition proceeds, which cannot be accommodated for in the fixed local probe simulation. In Aviii, the local probing technique fails to capture the MEP, as two dimples form in the same transition in the MEP, but not the probing simulation. As can be seen at the end of Supplementary Movie 1, the start of the transition sees one dimple partially form, then the other grows to match it, before both grow together to complete the transition. For path B in Supplementary Figure 3c, the nine-point probing scheme departs from the MEP after only two transitions. This is because for this sequence of dimples, a probe applied at one point causes a large shift in the previouslyformed dimple locations. Successive probing therefore continues to deviate further from the MEP. The same observation is made later in the pathway for path C in Supplementary Figure 3d.
Overall, this multiple local probing scheme is able to closely estimate the minimum energy barrier, but only when: (1) The dimpling event is single-dimpled; (2) the probe does not cause a large shift in previously-formed dimples; (3) the dimple formed by the probe does not undergo a large shift on removal of the probe. This latter point would be particularly challenging to maintain for irregularly dimpled structures, where a best estimate of the probe location may not be readily available.

Supplementary note 3: Elastic property effects on the free energy landscape
In Fig. 3c (main text), it was observed that for cylinders of high aspect ratio and thin shell thickness, at compression ratios where the unbuckled, singly dimpled and multiply dimpled states coexisted, the energy landscape became highly complex. This was observed for the cylinder of aspect ratio A 0 = 3.0, elastic control ratio 2.5×10 5 , and compression ratio λ=0.999. This complexity was manifest both in the number of minima, O(10 3 ), and multiple energy scales over which the landscape was 'rough'. Here, we observe the effect of decreasing the elastic control ratio (equivalent to increasing the shell thickness) on the landscape complexity. In Supplementary Figure 4, we illustrate the energy landscape for a cylinder of the same aspect ratio, A 0 = 3.0, and use a compression ratio where the unbuckled, singly dimpled and multiply dimpled states also coexist, λ = 0.996, but decrease the elastic control ratio to 2.5×10 4 . This decrease in the elastic control ratio, equivalent to a thickness increase of ≈ 3×, has a marked, simplifying affect on the energy landscape. The landscape is transformed from being rough and glassy (O(10 3 ) minima, bin population variance = 93), to being funnel-shaped with several deep states(O(10 1 ) minima, bin population variance = 0.9). The global minimum is a (2 × 4) diamond pattern, but a deep (1 × 4) state also exists. At low energies the landscape is dominated by systems with 7 to 8 dimples. The dominating reason for this pronounced simplification is due to the size of each dimple relative to the cylinder area: whereas 8 dimples pack round the cylinder at the global minimum here, when the elastic control ratio was 2.5×10 5 , 12 were able to do so. The number of dimple arrangements possible on the cylinder here is therefore markedly reduced. It is also instructive to compare Fig. 3b (main text) where A 0 = 0.8, elastic control ratio 2.5×10 5 , and compression ratio λ=0.998, with the disconnectivity graph here. In Fig.  3b (main text), the strong confinement introduced by the dimples interacting with the fixed ends caused numerous deep states to exist, as well as for the majority of states to occur in the low-energy region of the landscape. Both of these observations are also made here in Supplementary  Figure 4, suggesting that by decreasing the elastic control ratio, the dimples are confined by the fixed ends. This is in contrast to the cylinder with A 0 = 3.0, elastic control ratio 2.5×10 5 , in which the dimples were weakly confined, giving rise to the glassy landscape.      Table 3. The deep states exhibited by the A 0 = 0.8, λ=0.998; and A 0 = 3.0, λ=0.996 systems mean the range of barrier energies for both span approximately 3 decades. This is compared to the 2 decade range in energy barriers for the systems lacking deep states: A 0 = 0.8, λ=0.9986; and A 0 = 3.0, λ=0.999.
Next, we partition each landscape minimum energy range into 100 bins of equal width, and count the associated number of states within each. These state density distributions are plotted in Supplementary Figure 5 (lower panels). Uniquely in our study, for each dimpling configuration in the A 0 = 3.0, λ=0.999 system, a number of similarly structured minima exist with small energy barriers to their interconversion. For this reason, we clustered states with the same number of dimples and energy barrier less than 10 −3 . The uniformity of each distribution is quantified by calculating the variance in the bin populations, listed in Table 3. The glassy landscape of the A 0 = 3.0, λ=0.999 system is readily distinguishable from the large population variance, 93. The two roughness scales of the energy landscape also become apparent here: within clusters, barriers of less than 10 −3 are present, associated with local movement of dimples in the same configuration, whereas the typical energy barrier between clusters is shown in Supplementary Figure 5c  The sequence of steps required for landscape biasing is shown in Supplementary Figure 6. To begin with, a reference state is chosen which undergoes the target buckling 5 response. The example reference shown in Supplementary  Figure 6a is the unbuckled state for the A 0 = 0.8, λ=0.9986 system, with the example biasing target shown as the 0-1 dimple transition state in Supplementary Figure 6b.
Next, the elastic potential energy associated with each node, i, in the triangulated mesh is calculated for both the biasing target, E target (i), and the reference E ref (i). The difference, is shown for the biasing example in Supplementary Figure  6c.
Numerous choices now exist in using the local energy difference to inform the local thickening. We now demonstrate two example methods: one was used to bias against the 0-1 transition state, and the other for the (1×8) minimum.
To bias against the 0-1 transition state, ∆E(i) is shifted by the minimum value of {∆E(i)}, then normalised, so that the weighting of element i is expressed: . ( The local weights are then used to simulate a local thickening of the elastic shell, so that the ratio of the new thickness of element i, t i , to the local thickness of the reference state, where α is the biasing amplitude, and represents the prescribed fractional increase in mass of the shell, if we make the reasonable assumption that the mass increases in proportion to the simulated shell thickness. For the 0-1 transition state example, we make a preliminary step in which the ∆E(i) are averaged about the circumference of the cylinder. The resultant local thickening of this treatment, using a biasing amplitude α = 0.01, yields the thickening map illustrated in Supplementary Figure 6d. The inset shows the z-variation in the local thickening profile.
To bias for the (1×8) minimum, the local weighting w i is obtained via shifting ∆E(i) with respect to the average value of {∆E(i)}, and rescaled by the maximum value of {|∆E(i)|}, such that, The local thickening ratio, t i /t o is then expressed, This treatment ensures that the total mass increase of the biased cylinder is zero, as i w i = 0, and the magnitude of the greatest local fractional change in thickness is the biasing amplitude α. The final step is to transform the local elastic constants to effectively simulate the local thickening. The new local stretching constant between nodes i and j, k stretch ij , relative to the local reference constant k stretch Similarly, the new local dihedral bending constant associated with nodes i, j, k, and l, k bend ijkl , relative to the local The example result of the local thickening treatment for the 0-1 transition state is shown in Supplementary Figure  6e: the transition state has been shifted off-centre, increasing the minimum energy barrier.

Interpolated biasing
Supplementary Figure 7: Impact of the biasing localisation on the displacement of the dimple C from the centre of the cylinder length (red data, left axis), and on the energy barrier ∆E B (black data, right axis).
Upon biasing against the unbuckled-single dimple transition state for the A 0 =0.8, λ=0.9986 system, the biased transition state is observed to form off-centre. This displacement is caused by the energetic penalty associated with deforming the locally thickened central region of the cylinder. As it is also energetically unfavourable to deform the cylinder close to the simply supported ends, the dimple is expected to form at a location C between 0 and L/2 from the centre of the cylinder. We test the hypothesis that symmetry breaking is caused by local thickening about the centre with the following set of simulations (all at the prescribed 1% mass increase). We begin by testing the energy barrier ∆E B on the uniformly thickened cylinder, which we assign an 'interpolation fraction' of 0. We then obtain the energy barriers a series of cylinders with interpolated thickening fields between the uniform cylinder, and the fully biased cylinder (which we assign an interpolation fraction of 1). The dimple position and energy barrier of these interpolations are shown in Supplementary Figure 7.
As anticipated, the more localised the biasing is to the centre, the larger the energetic penalty for deforming the centre, so the further away the dimple forms. However, interestingly there is an optimum interpolation fraction between local and uniform, at approximately 0.6. We can rationalise this observation by considering that to maximise the barrier to dimpling, the interpolation fraction should be large enough to force the dimple off-centre, but not so large that the dimple forms in a region of low thickness. In Supplementary Figure 8, the influence of the triangulated mesh node density on the computed energy is tested. To perform the test, a multiply dimpled state on the A 0 =0.8, λ=0.9986 cylinder was chosen. The energy for this state was then calculated using different numbers of nodes, N nodes . Supplementary Figure 8 shows for N nodes > 10 4 (filled data), the accuracy of the computed energy cannot be further improved. By averaging the energies of systems with N nodes > 10 4 , we show that all meshes produce energies consistent to within ≈1% of this average. The resolution we choose to perform further simulations at is highlighted with a blue circle, balancing accuracy with computational efficiency. To verify the triangular lattice model, we first conducted simulations of indenting a cylindrical shell under axial compression [1] using the triangular lattice model and ABAQUS/Explicit [2]. The cylindrical shell has radius R 0 = 28.6 mm, length L 0 = 107 mm, and thickness t=0.104 mm. It is made of linear elastic material with Young's modulus Y = 71000 N mm −2 and Poisson's ratio ν = 0.3. The indenter is a rigid sphere with diameter D = 4.7 mm. The cylindrical shell is modelled with 4node shell element (S4R) in ABAQUS. The shell is first compressed axially by assigning a constant velocity at the bottom and top edges. This corresponds to the "simply supported" boundary conditions for thin shells. After a targeted axial compressive strain is reached, the velocities of the top and bottom edges are set to zero. The indenter is then approached to the cylinder with a constant velocity perpendicular to the shell surface to generate local probe deformation. (Supplementary Figure 9a). Different loading rates are checked to make sure it is a quasi-static process. The indenting force-displacement curves of the triangular lattice model and the shell element in ABAQUS are in agreement, as shown in Supplementary Figure 9b. The displacement maps of the dimpled structures are also com-pared in Supplementary Figure 9c We then simulated bending of a slender cylindrical shell under an axial compressive strain slightly above Euler's buckling load of a beam using the same boundary conditions. The comparison between our triangulated lattice model and ABAQUS is presented in Supplementary Figure 10. Although there are no constraints for shell edge rotating, the bottom and top surfaces remain flat due to the finite size of the cross section. Therefore we need to apply the "clamped boundary" in the framework of beam theory. The parameters for cylindrical shell were chosen to be R 0 =28.6 mm, L 0 =858 mm, t=3.575 mm, Y =71000 N mm −2 and ν=0.3. This makes sure its global (beam bending) buckling strain global cr = π 2 I A(0.5L) 2 ≈ 0.0220 [3] is much smaller than its local (plate bending) buckling strain [4], where I and A are the second area moment and area of cross section. The shell is first compressed until its axial strain reaches 0.024, slightly above global cr . After applying a small indenting force to initiate global buckling, the cylinder relaxes to an equilibrium post-buckling shape, shown in Supplementary Figure 10a and b. The lateral deflection fields of both models are in agreement, as shown in 10a and b, with a difference of less than 1% in maximum values.

Supplementary note 6: Resolution tests
Step number Supplementary Figure 11: Example convergence of a perturbed cylinder to a buckled state using the L-BFGS algorithm.

Supplementary note 8: L-BFGS convergence
We show an example convergence plot in Supplementary Figure 11 for the minimisation of a perturbed cylinder to a buckled state using the L-BFGS algorithm. Typically, convergence occurs when the energy change on successive steps is less than 10-11% of the energy of a minimum.
Supplementary note 9: Generalising the model In the continuum limit, it has been shown that the effective bending energy of the perfect triangular lattice can be expressed as [5,6], where H is the mean curvature, and K is the Gaussian Curvature. For a continuum sheet, the bending energy can be written as [6,7], where κ =  Figure 12: a Parameters and boundary conditions in the spherical cap problem. The parameters chosen in our study are Y =71000 MPa, ν=-0.3, t=0.7984 mm, R=151 mm and α=0.2 rad. b Indenter force-apex position curves. c Free energy profile along the minimum energy path (MEP) connecting the initial state A and inverted state C and passing through the transition state B, compared with the free energy profile along the probe paths. d Deformation profiles of the cylindrical cap corresponding to the three states on the MEP. By comparing Eqs. (8) and (9), the effective Poisson's ratio of the triangular lattice is ν bending = −1/3. It should be pointed out that the in-plane stretching energy and outof-plane bending energy in the triangular lattice model are decoupled. Therefore the effective Poisson's ratio derived from bending energy has no effect on the Poisson's ratio in the in-plane deformation, which has been shown to be ν stretching =1/3 [8].
To show our method can be generalized to other geometries, we consider the buckling of a spherical cap shell shown in Supplementary Figure 12a. The undeformed zeroenergy shape of the shell is cut from a sphere with radius R, Young's modulus Y , Poisson's ratio ν and thickness t, with angle α between the pole and the edge, represented by the dark blue arc. Its edge is free to rotate but restricted to move within the z=0 plane (u z =0), and its apex is fixed in the x − y direction (u x =u y =0). For a spherical cap with large enough R/t, when probed with an indenter moving down along the z axis, the shell will snap to an inverted local minimum, represented by the light blue arc, once the indenter displacement reaches a critical value [9].
For the triangular lattice model, we mesh the spherical cap with arbitrarily shaped triangles instead of equilateral triangles in cylindrical shells, as spherical shells cannot be isometrically mapped to simple planar shapes. The corresponding bending energy elastic constants for general tri-angulated meshes are [10]: where A α ,A β are the triangular areas and l αβ is the length of the shared edge of two triangles. The equilibrium values of the dihedral angle θ 0 αβ are also computed based on the initial coordinates of the nodes on the spherical cap. To verify the results of triangular lattice model, we also conduct finite element (FE) simulations using shell elements in ABAQUS. It should also be pointed out that we take Poisson's ratio ν to be -0.3 in the FE model to match the ratio between Gaussian and bending rigidities in the triangular lattice model. This issue is not a concern in the main text as Gaussian rigidity does not play a significant role in the buckling of cylindrical shells, whose Gaussian curvature is zero. But in our preliminary indenting simulations of spherical caps, the force-displacement curves vary significantly after choosing different Poisson's ratios in the continuum model. In principle, the in-plane stretching energy will be also changed when varying Poisson's ratios in FE simulation, as the bending and stretching energy are coupled in the continuum model. However, the in-plane stiffness has little influence in the current study because buckling of spherical caps is dominated by bending energy. The parameters chosen in this study are Y =71000 N mm −2 , ν =-0.3, t=0.7984 mm, R=151 mm and α=0.2 rad, so that the spherical cap has an inverted local minimum [9]. We first simulate probing the shell with a hemispherical indenter modeled with triangular lattice and finite element in ABAQUS. Supplementary Figure 12b compares the indenter force apex position z a curve of the two methods, and marks in black star the undeformed local minimum A, transition state B and inverted local minimum C. Supplementary Figure 12c compares the free energy profiles along the probe paths and the minimum energy path (MEP) connecting A and C obtained by string method. Similar to the unbuckled single dimple transition of cylindrical shells, probing meets the MEP at the transition state B, but does not access the MEP generally. Supplementary Figure 12d compares the deformation profiles of the three states obtained by probing and string method. This preliminary work demonstrates that the methodology of integrating triangular lattice model and string method can capture the snap-through buckling of a spherical cap shell. Future work needs to be done to explore the energy barriers and MEP for more complicated buckling states, such as non-axisymmetric modes [9]. In addition, the current triangular lattice has a fixed ratio between Gaussian and bending rigidities, corresponding to a Poisson's ratio ν=-0.3. This is also understandable because there is only one parameter associated with bending in the current triangular lattice model, which cannot produce two different bending rigidities independently. It will be of great interests in developing a more general lattice model capable of describing materials with a wide range of Poisson's ratios.