Negative extensibility metamaterials: phase diagram calculation

Negative extensibility metamaterials are able to contract against the line of increasing external tension. A bistable unit cell exhibits several nonlinear mechanical behaviors including the negative extensibility response. Here, an exact form of the total mechanical potential is used based on engineering strain measure. The mechanical response is a function of the system parameters that specify unit cell dimensions and member stiffnesses. A phase diagram is calculated, which maps the response to regions in the diagram using the system parameters as the coordinate axes. Boundary lines pinpoint the onset of a particular mechanical response. Contour lines allow various material properties to be fine-tuned. Analogous to thermodynamic phase diagrams, there exist singular “triple points” which simultaneously satisfy conditions for three response types. The discussion ends with a brief statement about how thermodynamic phase diagrams differ from the phase diagram in this paper.


Introduction
Mechanical metamaterials exhibit mechanical properties not observed in most natural and artificial materials. Properties are the result of an engineered internal structure rather than the chemical composition of the bulk material. They feature non-natural elastic constants including negativity in the elastic modulus, shear modulus, bulk modulus or Poisson's ratio [1][2][3][4]. Due to advances in rapid prototyping it is possible to fabricate 2D and 3D lattice structures composed of simple microstructural spring and bar elements [5][6][7][8]. The field of mechanical metamaterials is concerned with understanding how an unconventional macroscopic response is influenced by the deliberate structuring of elastic elements.
The mechanical property under investigation is negative extensibility. Negative extensibility occurs when a material contracts in length so as to oppose the line of increasing external tension. The concept is illustrated using a one dimensional elastic bar in Fig. 1a. In the initial state the bar expands smoothly in tension. At a critical destabilizing load the bar undergoes a marked contraction corresponding to the forward More accurately, it would possess a negative "incremental" stiffness since the slope of the stress-strain curve shifts from positive to negative after an initial strain. Negative incremental stiffness has been observed experimentally in pre-strained materials [9] and cellular structures under displacement control [10,11] and realized theoretically through a nonlinear force potential between particle constituents [12]. Elasticity theory does allow for negative incremental stiffness, in particular when the object is constrained [13][14][15].
A reoccurring motif in metamaterials research is that properties are tunable [16][17][18]. Certain classes of mechanical metamaterials display multiple stable states of equilibrium [19][20][21]. Tuning capabilities are achieved during the process of state switching, which often occurs by reconfiguration of the lattice or unit cell geometry [22,23]. In the previous paper in the series, Karpov et al. [24] developed a geometrically nonlinear unit cell structure that exhibits negative extensibility for the external degree of freedom. The unit cell is shown in Fig. 1b. Contraction in tension is made possible because the unit cell is bistable. Moreover, several different types of mechanical responses are observed depending on the choice of member stiffnesses and unit cell dimensions. Each mechanical response is associated with a distinct pattern of state switching over a load-unload cycle. Of particular interest is the negative extensibility response, which is characterized by a "pinched" hysteresis loop. During the contraction, the direction of the displacement opposes the applied tensile force. The unit cell does negative work. Energy is released from the system to the surroundings. Thermodynamically, this is a solid-to-solid condensation reaction. When integrating F ·du over a load-unload cycle there is an interval where the displacement during unloading lies above the displacement during loading. The area enclosed by this region represents negative work. Overall, the net work of the cycle remains positive, which is consistent with other dissipative materials such as shape memory alloys [25,26]. However, since the hysteresis loop is pinched there is a small, but significant region where the area under the force-displacement curve is negative.
A systematic analysis of the strain energy function is presented. The analysis is quasistatic. The previous work [24] invoked a simplifying assumption regarding the strain energy stored in the elastic bars. In this paper, no approximations for the strain energy function are used. A phase diagram is cal-culated which maps the mechanical response of the unit cell to regions in the diagram. The axes of the phase diagram are system parameters that define member stiffnesses. Boundary lines mark the onset of a particular mechanical response. The critical boundaries that divide monostability from bistability are determined using a set of constraints, which is informed by the principles of catastrophe theory [27][28][29]. Analogous to thermodynamic phase diagrams, there exist singular "triple points" where the system simultaneously satisfies conditions for three different mechanical responses. Contour lines in the diagram specify the width of the hysteresis loop and the intensity of the strain change over the forward transformation.
The response can be fine-tuned to control various material properties such as the change in length during the contraction or the critical load at the onset of the contraction. A step-bystep example is carried out to demonstrate how to use the phase diagram in order to design a unit cell with specific material properties.

Strain energy function
The unit cell in Fig. 1b is composed of five linear elastic elements: a middle bar of stiffness k 1 , top and bottom bars of stiffness k 2 and springs on either side of stiffness k 3 . The initial dimensions are uniquely determined by a horizontal length L, vertical height H and skewness h. The skewness is the degree of offset from the horizontal. At h = 0 the structure is rectangular. The structure is pinned at the midpoint of the center bar. Due to symmetry there exists only two independent degrees of freedom, the displacements u and v. Rollers constrain their movement along the vertical direction. The displacement u is the external degree of freedom since the force F is applied at this node. Mechanical work is being done over this node. The displacement v is the internal degree of freedom. No work is being done at this node.

Balance of energy
The stored strain energies in the middle bar π 1 , top and bottom bars π 2 and the spring π 3 are linear elastic relations written as The potential or strain energy function is defined as the strain energies of the elastic elements minus the energy of the external load Upon substitution of the strain energies for each element, the fully dimensional potential for the unit cell is The potential is placed into dimensionless form so as to reduce the number of independent system parameters. This is accomplished by dividing the entire equation by an energy term associated with the side springs, k 3 H 2 . The dimensionless potential U is defined Quantities are redefined in terms of dimensionless parameters. There are four system or design parameters that specify unit cell dimensions and member stiffnesses Parameters a and b represent dimensionless stiffnesses scaled by geometry terms. In the previous work [24], the definitions of a and b were the same. Parameter s is the dimensionless skewness. Parameter r is the aspect ratio. In this paper, the aspect ratio is defined as the height H divided by the length L, which is the inverse of the standard definition of the aspect ratio. The reason for adopting this definition is because the flat structure with infinitesimally small height represents the limiting case when the aspect ratio approaches zero, r → 0. Two independent state parameters are defined, x and y, which are the dimensionless displacements for the two degree-offreedom unit cell The control parameter f is the dimensionless load The final form of the dimensionless potential is Although the material laws are all linear elastic, the potential is highly nonlinear as a consequence of geometry. Because of this nonlinearity, for certain combinations of element stiffnesses and dimensions it is possible for the solution to bifurcate or split such that two stable states of equilibrium exist for a single loading condition. The result is hysteretic bistability. In the subsequent analysis it will be possible to pinpoint what combinations of system parameters lead to bistability.

Force-response curves
The unit cell structure undergoes several qualitatively distinct types of force-response curves depending on the values of the system parameters a, b, r and s. The external degree of freedom x is of primary interest and is the point where the load is applied. The strain of the structure is defined There are five mechanical responses possible: monostable (MS), superelastic (SE), superplastic (SP), negative extensibility of the superelastic type (NESE) and negative extensibility of the superplastic type (NESP). The force-response curves corresponding to each type are shown in Fig. 2. Prior to generating a force-response curve, the four system parameters must be specified. Under quasistatic loading conditions, equilibrium is satisfied when the total potential energy of the system is a minimum [30]. Equilibrium conditions require U x = U y = 0. These two equations are sufficient to solve for the two unknown displacements, x and y, at discrete values of the load f . To simulate the load-unload cycle, the load is incremented gradually and decremented gradually using a small enough step size f . The Newton-Raphson method solves for the equilibrium displacements at each step. The converged displacements of the previous step are used as the trial solutions for the next step. The first characteristic response of the unit cell is monostability (MS), which is drawn as the purple curve in Fig. 2. The response is a smooth, continuous expansion and contraction of the elastic members during load-unload cycles.
The second characteristic response is superelasticity (SE), which is shown as the green curve in Fig. 2. In general, superelasticity is the ability of a material to accommodate large applied strains and recover to its original configuration when the load is removed. Superelasticity is characterized by the cyclic dissipation of energy, a process called hysteresis [31]. Initially, the unit cell deforms smoothly with increasing tension until it reaches a critical destabilizing load. At this point the structure experiences a discontinuous expansion or snap-through action in both degrees of freedom until it reequilibrates reaching a new stable configuration. The forward superelastic transformation A → B is comparable to a solidto-solid vaporization reaction since energy is absorbed by the system. Work is positive since the displacement is in the direction of the applied load. During load removal, the reverse transformation B → A occurs as an abrupt contraction but at a lower critical force than the forward transformation. The reverse superelastic reaction is analogous to solid-to-solid condensation since energy is released from the system to the surroundings. Work is negative for the reaction. During cycles of load-unload, the net work of the cycle is positive. Although a net amount of energy is absorbed by the system, the absorbed energy would be dissipated as heat. The total energy dissipated in the process is proportional to the area enclosed by the hysteresis loop [32]. Shape memory alloys show qualitatively a similar superelastic hysteresis effect in their stress-strain curve due to the stress-induced martensitic transformation [33].
The third characteristic response is superplasticity (SP), which is shown as the gray curve in Fig. 2. Superplasticity occurs when a material accommodates large applied strains and remains in a highly deformed state after load removal. Superplastic materials are hysteretic. Like the SE response, the SP response exhibits a critical point of destabilization associated with the forward phase transformation. However, the fundamental difference between SE and SP is the location of the critical force for the reverse transformation. For the superplastic response when the load is removed, the unit cell remains deformed in state B. In order to restore the structure to its original configuration, the load must be applied in the opposite (negative) direction. Only then will the structure overcome an energy barrier, which allows it to snap back to its original state. The superplastic response of the unit cell parallels that of a shape memory alloy. For instance, if a shape memory alloy is stressed beyond a certain threshold or the temperature is too low then the larger volume martensitic phase persists and the metal remains in a superplastic state after load removal [34]. Compressive stresses or heat is used to restore the superplastic alloy back to its original position.
The last two characteristic mechanical responses of the unit cell are the negative extensibility responses. They feature pinched hysteresis loops. They are defined by the negative extensibility transition, which is the marked contraction in the external degree of freedom during the forward phase transformation. Like the other bistable responses, during the forward transformation the structure snaps-through as the middle bar rotates past the horizontal configuration. However, due to the relative stiffnesses of members in the unit cell, a large expansion of the internal degree of freedom actually facilitates the contraction of the external degree of freedom when equilibrium is restored in state B. In other words, the contraction is made thermodynamically possible through an expansion of the internal degree of freedom. During unloading, if the structure fully recovers to its original configuration then the response is termed negative extensibility of the superelastic (NESE) type, which is shown as the red curve in Fig. 2. If the structure remains in a deformed state requiring compressive forces to restore it to the original configuration then the response is defined as negative extensibility of the superplastic (NESP) type, which is shown as the blue curve in Fig. 2. The two negative extensibility responses reveal a secondary hysteresis loop when the load is increased beyond the point of contraction. The secondary loop is associated with the fourfold switching pattern, In Fig. 3, a dynamic simulation is performed in order to determine whether or not the NESE response is maintained when the system is excited with a time-dependent load. The simulation uses the second-order Verlet numerical integration scheme [35]. During the transient process of state switching, the snap-through action causes the system to oscillate and then decay as vibrational energy is released into the surroundings. Compared to the quasistatic case, the dynamic case displays a similar effect magnitude of the contraction. In general, if the load grows slowly enough with time and the damping is large enough then the intensity of the contraction in the dynamic case should be about the same as for the quasistatic case. In order to preserve the NESE effect at higher load rates, a practical material may consist of two phases including a bistable core lattice and a viscoelastic matrix to dissipate kinetic energy.

Mapping response to phase diagram
The phase diagram maps the mechanical response of the unit cell to values of the dimensionless system parameters. Each point on the diagram corresponds to a specific unit cell design -i.e., a unit cell with specific dimensions and member stiffnesses. The aspect ratio r and skewness s are fixed for a given diagram. Fig. 4 is the phase diagram for r → 0 and

Phase diagram calculation
A complete phase diagram consists of five boundary lines and two sets of contour lines.

Boundary lines
The boundary line Γ E divides superelasticity (SE or NESP) from superplasticity (SP or NESP). Conditions for Γ E define a point of destabilization associated with the reverse transformation B → A with the added constraint that the critical force occurs at zero load, f c = 0. The structure is neither definitively superplastic nor superelastic. Conditions are written as The functions g 1 = U x , g 2 = U y and g 3 = det H are set to zero. They take the five arguments in the parentheses. Aspect ratio and skewness are not shown explicitly in the system of equations since they are fixed for a given phase diagram. The determinant of the Hessian is defined When the determinant of the Hessian is zero and equilibrium conditions are satisfied in the vicinity of the critical point, these are the criteria for a point which is simultaneously an inflection point and a stationary point -i.e., a saddle point. Physically, this is a point of structural destabilization. Points of destabilization mark the onset of a vertical jump discontinuity in the hysteresis loop. Γ E is the special case when the point of destabilization occurs at f c = 0, which fixes the reverse transformation to occur at zero load. This added constraint on the critical force reduces the number of unknowns to four: a, b, x and y. Using one of the unknowns as a running variable reduces the number of unknowns to three. The three equations are now sufficient to generate Γ E provided the initial guess to start the Newton-Raphson scheme leads to physical solutions for the design parameters a and b. Points along Γ S are considered cusp points. The conditions are based on two sets of destabilization criteria (6 total equations). Each set has different variables for the displacements x 1 , y 1 and x 2 , y 2 . However, the critical force is the same for both sets, f 1 = f 2 . When a hysteresis loop initially forms, critical points of destabilization corresponding to the forward and the reverse transformations are nearly touching (the sides of the loop are nearly touching). Only at a cusp singularity will the two points of destabilization meet exactly at a single critical force. The hysteresis loop is reduced to a point. The conditions are written as 1 (a, b, x 1 , y 1 , f 1 For the Γ S conditions to work a single running variable is taken, usually a or b, reducing the number of unknowns to six: x 1 , x 2 , y 1 , y 2 , f 1 and either a or b. The system is fully specified. Figure 5 (left) shows a force-response curve generated at a converged Γ S point. A hysteresis loop is observed once inside the region bounded by Γ S . The boundary lines Γ M , Γ N and Γ O all use the same set of conditions. The trial solution in the Newton-Raphson algorithm determines what curve will be accessed. Although the conditions are the same, the nature of the bifurcation associated with each boundary line is different. For Γ M the conditions specify the onset of a pinched hysteresis loop. This happens when a point of destabilization defined by g 1 , g 2 and g 3 occurs at the same critical load and critical displacement as a stable solution of equilibrium defined by g 1 and g 2 . Therefore, f 1 = f 2 and x 1 = x 2 where x 1 , y 1 and  1 (a, b, x 1 , y 1 The force-response curve generated at Γ M is depicted in Fig. 5 (middle, top). At a constant critical force during the forward transformation A → B, the system destabilizes and then re-equilibrates with no net change in the displacement, Although not shown, there is still a large expansion in the internal degree of freedom y, y 2 = y 1 . Fig. 5 (middle, bottom) shows that when crossing over Γ M the response becomes NESE. In general, the boundary Γ M divides the usual SE or SP hysteretic response from the more interesting NESE or NESP pinched hysteretic response. The boundary line Γ O , referred to as the coalescence curve, divides the regions of NESE and NESP from a region of "apparent monostability." Apparent monostability refers to the fact that in the region above Γ O a second mathematical solution of equilibrium exists but it is not realized under normal loading conditions. The bistable response is inacces-sible. Figure 5 (right) shows that when crossing over Γ O the two hysteresis loops merge. The response shifts abruptly to monostability. The reason the above conditions work for Γ O is related to the merging of the two hysteresis loops. The initial point of destabilization at the onset of the negative extensibility transition A → B is stabilized when the two loops merge, x 1 (unstable) = x 2 (stable) at Γ O . Coalescence points along Γ O are singularities where a solution of equilibrium, which satisfies g 1 and g 2 , becomes a point of destabilization such that it now satisfies g 3 in addition to g 1 and g 2 . This higher-order bifurcation causes two hysteresis loops to emerge in the force-response curve.
The boundary line Γ N , referred to as the nucleation curve, divides monostability from apparent monostability. The logic in explaining why Γ N uses the above set of conditions is identical to the logic in explaining why Γ O uses it. A monostable point, which has one solution of equilibrium satisfying only g 1 and g 2 , shifts to a destabilization point such that now satisfies g 3 in addition to g 1 and g 2 . Mathematically, a second stable solution bifurcates from the destabilization point although it cannot be accessed under physical loading. Hence, the system is apparently monostable. To put another way, an inaccessible hysteresis loop develops when crossing over Γ N .
The effect of Γ N is more readily observed on the stability diagram shown in Fig. 6 (right). The stability diagram pertains to the vertical slice of the phase diagram at b = 5.21. It plots the system parameter a as a function of the critical destabilizing forces associated with a phase transformation. The peak of the top black curve corresponds exactly to a point of Γ N . As a is decreased, meaning the stiffness of the top and bottom bars is decreased relative to the side springs, physical bistability becomes accessible at Γ O . Below Γ O there are two hysteresis loops with a total of four points of destabilization. This particular stability diagram does not reveal a point of Γ M . This is because Γ M is different from the Γ O and Γ N bifurcations that separate monostability from bistability. Instead, Γ M alters the nature of the already existing bistable response causing it to become pinched.
Finally, the stability diagram reveals a strong cusp associated with Γ S at a = 0.04. This cusp point corresponds to the formation of the secondary hysteresis loop, which occurs along the dotted line Γ S in the phase diagram. Just above (Γ S ) there are two hysteresis loops and four-fold switching. Just below it there is only one loop and two-fold switching. At approximately a = 0.007 for a large critical force, f c = 21, there is another point of hysteresis loop formation associated with the secondary dotted Γ S . The four-fold switching pattern is again realized when a is decreased below this point.

Contour lines
Two sets of contour lines are plotted on a complete phase diagram: the force ratio φ and strain intensity I SE contours. The force ratio φ gives the relative width of the hysteresis loop. It is defined as the critical force for the reverse transformation, f B→A = f 2 , divided by the critical force for the forward transformation, Figure 7 (left) shows the definition of φ in relation to a force-response curve. The critical force at the forward transformation f 1 is the denominator because it will always occur at loads greater than zero. In contrast, the critical force at the reverse transformation f 2 shifts to zero at Γ E and then negative for the superplastic (SP or NESP) responses. Therefore, positive φ implies an SE or NESE response while negative φ implies an SP or NESP response. As the force ratio increases in the positive direction, the difference between the critical force at the forward transformation and reverse transformation grows less and less and the hysteresis loop shrinks in width. On the other hand, decreasing the force ratio widens the gap between the two transformations. The most negative φ corresponding to the most superplastic responses have the widest relative hysteresis loops. The conditions used to generate the force ratio contours are related to Γ S except that the value of f 2 is specified relative to f 1 based on a pre-defined value of φ. The conditions fix the width of the hysteresis loop based on the ratio of forces at the forward and reverse transformations. This is accomplished by first rearranging Eq. (16) to solve for f 2 The result is substituted into the set of destabilization criteria associated with the reverse transformation, Eqs. (18d)-(18f). Like Γ S , the force ratio conditions are composed of two sets of destabilization criteria The system parameters a and b are the phase diagram axes. They are frequently used as running variables in the Newton-Raphson method. The value of φ is specified at the beginning of the program in order to fix the value of f 2 relative to f 1 . The six equations are now sufficient to solve for the six unknowns: x 1 , x 2 , y 1 , y 2 , f 1 and either a or b. The force ratio conditions are related to the boundary lines Γ S and Γ E . When φ = 1 this implies f 2 = f 1 reducing the conditions to Γ S . When φ = 0 this implies the critical force for the reverse transformation f 2 = 0 reducing the force ratio conditions to the three independent equations associated with Γ E . The strain intensity I SE is a measure of the strain change over the forward transformation A → B. It is defined as the strain reached after the transformation ε 2 minus the critical strain at the onset of the transformation ε 1 normalized by the critical strain ε 1 Figure 7 (right) shows the definition of I SE on a forceresponse curve. Since the response is NESE, the forward transformation is a contraction. In the figure, the change in strain over the transformation ε c = ε 2 − ε 1 is negative, which gives a negative I SE . In general, a larger magnitude I SE in either the positive or negative direction means the expansion or contraction is relatively more intense -i.e., it is associated with a sooner onset (lower critical strain ε 1 ) and a larger change in strain. The most intense contractions are realized for low a and high b in the NESE or NESP regions. At the intersection of Γ E and Γ O is where most negative I SE occurs while still being NESE. This intersection point is referred to as the triple point, P E O . Figure 5 (right, top) shows qualitatively the force-response curve at the triple point. The which is equivalent to The equilibrium displacement after the transformation x 2 is substituted into the equilibrium criteria, Eqs. (22d) and (22e). The strain intensity conditions are written as The value of I SE is fixed at the beginning of the program. As an added constraint it reduces the number of unknowns by fixing the value of x 2 relative to x 1 . In the Newton-Raphson method, usually a or b is used as a running variable. The five equations are now sufficient to solve for the five unknowns: x 1 , y 1 , y 2 , f 1 and either a or b. As a caveat, the subscripts 1 and 2 in the I SE conditions and Γ M conditions pertain to points at the beginning and end of the forward transformation. In contrast, the subscripts 1 and 2 in the φ conditions and Γ S conditions pertain to the locations of the forward and reverse transformation.

Complete phase diagrams
Four complete phase diagrams are generated in Fig. 8a through d. The unit cell for each phase diagram is defined by its aspect ratio and skewness. The force ratio φ contours are drawn as dashed lines varying from black (wide hysteresis loop) to green (narrow hysteresis loop). The strain intensity I SE contours are drawn as solid lines that vary from cyan (contraction) to blue (small expansion) to red (large expansion). Increasing aspect ratio vertically shifts the NESE region to higher a while increasing skewness horizontally distorts the NESE region so it expands over a greater range of b. The unit cell for Fig. 8d was chosen for its symmetry. At s = 0.5 the middle bar is the same length as the top and bottom bars.

Design example
To illustrate how to use the phase diagram for the design of a unit cell with specific material properties the following problem statement is given: Design a unit-cell with aspect ratio r = 0.25 and skewness s = 0.25 that exhibits a NESE response with a contraction of 1 mm, strain intensity I SE = −0.01, force ratio φ = 0.2 and critical load at the forward transformation F 1 = 1000 N. After specifying r , s, I SE and φ the designer is allowed to specify two additional dimensional quantities. These quantities may give information about hysteresis loop, unit cell dimensions or member stiffnesses. This is because the dimensional potential, Eq. (5), has 10 variables whereas the dimensionless potential, Eq. (10), has 8 dimensionless parameters. This leaves two free dimensional quantities to be specified. In this case, the length of the contraction and the value of the critical load at the onset of the contraction are chosen. Specifying r and s determines the phase diagram to be used. Next, specifying I SE and φ fixes the dimensionless system parameters a and b. On the phase diagram for r = 0.25, s = 0.25, Fig. 8c, the intersection of I SE = −0.01 and φ = 0.2 occurs at b = 8.7842, a = 0.07851. Since the four system parameters are known the force-response curve can be generated. The following dimensionless quantities are obtained from the force-response curve: f 1 = 1.2371, ε 1 = 2.0219 and ε 2 = 2.0017. These quantities are not free variables. They are properties of the hysteresis loop which are determined based on the dimensionless system parameters.
For this particular problem, there are 9 dimensional quantities to be defined Of these 9 quantities, the stiffnesses k 1 , k 2 and k 3 and dimensions H , h and L are the 6 quantities needed to actually construct the unit cell. The phase diagram is essential because it relates properties of the mechanical response such as the critical force F 1 and the displacements before and after the contraction, u 1 and u 2 , to the design of the unit cell. Based on the statements in the problem, the following 9 variables are known: r = 0.25, s = 0.25, a = 0.07851, b = 8.7842, f 1 = 1.2371, F 1 = 1000 N, u = 1 mm, ε 1 = 2.0219, ε 2 = 2.0017. The definitions of the above variables are the 9 independent equations used to solve for the 9 unknown dimensional system parameters (24d) The system is fully determined. After performing algebraic substitution to solve for the unknowns, the requirements are satisfied with the following design:

Effect of changing aspect ratio and skewness
The four phase diagrams in Fig. 8 show roughly how the regions and contours change with varying aspect ratio and skewness. In this section, the maximum intensity of the contraction in the NESE region, which occurs at the triple point P E O , will be used as a metric to compare unit cells of different aspect ratio and skewness. Furthermore, diagrams will be developed to understand the distortion of the NESE region with increasing aspect ratio and skewness.

Intensity of the contraction at the triple point
The triple point P E O is the intersection of Γ E and Γ O . At this point the structure simultaneously satisfies criteria for The force-response curve at the triple point appears as Γ O like in Fig. 5 (right, top) where the two hysteresis loops merge.
The Γ E conditions fix the location of the reverse transformation at zero load, which is shown qualitatively in the same figure as the black dashed line. For an individual phase diagram with specified r and s, the a and b at the triple point P E O correspond to a unit cell design with the most intense contraction while still being a NESE response. Therefore, the following definition is introduced The definition of I N E SE is the same as I SE , Eq. (19). The subscript implies that only the NESE response is considered. More intense contractions are possible in the NESP region, especially when approaching Γ O with decreasing a. The maximum intensity of the contraction for the NESE region is used as a metric for comparing unit cells of different aspect ratio and skewness. Table 1 shows the effect of varying aspect ratio at constant s = 0. Decreasing the aspect ratio increases the maximum intensity of the contraction at the triple point. An infinitesimally small aspect ratio, r → 0, leads to a structure with small height compared to its length.  This flat structure represents the theoretical limit by which to compare other unit cell arrangements. Figures 9 and 10 visualize the effect of changing both the aspect ratio and skewness on the intensity of the contraction at the triple point. The plane curves in Fig. 9 are cross-sections of the three dimensional r -s-I N E SE surface in Fig. 10. The global minimum corresponds to rectangular unit cells, s = 0, with infinitesimally small aspect ratio, r → 0.

Evolution of the NESE region
The NESE region is defined by the three boundary lines Γ E , Γ M and Γ O . The region takes the shape of a distorted triangle. The intersection points P E O , P E M and P M O are the The points are plotted as aspect ratio is varied from r → 0 to r = 0.5 with a step size r = 0.05. Four complete NESE regions are drawn. For a linear increase in aspect ratio there is an exponential vertical shift in the location of the NESE region. A thinning and elongation of the NESE region is also associated with increasing aspect ratio corners of triangle. Like the triple point, P E O , the other two intersection points are computed precisely using combined conditions for the corresponding boundary lines. This process requires careful choice of the trial solutions for each of the unknowns. However, once a solution is found for an intersection point there is less trial and error if the same type of intersection point is desired for different aspect ratio and skewness. An iterative program was constructed to gradually increment r and s and compute each of the intersection points at each step. Using this technique, it is possible to track the evolution of the NESE region with changing r and s. Figure 11 shows the drift of the NESE region with increasing aspect ratio at constant skewness. A larger aspect ratio shifts the NESE region vertically. Figure 12 shows the drift of the NESE region as a function of increasing skewness at constant aspect ratio. Skewed structures see a horizontal distortion of the NESE region. The NESE region doubles in area when skewness is increased from s = 0 to s = 0.25 at constant r = 0.25.

Conclusion
The unit cell response is governed by a non-convex, nonharmonic and consequently highly nonlinear potential, Eq. (5). Having determined the sets of equations needed for each Gamma line calculation gives insight into systematic numerical analysis of concave nonlinear potentials. Gamma lines that mark the onset of bistability as well as Gamma lines that reveal the nature of the mechanical response are computed precisely using combined conditions of equilibrium and destabilization. Plotting the Gamma lines and contour lines on a phase diagram makes it possible to design the unit cell to achieve a specific response. In essence, the phase diagram is a visual representation of the mechanical response as a function of design parameters that specify unit cell dimensions and elastic member stiffnesses. However, in order to calculate the Gamma lines, contour lines and intersection points requires solving systems of nonlinear equations. In this paper, the Newton-Raphson method is used. The major downside to this method is that convergence to physically significant values depends on the initial guess for the unknowns. A poor choice of trial solutions for the unknowns may lead to divergence or non-physical roots. The upside to the Newton-Raphson method is that if the functions are smooth and continuous and the initial guess is close enough to the roots then the method will converge.
A phase diagram differs from a normal graph in which one variable is plotted as a function of another. In a phase diagram, the coordinate axes all represent independent variables and the coordinate space shows the state of the system at equilibrium [36]. While it is interesting to draw comparisons between the phase diagram in this paper and microstructural phase diagrams used in materials science, the two diagrams are fundamentally different. For instance, the Gibbs phase rule does not apply to the diagram presented here. Regions do not correspond to microstructural phases in thermal equilibrium. Instead, a "phase" is a distinct mechanical response, one of either MS, SE, SP, NESE or NESP. Furthermore, the design parameters a, b, r and s are not chemical components that constitute a thermodynamic phase. Nevertheless, it may be instructive to think of them as independent constituents that when specified lead to a particular mechanical response. In the future, the phase diagram methodology will be extended to other systems whose response is governed by a nonlinear potential in hopes of discovering new metamaterial phenomena.