Bespoke extensional elasticity through helical lattice systems

Nonlinear structural behaviour offers a richness of response that cannot be replicated within a traditional linear design paradigm. However, designing robust and reliable nonlinearity remains a challenge, in part, due to the difficulty in describing the behaviour of nonlinear systems in an intuitive manner. Here, we present an approach that overcomes this difficulty by constructing an effectively one-dimensional system that can be tuned to produce bespoke nonlinear responses in a systematic and understandable manner. Specifically, given a continuous energy function E and a tolerance ϵ > 0, we construct a system whose energy is approximately E up to an additive constant, with L∞-error no more that ϵ. The system is composed of helical lattices that act as one-dimensional nonlinear springs in parallel. We demonstrate that the energy of the system can approximate any polynomial and, thus, by Weierstrass approximation theorem, any continuous function. We implement an algorithm to tune the geometry, stiffness and pre-strain of each lattice to obtain the desired system behaviour systematically. Examples are provided to show the richness of the design space and highlight how the system can exhibit increasingly complex behaviours including tailored deformation-dependent stiffness, snap-through buckling and multi-stability.


Introduction
Historically, elastic nonlinearities have been avoided by structural engineers. This is often due to the need to simplify the design process and limit the sensitivity  to design parameters in order to create robust, repeatable behaviour. Moreover, elastic nonlinearities are often synonymous with sudden catastrophic failure. However, recent advances in material construction processes, e.g. three-dimensional printing, and the availability of advanced computational/experimental analysis methods have permitted a relaxation of this design conservatism [1][2][3][4][5]. In particular, computation has permitted the probing of complex nonlinear design spaces with resultant improvements such as reduced structural weight, increased robustness or additional functionality [2,6].
In general, the potential for multifunctional designs to reduce system part count is particularly attractive from the perspectives of structural weight and minimal design philosophy. However, a trade-off ensues. Conventional components are stiff and movement is achieved through mechanical joints, such as hinges. To eliminate mechanical joints but retain operational loadbearing capacity, greater control of stiffness and form is required at both global and local structural levels [7]. In this context, nonlinear elastic deformations are an attractive means for robust and repeatable shape and stiffness adaptation. The desire to explore and exploit nonlinearity has motivated more complex descriptions of the underlying mechanics. For instance, with structures capable of large deformations, geometrically nonlinear stress analysis is often involved in the design process [8][9][10][11].
An emerging and somewhat more challenging task than describing the nonlinear response of a structure under loading is its inverse, i.e. designing a structural topology and material system to exhibit a specific nonlinear response. The need for bespoke elastic responses has arisen in many fields including energy absorption [12], robotic and MEMS actuation [8,13,14], extreme thermoelastic devices [15] and flow regulation valves [6,16]. The strength of the nonlinearity in these bespoke responses is often limited so it can be realized and manufactured. In cases where the nonlinearity is unconstrained, a discussion of how such responses can be realized with current design and manufacturing technologies is often omitted. In the instances where physical systems have been fabricated, they are of little use outside of the intended application [12]. Therefore, a universal design framework able to achieve an arbitrary nonlinear response in a systematic manner would be immensely valuable to the wider engineering community.
Bespoke nonlinear elasticity has already been investigated by some authors. Typically, the focus has been on finding and coupling suitable problem parametrizations with nonlinear finiteelement analysis (FEA) solvers and optimization algorithms [8][9][10][11]. Due to the complexity of the underlying structural models, the process of tailoring nonlinear responses has often been reported as computationally expensive. As an example, Jutte [17] successfully used topology optimization, coupled with a nonlinear FEA solver, on a two-dimensional network of geometrically nonlinear curved beams. The level of achievable nonlinearities was restricted though, with buckling or snap-through behaviours excluded. A variety of nonlinearities were observed, including constant and nonlinear strain hardening and softening curves.
In the context of nonlinear lattice structures, Pirrera et al. [18] inspired by the behaviour of the virus bacteriophage T4, explored the design space of geometrically frustrated helical lattices (figure 1d) demonstrating robust, elastic multi-stability. Similar structures, with a focus on the handedness of their auxetic properties, were investigated by Lipton et al. [19]. Also of note is the work on gridshells by Baek et al. [20,21], where the authors investigate the shape and rigidity of shell-like structures arising from the buckling of initially planar grids of rods.

(a) Main contributions of this article
To the authors' knowledge, there are currently no tools capable of designing a one-dimensional system to exhibit arbitrary nonlinear extensional responses. We address this deficiency by presenting a system that can be constructed as an assembly of well-understood nonlinear, springlike 'units'. The units' behaviour is described by means of a simple analytical model, which provides insight into the system's mechanics while requiring minimal computational effort for design. As units are added to the assembly, new aspects of behaviour may be observed. In essence, as we increase the complexity of the system, richer responses can be achieved without increasing Sketch of a typical lattice comprising four helical strips (i.e. N = 2) on the surface of its constraining cylinder and its planar representation. Adapted from [18].
the complexity of the description of the system. Following this mechanism, we are able to design nonlinear spring-like systems that can achieve a plethora of desirable response characteristics.
In the present work, the fundamental units are helical lattices like that illustrated in figure 1. This structure was first presented in [18], where robust nonlinear features, such as multi-stability, were observed for various ranges of lattice geometry, stiffness and pre-strain. This lattice acts as an effective one-dimensional spring, the behaviour of which may be described succinctly using its energetic profile under extension. Here, multiple helical lattices are coupled through rigid connections, in arrangements similar to those presented in [19] whilst preserving the kinematic assumption outlined in [18]. We demonstrate the appearance of new response characteristics that are unobtainable via a single lattice in isolation. Besides, further behaviour continues to emerge with additional geometrically distinct lattices.
Mathematically, linking lattices changes the system's energy in a manner that resembles a series expansion. We demonstrate that this offers sufficient design freedom to approximate any polynomial energy, and thus to approximate any continuous energy to any desired accuracy in L ∞ . In particular, we can prescribe the positions and stability of the system's equilibria, as well as its stiffness as a function of extension. In short, we can design a bespoke nonlinear elastic response. We present a systematic approach to obtaining the design parameters for individual lattices and the overall system assembly requirements.

(b) Outline
The remainder of paper proceeds as follows. In § §2 and 3, we first recall and then recast the results of Pirrera et al. [18] for a single lattice. In §4, we extended these results to coupled systems of multiple lattices. In §5, we prepare the ground for the main result which is then proved in §6. Ancillary results are presented in §7. Several examples of engineered behaviours illustrating the robustness of the design space and the algorithm are presented in §8. Finally, conclusion is drawn in §9.

Kinematics and elastic energy of single lattices
Let us consider a generic helical lattice-as depicted in figure 1a- [18] developed an analytical formulation to describe the mechanics of such a system, based on the elastic energy of the helical strips. For the reader's convenience, we repeat the kinematic analysis [15] here.
The lattice can extend or compress as a result of loading at its extremities and along the longitudinal axis (i.e. the X-direction). As the lattice deforms, we assume [15,18] that: (i) all helices remain on a single cylinder (whose radius and length can change); (ii) the helical strips are hinged where they overlap allowing scissoring motion only; (iii) the change in length of the helical strips is negligible; and (iv) the helical strips are sufficiently slender that the dominant contribution to the energy arises from bending and twisting with respect to the strips' principal dimension.
Furthermore, consider the '+' strip that passes through the origin in figure 1b. In the ζ -plane, this strip has the parametrization (ζ 1 , ζ 2 ) = s(sin(θ ), cos(θ)), where s is the arc-length parameter and θ is illustrated in figure 1c. From the assumptions above, this strip maps to the helix in the (X, Y, Z) coordinate system used in figure 1a, where a hat is used to indicate unit vectors along these axes. Q(φ) ∈ SO(3) is a rotation aboutX through an angle φ. The unit tangent t(s), normal n(s) and binormal b(s) are defined as where dash indicates differentiation with respect to the arc-length parameter, s. The unit cell extension, h ∈ [0, l], (2.2) as shown in figure 1c, is used to uniquely describe the state of the system. Accordingly, the curvature κ x , and torsion κ xy , can be defined in terms of the unit cell extension, such that and These definitions hold for the '−' strip as well, except that the sign of the torsion is reversed. Since we have assumed that the strips' behaviour is dominated by bending and twisting deformations, i.e. the stretching energy can be neglected, as in [15,18], we obtain a quadratic bending energy of the form where the curvature components κ and associated reduced bending stiffness matrix d are described by Classical Laminate Theory (CLT) [22] capturing stiffness anisotropy. The use of CLT is for notational convenience alone and does not imply that the system must be composed of composite strips. For example, anisotropy could be designed into the system by exploiting cross-sectional geometry or functionally graded materials.  Figure 1 illustrates the lattice geometry schematically, its development into a plane, and the representative unit cell upon which the lattice behaviour can be modelled. Herein, we limit right-and left-handed helices-denoted by subscripts · + and · − , respectively-to have equal and opposite pitch, θ . This restriction ensures unit cells are rhombi and prevents the lattice from twisting upon extension [18]. In addition, and conveniently for the current purposes, it creates a bijective correspondence between R and h [18], permitting the lattice extension to be used as the single degree of freedom to characterize the system [23].

Unit cell formulation
(a) Total potential energy of a unit cell With the unit cell representation, it is sufficient to consider the total potential energy of the constitutive cell to describe overall lattice behaviour. For simplicity, but without loss of generality, it is assumed that each right-handed helix has identical dimensions, stiffness and mechanical prestrain and that these properties are constant along the length of the helices. A similar commonality is shared between left-handed helices. Using (2.4), the energy of the representative unit cell can then be written as where (2π l/N) is the cell's side length and w ± the strip's width. The model by Pirrera et al. [18] considers axial, · x , and twist, · xy , components of curvature. Moreover, the total curvature comprises two additive terms: the curvature developed by the helices upon deformation while conforming to the surface of the cylinder, κ ± (h), and the mechanical pre-curvature (tooling), υ ± . Hence, which is included in the energetic formulation with the corresponding terms of the reduced bending stiffness matrix, d ± 3) The curvatures due to lattice extension may be stated using the Weingarten map [24] κ ± (h) = 1 R where, from figure 1 Equation (3.1) now may be expressed in terms of non-dimensional parameters.

(b) Non-dimensional form
To aid analysis, we make the following normalizations, indicated by an overbar: This allows the energy of the representative unit cell, (3.1), to be recast as for h ∈ [0, 1], with For feasible systems, the stiffness coefficients must guarantee the positive definiteness ofd ± and thus we requireφ Although the definition of the strain energy used herein only considers curvatures due to mechanical loading and pre-strain, the description can be readily modified to include additional fields, such as thermal, electromagnetic, piezoelectric, pH and moisture [15]. The inclusion of such effects is straightforward, but beyond the scope of the current analysis. Inclusion of these effects can improve prediction of manufactured performance of similar structural architectures [25,26] and field-dependent curvature changes offer additional avenues to manipulate the energy landscape of the system. In other words, the fields can be used, either in isolation or synergistically, to enhance the tailorability of the lattice. For instance, both the equilibria and their characteristics of stability can be tuned to change with temperature, creating nonlinear cyclic thermal actuators or highly tuned coefficients of thermal expansion beyond what is achievable by traditional materials [15].

Coupled lattice system
Although, as demonstrated in [18], a single helical lattice may be tuned to exhibit many desirable nonlinear behaviours, it does not provide the designer with complete freedom to develop a complete set of smooth energy profiles as a function of lattice extension. To overcome this limitation, it is advantageous to consider the behaviour of multiple lattices coupled together. Through the coupling of unique lattices, new responses are obtainable that cannot be achieved by a single lattice in isolation. Developing this hierarchical system shows how, depending on length scale, the tuning of nonlinear substructures can be considered analogous to the development of architectured or meta-materials (see, for instance [27]).
Then, the normalized extension of the ith lattice,h i (see (3.6a)), can be related to the normalized extension of the system,H, through a parameter ψ i : where, using (4.1),

(ii) Lattice energy
The energy of the ith lattice, ξ i , may be obtained by multiplying the energy of each unit cell, given by (3.7), and the number of unit cells, 2N i M i (cf. figure 2a): where we have used (4.2b). This can be written as and the coefficients a 0i , a 1i , a 2i , b 1i and b 2i are related to the material parametersδ i± ,¯ i± ,ῡ ix± ,ῡ ixy± andφ i through (3.8). Note that 3), the system energy, Ξ , can be written as with coefficients and and functions f ψ , g ψ : [0, 1] → R + given by and The A coefficients add no new behaviour to the lattice system compared to a single lattice since the summation merely replaces one constant with another. Functions of the form B i f ψ i (H) and C i g ψ i (H), however, contribute to the system energy with distinctive terms that a single lattice cannot possess. As we shall see, a generic lattice system will feature behaviour that cannot be reproduced by an individual lattice.

(iv) Stiffness and equilibria
Differentiating (4.7) with respect toH, we obtain Then, from (4.5), the first two derivatives of the normalized energy,Ξ , with respect to extension are The stiffness, equilibria and characteristics of stability of the system can be deduced from these expressions.

(b) Concentric lattices
We can restrict the behaviour of the system by arranging the lattices in a concentric manner, as illustrated in figure 3. As the lattices are now concentric they are not permitted to intersect one another. In other words, when the lattices are numbered outwards from the centre, R i R i+1 . Since, from (3.5), (4.1) and (4.2), the radius of the ith lattice can be written as for i = 1, . . . , I − 1. This is equivalent tō where we have used (4.1a) and (4.2c). Enforcing (4.10) for i = 1, . . . , I − 1, we obtain the range (which might be empty) ofH for which radial self-intersections do not occur to be ⎡

(i) Absence of radial self-intersections
Radial self-intersections do not occur when for all i = 1, . . . , I − 1, since under these conditions the inner maximum on the right of (4.11) is vacuous, and the inner minimum on the left is non-positive. By a similar reasoning radial selfintersections do not occur also when More comprehensively, (4.11) enables a characterization of the conditions on (l i+1 /l i ), (M i+1 /M i ) and (N i+1 /N i ) under which radial self-intersections do not occur. This would allow optimization to be conducted in two stages: first, the continuous design parameters for achieving the desired response can be obtained, assuming maximal range of motion. Second, banding numbers N i , M i and unit cell lengths l i are chosen to avoid radial self-intersections. See §8 for examples.

Taylor expansion of the system energy
In this section, we derive the Taylor expansions of the system energy (4.5) and present related results that are used in §6.
We denote the non-negative integers 0, 1, 2, . . . by Z * , the L ∞ -norm of a continuous function f : 1] |f (x)|, and the smallest/largest singular value of a square matrix T by σ ± (T). and for j ∈ Z * . It is convenient to also define, for J ∈ Z * , and We also list some properties of c 2j , c 2j+1 , C 2J and C 2J+1 which will be needed later: Proof.
(b) Truncated Taylor expansions of f ψ and g ψ Returning to (5.1), the finite truncations of the Taylor expansions of f ψ and g ψ are given by and Note that even though c j < 0 unless j = 0, from (5.1), the polynomials above are positive forH ∈ [0, 1]. By Taylor's theorem, we can truncate the expansions in (5.1). Given > 0, we can find where R 2J P (ψ·), which is the remainder in the Taylor expansion of f ψ to order 2J P , satisfies Then, instead of (5.4b). This implies which is the uniform version of (5.4e).
In §5d, we analyse the system energy using the expansions and bounds above. However, before we do that, we present some properties of the Taylor polynomials in (5.3) that will be required there.

(c) A basis for polynomials
Our first result is that the Taylor polynomials in (5.3) form a basis for the polynomials: To prove the lemma, it is convenient to set, for I ∈ N, which is invertible when ψ i = 0, i = 1, . . . , I; and Now, the polynomials x 2j , j = 0, . . . , J, form a basis for the vector space of even polynomials of order at most 2J. Since both matrices in (5.7a) are invertible it follows that the same is true for the polynomials P 2J (ψ i ·), i = 1, . . . , J + 1.
The proof of the second assertion for odd polynomials is similar and follows from: and the invertibility of the matrices in (5.7b).

(d) Truncated Taylor expansion of the system energy
From the Taylor expansions in (5.4), given J ∈ Z * , we have a Taylor expansion for (4.5): The terms of order 0, 1 and 2 have contributions from A 0 , A 1 and A 2 , respectively, but the higherorder odd and even terms follow the pattern of the terms of order 3 and 4, respectively. From (5.4), the L ∞ -error in (5.9a) is and As our last result for this section, we write the Taylor expansion (5.9) in the P-basis (cf. corollary 5.4 and (5.3)).
For the even coefficients of the polynomial in (5.9b), using (5.6), we obtain, Similarly, for the odd coefficients of the polynomial in (5.9b) we have

Lattice systems approximating a prescribed energy functional
In this section, we prove our main result, namely, that given > 0 and an energy E ∈ C([0, 1], R), a lattice system can be constructed whose energy, is approximately E up to an additive constant, with L ∞ -error no more than . In fact, it suffices to show this for polynomial energies since, from Weierstrass approximation theorem [30], we can find a polynomial E P : In the light of this, we will show that, given > 0 and a polynomial energy E P : [0, 1] → R, a lattice system can be constructed whose energy, Ξ , is approximately E P , up to an additive constant, E 0 , with L ∞ -error no more than ( /2): Then, it would follow from (6.1) and the triangle inequality that max We will prove this assertion in two steps: First, we will show in theorem 6.2 that, given a polynomial energy E P : [0, 1] → R, of the form for some N ∈ Z * , there exists an appropriate choice of lattice system parameters such that the system energy (4.5) differs from E P by no more than ( /2) in L ∞ . Then, we will show in theorem 6.3 that lattice systems with these parameters, with the possible exception of the constant term A 0 in (4.5), can actually be constructed.
Remark 6.1. For convenience, we will allow for the possibility that e 2N+1 = 0 in (6.2). This allows us to assume, with no loss of generality, that the highest power of the polynomial E P is odd.
(a) Approximation of polynomial energies Theorem 6.2. Given > 0 and a polynomial E P , there exist I ∈ N, (6.3a) such that the lattice system energy (4.5) differs from E P by no more than ( /2) in L ∞ .
Proof. From remark 6.1, we may assume that the highest power ofH in E P is 2N + 1, for some N ∈ Z * . Set I = N + 1 and pick D > 0 such that max sup where R 2N P and R 2N P +1 are defined in (5.4a,c); this is possible by remark 5. indeed, we can write it as a (non-unique) combination of P-and monomial-bases. To do so, we would find (non-unique) solutions to the affine equations The existence of solutions to these equations follows from the invertibility of the four matrices in (6.5), see remark 5.1(ii) and remarks following (5.6). The solution space is three-dimensional as can be seen by observing that solutions exist for any A 0 , A 1 , A 2 , and are unique once these are fixed (again, by invertibility of the four matrices that appear in (6.5)). Now choose A j , j = 1, 2, 3, and let and We wish to solve (6.5) under the constraint that Now choose ψ i ∈ (0, 1], i = 1, . . . , I, distinct so as to satisfy (6.7); this is possible since the expressions on the left in (6.7) can be made arbitrarily large by choosing as many ψ i s as necessary to be sufficiently close to each other. Finally, let B i , C i , i = 1, . . . , I solve (6.5). By the construction above (6.6) is satisfied.
The lattice system with the parameters chosen above has energy given by (4.5). By (5.9a,b) and (5.11) and (6.6), the Taylor expansion of this energy to order 2N + 1 is precisely the given polynomial E P . From (5.10), (6.4) and (6.6) the error introduced by the Taylor approximation is bounded by ( /2).

(b) Existence of lattice systems with required energy
Next we show, in theorem 6.3, that a lattice system can be constructed that attains the energy in theorem 6.2, up to an additive constant. We do this by showing that, for each lattice in the system, the material parametersδ ± ,¯ ± ,ῡ x± ,ῡ xy± ,φ and the geometric-material parameter Υ can be chosen so as to yield the coefficients (6.3) required by theorem 6.2, except possibly for A 0 . In fact, we shall show that this is possible even when we impose the additional constraints δ + =δ − , (6.8a) andῡ xy+ =ῡ xy− , (6.8d) on each lattice in the system.
Note that A 0 is omitted from the list of coefficients in (6.3) since the energy of the lattice system in theorem 6.3 is allowed to differ from the energy of the lattice system in theorem 6.2 by an additive constant. Similarly, in the proof below we will omit (3.8a) from the (3.8).
Next, for each lattice, given the coefficients a 1 , a 2 , b 1 , b 2 ∈ R we shall construct a solution to (3.8b-e), which satisfies (3.9).
To do this, we first pickφ ∈ R to satisfȳ ϕ > max (0, −1 − a 2 ) (6.10a) This is possible since the coefficient of the highest power of the polynomial above is positive, guaranteeing that a sufficiently largeφ will satisfy both inequalities. Note that (6.10a) implies (3.9a).

Next we pickδ
It is easy to see that (3.8c,e) are satisfied. Moreover, from (6.10a), (3.9b) is also satisfied. We now calculate¯ 2 and an algebraic calculation shows that this is negative precisely when (6.10b) is satisfied; thus (3.9c) is satisfied. Finally, we pickῡ x± andῡ xy± to satisfy which is simply (3.8b,d) written as a linear system. This has a unique solution since, by (3.9c), the matrix is invertible. This completes the proof.
Note the considerable freedom in the choice of coefficients and geometric and material parameters in the constructions in the proofs of theorems 6.2 and 6.3. In particular, as the computational examples in §8 show, the lattices in the system need not satisfy (6.8). The choice of (distinct) ψ i ∈ (0, 1], i = 1, . . . , I, so as to satisfy (6.7) can be simplified by the following two remarks: Remark 7.1. The bounds in (6.7) can be simplified (but made less tight and therefore more restrictive) using the observation that

(b) Preservation of local extrema
While the construction presented in §6 guarantees a lattice system whose energy is close to a prescribed function, there is no restriction on how the system energy behaves within the tolerance band. In particular, the system energy could eliminate local extremizers of the prescribed energy. However, for many applications, it is desirable that the system energy retain these local extremizers.
We now show that this can be easily incorporated into our framework, provided that the number of local extremizers is finite.
Given be the sets of local minimizers and maximizers of E, respectively. We assume that these sets are finite. For the system energy, Ξ , to be able to distinguish between adjacent local extremizers m ∈ m E and M ∈ M E we need the tolerance to satisfy In other words, we need Thus, Ξ preserves the local extremizers of E in the sense that Ξ has a local minimizer/maximizer between any two adjacent local maximizer/minimizer of E, respectively. Note, however that it is possible that Ξ has more extremizers than E.
The following counterexample shows that the condition that E have only finitely many local extremizers is essential. is continuous (in fact it is absolutely continuous and thus differentiable a.e., cf. e.g. [31, p. 179, 209]). However, it is clear that it has infinitely many local extremizers, only finitely many of which can be distinguished for any tolerance > 0.

(c) Specification of local extremizers
It might be that, rather than approximating an energy, we are concerned only with approximating its extremizers. For example, we might desire that a system have local minimizers at specified values ofH, with the precise energy being unimportant. Suppose we want an energy with local minimizers at 0 m 1 < m 2 < · · · < m K 1, for some E P 0 ∈ R. Then, the local minimizers of E P are precisely m 1 , m 2 , . . . , m K . To see this, note first from (7.5) that m 1 , m 2 , . . . , m K are extremizers of E P . From (7.5a), On the other hand, forH ∈ (m 1 , M 1 ), from (7.5a), F P (H) > 0 since it is the product of the positive factorH − m 1 and 2(K − 1) negative factors. From this, we deduce that m 1 is indeed a local minimizer of E P . It is clear from (7.5) that the same is true for m 2 , . . . , m K . (When K = 1, we choose F P (H) =H − m 1 ; it is clear that the resulting E P has a unique minimizer at m 1 .) We can now approximate the polynomial energy E P as per §6a,b. Note the K degrees of design freedom that results from the ability to choose M k , k = 1, . . . , K − 1 and E P 0 in (7.5b). From the reasoning above, it is clear that an energy with local maximizers at (7.4) can be obtained by changing the sign of (7.5a).

(d) Approximating derivatives of the energy
In reality, one is rarely concerned with prescribing a system's energy. Instead, one typically desires to specify equilibria, reaction forces and stiffnesses. In other words, it is not quite the energy landscape but rather its first and second derivatives that are prescribed. However, as the following simple example illustrates, the derivative of an approximate energy need not approximate the derivative of energy. Consider the energy for some δ > 0, and suppose we wish to approximate it to error > δ. Since approximates this energy to the desired accuracy. From theorems 6.2 and 6.3, we can construct a lattice system whose energy approximates E up to an additive constant. (In this case, the lattice 'system' is trivial since there is only one lattice, i.e. I = 1; see [18,Section 5]). So we have for some constant E 0 . Thus, the derivatives of both the polynomial approximation of the energy and the lattice energy are However, and, when > 1, it is not true that Thus, the derivative of an approximate energy need not approximate the true derivative of the energy. Moreover, this is true even if we allow for different tolerances for the energy and its derivative. To see this consider the following slight modification of example 7.3: Thus, in order to approximate the first or second derivative of the system energy one would have to repeat the analysis in § §5 and 6 for the expressions in (4.8), which analysis we shall present elsewhere. In the meantime, we remark that examples such as the two presented above are pathological and the derivatives of reasonable energies can indeed be approximated. We demonstrate this computationally in §8.

Examples of bespoke system design
We now proceed to demonstrate the ability to tailor the energy landscape by presenting four examples of bespoke nonlinear elastic responses designed as described above: (i) A monostable system with specified nonlinear stiffness (Example 1 in §8c). (ii) A bi-stable system with specified local minima (Example 2 in §8d). (iii) A multi-stable system with specified local minima (Example 3 in §8e). (iii) A multi-stable system with specified snap-through loads (Example 4 in  §8f).
Although the analysis is valid forH ∈ [0, 1], for the computational examples, we choose to limit the extension toH ∈ [0, 0.8] for the following reason: Because the truncated Taylor expansions in §5b are less accurate whenH → 1, the designer can account for this by a priori choosing to work with less than the maximum possible extension. This decreases (in L ∞ ) the remainder terms in (6.4) and consequently allows a larger D to satisfy that inequality. From (6.7), this increases the minimum spacing between the ψs and thus increases the robustness of the design. Here we have pickedH ∈ [0, 0.8] as an illustrative example.

(a) Design
The system is designed as follows. First, for each lattice, the continuous non-dimensional geometric and material variables and the continuous dimensional parameter, Υ , are determined subject to the restrictions (Note that (8.1a) does not apply to ψ 1 , which is not a design variable and is chosen to be 1.) These restrictions, which incorporate (3.9a,b) and (4.4b), ensure that the lattice system is physically feasible. However, the precise numerical values in (8.1) are arbitrary. Note that at this stage the system is independent of length scale and material. The nondimensional variables can be steered towards preferred regions in the design space, e.g. a preference that helices of different handedness exhibit similar bending stiffnesses, desiringφ ∼ 1 or requiring that the pre-strain magnitudes are similar between lattices. However, no preference has been imposed on the following examples; the restrictions in (8.1) are deemed sufficient.
In a second stage, the dimensional geometric variables, N, M and l are recovered to satisfy ψ (determined in the previous step) in accordance with (4.2c), subject to the restrictions Length-scale dependence is introduced through l, and material dependence through the remaining physical variables d 11± and w ± . These are chosen to satisfy the known system parametersφ and Υ for each lattice in accordance with (3.6e) and (4.4a). The material coefficients are restricted to control system feasibility; the unbounded nature of Υ is used to determine appropriate materials and also gives an indication of design feasibility at the interested length scale. A modified Newton's method is used to optimize over the continuous variables (Stage 1) and a genetic algorithm is used to optimize over the discrete variables (Stage 2). It should be noted that the algorithm can also be run in reverse order, i.e. given an envelope on banding numbers or/and material length or/and width restrictions, the lattice material parameters can be tuned to achieve the required system response.

(b) Implementation
Next we highlight some details of the implementation. As demonstrated in §5, the strain energy function of a lattice system can be approximated by a truncated Taylor polynomial. The accuracy of the approximation depends on the order of the expansion, which, in turn, determines the minimum number of lattices required, as detailed in §6. Systems exhibiting complex behaviours require higher-order polynomials to approximate them to within the predefined tolerance, and consequently require more lattices. In the presence of restrictions (8.1) and (8.2) on the geometric and material parameters, the system will often need to be expanded to an order higher than the target polynomial so that the system response is within the predefined tolerance and the parameter restrictions are satisfied. We describe such a system as over-expanded. Example 2 ( §8d) is an exactly expanded system, Example 3 ( §8e) is under-expanded. Examples 1 and 4 ( §8c,f, respectively) are over-expanded.
To highlight the errors associated with the framework presented in this paper, as opposed to the approximation method, we use Lagrange interpolation to generate the target polynomials. The nature of the numerical framework we have presented allows for discrete implementation of the approximation constraints; the target polynomial may be replaced by discrete target points, assuming the behaviour at locations other than the target points is not of interest. See, [32,33] for a comprehensive discussion on interpolation and approximation techniques of smooth functions by polynomials.
When it is not the energy per se but rather its first or second derivative (i.e. reaction force or stiffness, respectively) which are prescribed, the highest prescribed derivative is chosen to be matched where integration constants are inherently determined from the chosen material coefficients (cf. (3.8)).
We are now ready for the examples. Each of these represents only a single point in a solution space whose richness allows the designer to explore and optimize designs in accordance with requirements. The non-uniqueness of the solution space maximizes both the range of applicability of the algorithm and the physical plausibility within the design space itself.
(c) Example 1: monostable system with prescribed nonlinear stiffness   (kN mm -1   ) target stiffness, E≤ system actual stiffness, X ≤ system Taylor stiffness, E P ≤ system actual force, X ¢ system Taylor force, E P ¢ system actual energy, X -X 0 system Taylor energy, E P -E P 0 reaction force (MN) to a tolerance of ±0.2 kN. The system is exactly expanded using two lattices, the resulting second-order bi-stable target reaction force, actual reaction force, polynomial reaction force and the absolute error in the approximation are shown in figure 5a. The system's corresponding fourth-order bi-stable energy profiles are shown in figure 5b.
As discussed in §7c, the even roots of the stiffness polynomial are stable, i.e. two equilibria exist for the quadratic reaction force, one stable and one unstable, as shown in figure 5. It should be noted that due to the relatively benign nonlinearity in the system only two lattices with an expansion of O(H 4 ) exactly matching the order of the target polynomial is required to capture the (e) Example 3: multi-stable system with specified local minima Extension of the requirements to design a bi-stable system to be multi-stable is trivial within the current framework. In this example, two internal and one boundary stable equilibria are selected in accordance with equation ( to a tolerance of ±0.05 kN, with stable equilibria specified atH ∈ {0, 0.32, 0.64}. The system is composed of four lattices, the resulting fifth-order multi-stable target reaction force, actual reaction force, polynomial reaction force and the absolute error in the approximation are shown in figure 6a. The system's corresponding multi-stable energy profiles are shown in figure 6b. As discussed in §7c, even roots of the stiffness polynomial are stable. The depth of the energy wells associated with the stable equilibria are particularly shallow in this example, hence the resulting transition force between the energy wells is also rather small. Consequently, a tighter tolerance has been enforced to ensure all equilibria of interest are captured. Within the current framework, it is trivial to adjust the constraints in this manner. The tighter tolerance coupled with the predefined bounds on the material parameters has required an over-expansion of O(H 8 ) to match the target reaction force O(H 6 ) to within the specified tolerance.
The geometric and material parameters of this system are presented in the electronic supplementary material, table S3. The effect of insufficient over-expansion can be observed at the extension limit,H = 1, and highlights the underlying assumption of the algorithm. Although the polynomial expansion can be tailored to match the target function to within a given tolerance, if the order of expansion is not sufficient to capture the actual nonlinear response over the entire domain the two responses will diverge; the polynomial response may remain within the tolerance band while the actual response does not. This behaviour is exacerbated with decreasing number of lattices (and thus expansion order) as well as over restrictive material parameter bounds, resulting in the point of divergence moving closer to the origin (expansion point), a common feature of low-order Taylor expansions. This example presents a multi-stable system where the transition force required to snap between stable regimes as well as their extensional position have been specified. Limit points, i.e. a local maximum in the reaction force, are defined by specifying both a root in the stiffness and the desired magnitude on the reaction force. Care must be taken to ensure this point is a local maximum.
In addition to prescribing the magnitude and location of limit points, three stable equilibria have (optionally) been specified in between them, one below the first snap-through load,H = 0, one (anywhere) between the two snap-through loads and one at the extension limitH = 1, such that the system's preferred stable configuration changes when each snap-through load is exceeded. Specifying monotonically increasing snap-through loads and their associated stable equilibria in this manner allows loading thresholds experienced by the system to be preserved upon unloading; the elastic system behaves as a reusable sensor.
The sixth-order target reaction force polynomial is defined as the Lagrange polynomial of the data points with the snap-through loads specified as 5 ± 1 kN and 20 ± 1 kN atH = 0.16 andH = 0.696, respectively. Accordingly, an over-expanded system of six lattices with an expansion of O(H 11 ) has been employed to capture this highly nonlinear behaviour over the entire domain to within a tolerance of ±1 kN in the reaction force, as shown in figure 7a. The system's corresponding multistable energy profiles are shown in figure 7b. Under the restrictions (8.1) and (8.2), a significantly higher level of expansion is needed to control the behaviour of the system to within the given tolerance over the entire domain, mitigating the effects observed in figure 6a. This demonstrates how increasing the nonlinearity through addition of limit points and/or raising the order of target polynomial by a single degree may require more than the addition of a single lattice so that the system is sufficiently over-expanded. Equally, a reduction in the number of lattices may be achieved through relaxation of the bounds on the material parameters. The geometric and material parameters of this system are presented in the electronic supplementary material,

Conclusion
Fully exploiting the paradigm shift from the linear-only to linear-and-nonlinear structural responses remains an open challenge. To address this challenge, and highlight the new opportunities it presents, we have developed a systematic approach to designing bespoke, elastically nonlinear, one-dimensional extensional responses. We achieve this by harnessing the behaviour of the helical lattice structure proposed by Pirrera et al. [18] and demonstrating how a hierarchical system comprising several lattices coupled together offers a route to obtain bespoke elasticity via a robust design process. We show how an increase in the geometric complexity of the system increases the richness in response by identifying how each unique lattice (i.e. one that possesses a unique geometric parameter ψ which is the scaling factor between local and global extension) added to the system extends the possible response types. This permits the formulation of a basis that spans polynomial energies to arbitrary accuracy in L ∞ . We use this framework to demonstrate how several desirable responses, such as multi-stability and snap-through buckling, can be achieved and explore the requirements for their existence. The design space for the helical lattice system proposed here is extremely robust. The energetic description of the system permits an exploration of the interplay between pre-strain, geometry and stiffness. The non-dimensional approach allows behaviours to be considered at various geometric length scales, characteristic material stiffnesses and system curvatures. This leads to a description that is applicable for both structural (from civil to MEMS length scales) and meta-material-like requirements.
Data accessibility. The datasets supporting this article have been uploaded as electronic supplementary material. Authors' contributions. A.P., I.V.C. and M.P.O. conceived the research idea and derived the model for the system