Shape Entropy of a Reconfigurable Ising Surface

Techniques arising from papercraft – origami and kirigami – inspire efforts to impart three-dimensional shapes to two-dimensional materials [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. Origami and kirigami both work by imposing one-dimensional modifications (folds or cuts) into surfaces. Here we investigate an alternative method of shape control that inserts zero-dimensional objects – isolated disclinations – into an atomically thin sheet. Inversion of bistable conical disclinations provides a means to control the shape of the surface, with tunable topgraphy down to the length-scale of protein secondary/tertiary structure and energy barriers against shape change that are tunable from below room temperature to values much larger. A tunable topography in a biologically relevant size regime suggests that reconfigurable shape-driven binding interactions may be possible with such sheets.

A disclination removes or adds a wedge of material in the lattice of a two-dimensional sheet, yielding a conical or hyperbolic local region; in a hexagonal lattice the least costly such defects are 5 and 7-membered rings. Since a cone (in a continuum approximation) has a local $C_{\infty}$ symmetry for rotations about its apex, its mechanical states are anticipated to be particularly simple: just “up” and “down”, related by a local inversion. The reduced rotational symmetry of the saddle suggests more complex behavior, as it may be able to assume more than one in-plane orientation depending on the local mechanical environment. Here we investigate a sheet containing multiple disclinations balanced between positive and negative (i.e. asymptotically flat) and determine the multiplicity and character of the shape metastability thereby instilled in the 2D layer. Anticipating our conclusions, the sheet’s conformational freedom is dominated by Ising-like bistable cone degrees of freedom; hence we call this means of shape control Isigami.

To create such a surface, we begin with a so-called Haeckelite [14, 15] structure and “inflate” each ring with a penumbra of hexagons, similar to inflation of larger fullerenes from $\text{C}_{\text{60}}$ [16, 17]. We choose a Haeckelite in which the pentagonal rings form a slightly deformed kagome lattice (Fig. 1) [18]. The shape entropy phenomena that we seek to study should not require an ordered lattice of disclinations; this choice is for computational convenience, and to explore whether ideal mechanical frustration may bring the set of possible shapes closer to mutual degeneracy. Other disclination patterns are also possible, and those formed by phase field crystal methods have been examined independently for their effect on sheet toughness in work by Zhang, Li and Gao [19]. We choose graphene as an archetype, since its mechanical response is well-studied and methods exist to functionalize the surface to control mechanical stiffness, interfacial energetics, and hydrophilicity [20, 21, 22]. Since the behavior described below is largely geometrical in origin, conclusions derived from graphene Isigami should generalize to other atomically thin two-dimensional materials.

We construct both periodic and finite-patch regions of material hosting twelve cones and twelve saddles on a kagome lattice as depicted in Fig. 1 and model the mechanical response using the AIREBO potential [23] as implemented in LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [24], as it is well-validated for the mechanical response of nanoscale sp² carbon [25, 26, 27, 28, 29]. We first examine a finite sheet whose disclinations are separated by six hexagons ($\sim$20 Å); the edges are terminated by hydrogen. This separation is large enough that even cones on the edge of the sheet have two well-defined metastable “up” and “down” orientations. Shape variation is most dramatic in a finite patch of material, as the periodic boundary constraint is absent. The patch is systematically forced into every possible set of cone orientations by mechanically inverting various sets of cones. To this end, a force of $\pm 1.7$ eV/Å is applied to the five atoms at the apices of various cones according to the desired up/down states and the system is then allowed to evolve for 500 femtoseconds while the atoms are constrained to move no more than $0.1$Å each timestep. A constraint force is also applied evenly across every atom in the system to zero-out the net force. The structure is then relaxed either through a thermal bath at 300 K for 50 ps followed by a linear ramp anneal over 30 ps from 300 K to 3 K, or (for the periodic systems described below) through a shorter 1 ps thermal bath followed by a series of conjugate gradient minimizations where the unit cell geometry is allowed to relax during alternate optimization steps. These methods both produce well-converged energy minima (converged to less than 100 meV in systems with thousands of atoms), but the anneal method performs better for the finite patch in avoiding local minima.

For every initial choice of the twelve up/down cone configurations, the system relaxes into a distinct metastable shape, i.e. all $2^{12}$ nominal Isigami configurations are accessible. Fig. 2 shows some examples of these shapes, labeled by the cone configuration: the different shapes that are metastably assumed by the same sheet can be dramatically different. Several of these configurations, including the five highest in energy, were simulated at 800 K in an NPT ensemble for one nanosecond, with the sheet maintaining its original Ising state. This high-temperature “challenge” anneal was also applied to a subset of configurations that were repeated to create a $2\times 2$ supercell patch, providing assurance that shape metastability is preserved for larger patch sizes as well.

Collapsing the reflection symmetry about the sheet midplane, we actually obtain $2^{11}$ distinct Isigami configurations for this 12-pentagon patch, with an additional approximate reflection symmetry about the short diagonal (broken very weakly by the detailed termination of the patch). Additional symmetries in the periodic case further reduce the number of symmetry-distinct configurations per unit area, but the conformational entropy of metastable Isigami states remains extensive.

As noted earlier, whereas cones have two distinct conformations (i.e. up and down), the metastable orientations of the principal axes of saddle-point disclinations are not naturally two-fold; their number depends on their surroundings. We do in fact find certain conformations that have the same set of cone orientations (i.e. the same Isigami state) but visibly different overall shapes, sometimes differing by over an eV in energy and stable in shape up to at least 800 K over half a nanosecond. Fig. 3 depicts two such states. While we have not exhaustively quantified the frequency of these non-Ising states, our observations suggest they contribute a number of configurations equivalent to roughly one additional Ising degree of freedom for a 12-cone sheet; i.e. their contribution to the overall configurational entropy is modest. Although these non-Ising states do sometimes correspond to different orientations of the saddle disclinations, we do not find a simple rule governing their appearance, as reflected by the examples shown in Figure S5. Similar non-Ising states are seen in the larger $2\times 2$ supercell patch. Simulations of this larger patch (and larger periodic systems) reveal a similarly pervasive metastability to Ising conformations, suggesting that the system has an extensive conformational shape entropy.

How are these states distributed in energy? The “density of shapes” (i.e. the number of shapes in a given energy range) for the finite 12-cone patch is strikingly close to Gaussian, as shown in Fig. 4. It is tempting to consider this outcome as an expression of the central limit theorem for multiple nearly independent energetic contributions from different disclinations, but the stronger deviations from a Gaussian shape seen in a similarly-sized periodic patch argue against a simple application of this notion. The distribution for the finite patch may be more Gaussian because its boundary condition produces multiple symmetry-distinct pentagonal disclinations within the patch and these heterogeneous local environments, under the action of elastic interactions that span the entire patch, yield something closer to the requirements of the central limit theorem. In contrast, every pentagonal disclination under periodic boundary conditions is symmetry-equivalent. For less well-ordered disclination networks we anticipate the periodic case to more closely approach the finite-patch distribution.

Mechanical intuition suggests that nearby cones will have pair-wise antiferromagnetic interactions, since opposing up/down orientations produce compatible sidewall slopes. In practice, we find that longer-range and many-body elastic interactions between disclinations are too strong to admit a strict frustration-derived degeneracy. While these non-idealities (and also non-Isigami states) are sufficiently important that the sheet is not a simple nearest-neighbor antiferromagnet, the energies of different sheet configurations do reveal an overall nearest-neighbor antiferromagnetic trend, as the energy varies systematically with the number of oppositely-oriented nearest-neighbor cones for both the periodic-boundary and finite-patch cases (Fig. 5). The effective $J$ of this trend is proportionately smaller for a sheet with 1-hexagon separation between the 5-fold rings (Fig. S1). As noted in Supplementary Information, it was not possible to describe the sheet energetics with a simple cluster expansion [30]; this was likely due in part to the long-range nature of elastic interactions amongst nearby cones (i.e. when even a single cone inverts, we see significant changes in shape across the entire patch, as shown in Figure S4) and in part to the presence of a modest number of additional states not fully specifiable by an Isigami state vector, as noted earlier.

These reconfigurable surfaces provide a means to access a multitude of distinct sheet topographies from a single patch of multiply disclinated 2D material: an obvious area of potential application for such an object is shape recognition and binding of biomolecules. Ideally, such a surface would have topography on the length scale of protein tertiary structure, a hydrophilic surface that does not denature biomolecules, and a tunable barrier against shape change so that both dynamic ($\sim k_{B}T$) and static ($\gg k_{B}T$) regimes can be accessed. The $\sim 0.5$ nm inter-pentagon spacing of the sheet depicted in Fig. 6 provides topography on a biologically relevant length-scale, as demonstrated by the comparison to insulin, a small protein, in the same figure. The use of graphene oxide [31, 32] rather than hydrophobic bare graphene [33, 34, 35] could introduce the needed hydrophilicity, and perhaps also increased flexibility [36].

The question of barrier height control requires more detailed study. The inversion barriers of isolated graphene cones have been simulated previously [37]; the barrier depends on the chemical character of the conical apex but is generally multiple electron volts for well-separated cones. Energy barriers against cone inversion will depend on the material of the sheet, the separation between conical disclinations, functionalization [37], or grain boundaries interconnecting disclinations [38, 39]. We observe in the periodic sheets that the stability of the Ising and non-Ising states indeed varies with cone separation, and also that the barriers against inter-conversion of the non-Ising states are typically much lower than those between Ising configurations. We monitored the energy during cone inversion for a few periodic sheets with different cone separations; although the distribution of inversion barriers could not be exhaustively quantified, the barriers all remained on the order of electron volts, likely lowered somewhat (compared to isolated cones) by irregularities in the local environment. The typical barrier height scales roughly linearly with the separation between pentagonal disclinations, being around 1 eV for a 1-hexagon separation and around 6 eV for the 6-hexagon separation. At the 1-hexagon separation a significant fraction of the Isigami states destabilize (see Figure S2) while almost all of the non-Ising states disappear, where our criterion for stability is surviving the room-temperature anneal mentioned earlier.

Since most of the cases considered above have inversion barriers much greater than $k_{B}T$, we now consider a means to “thermalize” these barriers even at larger disclination separations: oxidizing the sheet and etching away the central five-membered ring and one or more surrounding rings to reduce the mechanical deformation associated with nucleating inversion, terminating the resulting apical multivacancy with hydrogen atoms. To this end, we oxidized a six-hexagon-separation periodic sheet by random placement of oxygen atoms equally onto upper and lower surfaces followed by a structural relaxation and subsequent removal of any unreacted oxygen, yielding a functionalized sheet with 12.6 at% oxygen (Fig. S3). The apices of this sheet up to three surrounding rings were etched away by manual deletion, and the resulting edges terminated with hydrogen. These simulations used a version of the REAXFF reactive force field [40, 41, 42] designed [43] to describe carbon, hydrogen and oxygen [44] and updated to more precisely describe chemical properties of C/H/O/N systems [45, 46, 47, 48]. Neither the structural relaxation following etching (nor that which followed oxygenation) caused cone inversion, but during a subsequent MD simulation at 800K we observe rapid ($\sim 1$ ns) spontaneous inversion of cones; we expect a similar simulation without periodic boundary conditions to be at least as floppy. Although these simulations were performed in the gas phase, we anticipate that the qualitative conclusion of low, thermally accessible inversion barriers will be robust against the introduction of aqueous conditions (which may themselves further reduce the barrier). In constrast, the same system with un-etched apices shows no inversion over 1 ns at 800 K.

Shape reconfigurability at the nanometer-to-micron scale is a signature of biology, but one that is not straightforwardly expressed in non-biological systems. We describe a strategy towards this end that uses the bistability of disclinations to afford a large library of shapes from a single 2D sheet. Selection from this library could either happen spontaneously – for example, through binding interactions with targets under conditions of shape lability – or be imposed through mechanical deformation by a specified indenter. We demonstrate the second of these through large-scale molecular dynamics simulations of indented Isigami sheets. The simulated sheets are large enough to be compatible with experimental nano/micromechanical manipulation [49, 50] (i.e. hundreds of nanometers across) and make use of simple radial and annular indenter geometries, as depicted in Supplementary Videos S1-S3. The three sheets simulated all start from the same asymptotically flat Ising configuration. The indenters’ pressures force all the cones underneath to orient similarly, which creates linear structures of unbalanced (i.e ferromagnetically aligned) Ising state. The antiferromagnetic interaction between neighboring cones in the indented region then drives a sheet distortion that splays these cones away from each other to minimize their unfavorable co-alignment: the resulting crumples “bunch up” material under the indented region. For the radial indenter this bunching reduces the circumference at fixed radius from the center of the indenter, while the annular indenter does the opposite; thus they create large-scale Gaussian curvatures of opposite sign. The cross-shaped identers of different widths also produce different bend angles. We apply four-fold symmetric indenters to show that vertices defined by the intersection of four folds can be easily created, while a generic crumpled fold would involve only three. These deformations should be reversible by compressing between flat planes and re-indenting.

As this approach is purely geometrical, it should apply to any atomically thin sheet, which could be polar (hBN), hydrophilic (graphene oxide [31, 32]) or stiffer than graphene [51]. These phenomena should also occur in more irregular disclination networks that could be produced by e.g. growing a 2D layer on a roughened substrate, as strict degeneracies due to mechanical frustration on a regular lattice were neither observed nor required: mechanical bistability is a local property of each conical disclination, one whose antiferromagnetic near-neighbor interactions then favor a large shape multiplicity. The degree of metastability is tunable from highly stable to potentially room-temperature reconfigurable. Possible longer-term applications for 2D shape reconfiguration could involve diverse shape libraries exposed to binding targets under shape-labile conditions, selected for binding, and then utilized under shape-conserving conditions. On longer length scales, linear arrays of disclinations could recover phenomena akin to more traditional fold-based origami wherein mean curvature plays a more prominent role. The presence of disclinations (with associated geometrical phases [52]) and local strains (with associated pseudomagnetic fields [53, 54, 55, 56]) may also introduce unusual electronic and plasmonic [57] properties.