Self-limiting states of polar misfits: frustrated assembly of warped-jigsaw particles

Abstract

We study the ground state thermodynamics of a model class of geometrically frustrated assemblies, known as warped-jigsaw particles. While it is known that frustration in soft matter assemblies has the ability to propagate up to mesoscopic, multi-particle size scales, notably through the selection of the self-limiting domain, little is understood about how the symmetry of shape-misfit at the particle scale influences emergent morphologies at the mesoscale. Here we show that polarity in the shape-misfit of warped-jigsaw puzzles manifests at a larger scale in the morphology and thermodynamics of the ground-state assembly of self-limiting domains. We use a combination of continuum theory and discrete particle simulations to show that the polar misfit gives rise to two mesoscopically distinct polar, self-limiting ribbon domains. Thermodynamic selection between the two ribbon morphologies is controlled by a combination of the binding anisotropy along distinct neighbor directions and the orientation of polar shape-misfit. These predictions are valuable as design features for ongoing efforts to program self-limiting assemblies through the synthesis of intentionally frustrated particles, further suggesting a generic classification of frustrated assembly behavior in terms of the relative symmetries of shape-misfit and the underlying long-range inter-particle order it frustrates.

---

1 Introduction

Geometrically frustrated assembly (GFA) has recently emerged as an organizing principle for complex, self-organized structures in soft matter. In geometrically frustrated systems, such as in the frustrated ordering of magnetic spins,^1–3 a locally preferred arrangement of constituents is not compatible with the constraints on uniformly filling space. In soft matter assemblies, this takes the form of molecular or particulate subunits whose shape and interactions favor a local “misfit” with the longer-range inter-particle order. This basic paradigm applies to a broad range of soft matter systems, including chiral and polar phases of liquid crystals,^4–6 chiral membranes^7–12 and protein bundles,^13–18 colloidal crystals on curved surfaces,^19–22 and ill-fitting colloidal and nanoparticle assemblies.^23–31 The key effect of frustration is to introduce intra-assembly gradients of misfit that extend to multi-subunit scales, leading to elastic costs that accumulate with increasing domain sizes. This accumulation of elastic energy³² with domain size notably results in self-limitation, in which the equilibrium assembly states have well-defined, finite, multi-subunit sizes, a behavior that is not possible for nearly all other classes of equilibrium assembly.^33,34

These exotic and potentially useful behaviors of GFAs, in combination with an increasing array of synthetic techniques for fabricating geometrically-programmable assembly units,^{25,26,35–37} have inspired theoretical efforts that attempt to connect microscopic features of the subunits and their interactions to the emergent larger-scale morphologies.^{27,29,30,38,39} Such studies, which range from continuum theories to numerical simulations of discrete misfitting subunit models, address basic and still broadly open questions: given a “misfitting” subunit design, under what conditions does self-limited assembly occur, and what morphologies emerge from it? From the point of view of engineering self-limiting states, it is particularly important to understand what design features of “misfit”, such as particle shape and interactions, control assembly thermodynamics at the largest possible scales and ultimately, to translate those features to synthetically accessible modifications of subunits.

It is in this context that we consider a basic classification of distinct GFA mechanisms in terms of their misfit symmetry, shown schematically in Fig. 1 for the case of two-dimensional assemblies.† This notion applies to the cases where the degree of frustration may be continuously increased, through some experimental parameter, from the frustration-free cases where there is some type of long-range order such as in 2D hexagonal or square crystalline packing. The introduction of frustration takes the form a local shape distortion incompatible with the (unfrustrated) background order. The symmetry of this local distortion, and in particular its mismatch with the underlying frustration-free order, is what we refer to as the misfit symmetry. For example, for the case of crystalline assemblies on spherical surfaces,^19–22 shown in Fig. 1a, the underlying (frustration-free) order corresponds to a hexagonal lattice which uniformly tiles the 2D plane. Imposing positive Gaussian curvature on the assembly, say by forcing particles to assemble on a spherical drop, disrupts the uniform hexagonal lattice order, since equilateral triangles do not (in general) tile spherical surfaces. Since there is no particular directionality of this distortion of hexagonal order due the spherical curvature, we may refer to this as scalar misfit. This is distinct from the other examples of directional misfit shown in Fig. 1, in which the misfit is characterized by distortions along one or more specific axes in the 2D subspace of the assembly. For example, in the model of trapezoidal-shaped particles shown in Fig. 1b, frustration arises from angling opposite edges of a square tile away from each other, which generates a directionality or polarity between the long and short edges. In the model of wedge-shape particles studied in ref. 30, shown in Fig. 1c, frustration arises from angling opposite faces away from each other such that the shape and interactions favor a 2D square lattice assembly that is frustrated by a simultaneous preference for opposite curvatures in two orthogonal directions. Hence, the preference for negatively-curved, minimal surface shapes of an otherwise crystalline membrane can be characterized as a quadrupolar or tensorial misfit. Finally, in the model of irregular hexagon assembly developed in ref. 23 and 24, shown in Fig. 1d, frustration is introduced by shortening and lengthening alternate edges of an otherwise hexagonal elastic tile. We can characterize this deformation of an unfrustated 6-fold symmetric particle into a lower 3-fold symmetry particle as hexapolar misfit.

Fig. 1 Classification of geometrically frustrated assemblies in terms of the symmetry or directionality of the frustration due to shape distortion and its mismatch with the underlying frustration free inter-particle order. The top row shows unfrustrated assemblies and their preferred arrangement while the middle row shows the corresponding frustrated assemblies and their distortions (position, orientation, and shape) from those preferred arrangements. The bottom row shows example realizations of self-limiting assemblies. Scalar frustration (a) is represented by the crystallization on spherical surfaces as in ref. 19. Vectorial frustration (b) can be represented by assemblies of warped particles with a directed shape asymmetry.²⁹ Quadrupolar frustration (c) is represented by double wedge particles that prefer locally nonzero Guassian curvature from ref. 30. Hexapolar frustration (d) is represented by warped hexagonal particles with long and short edges from ref. 24.

A central question, which we address in this article for a specific model of polar frustration in so-called warped-jigsaw particles,²⁹ is how does the misfit symmetry of the building blocks impact the morphologies of self-limiting states that are selected by their assembly thermodynamics? Previous theoretical and computational studies of a diverse range of GFA models found that 2D minimal energy morphologies of self-limiting states tend to break the symmetry of the underlying background order.^{18,23,24,29,30,40} For example, while low energy shapes of hexagonal crystallites favor generally compact and six-fold symmetric shapes that optimize edge energies, the growth of elastic energy with the domain size of crystals on a sphere favors strip-like morphologies^19,41 that grow arbitrarily long in one direction but remain narrow in the other (see example self-limiting morphologies in Fig. 1, bottom). The basic mechanism that promotes domain anisotropy, sometimes referred to as “filamentation”, stems from the elasticity of frustrated domains, which requires inter-domain gradients in the elastic stress. From an energetic point of view, it is favorable to orient the strain gradients along the narrow direction of the domain, so that the range of strains (i.e. difference between maximal and minimal strains) remains small. Hence, such scalar mechanisms of GFA, like crystalline assemblies on spherical surfaces, favor strip-like morphologies that break the underlying six-fold symmetry of the unfrustrated assembly, leading to multiple degenerate states (e.g. degenerate strip morphologies associated with the three low-energy growth directions). A similar phenomenon occurs in spin lattice models of frustrated assemblies with scalar frustration⁴⁰ (which also exhibit similar branching morphologies as crystalline assembly on spherical surfaces), frustrated hyperbolic crystalline assemblies,³⁰ and irregular elastic hexagon assemblies,^23,24 in which optimal self-limiting morphologies break symmetry by remaining narrow in one direction of their 2D assembly and growing long in the other. Notably, for these examples (e.g. scalar, quadrupolar and hexapolar), the misfit symmetries are commensurate with broken symmetries of self-limiting morphologies. Hence, for these symmetries the possible resulting “ribbon-like” morphologies remain degenerate, and the physical assembly processes must spontaneously break symmetry to form self-limiting states.

In this article, we investigate the geometrically-frustrated assembly and self-limiting morphologies of a polar misfit model of planar warped jigsaw (WJ) particles, and show that the polarity or directionality of local misfit profoundly reshapes the energetic landscape of optimal morphologies. Specifically, we show that, unlike the scalar, quadrupolar, or hexapolar cases, the effect of misfit polarity on the elastic energy landscapes generically lifts the degeneracy between multiple competing strip-like morphologies. In general, we find that these non-degenerate morphologies are distinguished by the type of mechanical strain-gradients they develop along their narrow dimensions. We first consider the case where polar misfit is aligned with a lattice (binding) direction and find that the distinguishable mesoscopic morphologies have distinct patterns of local gradients, reflected in the polarity of the large scale assembly of either straight-polar or circular finite-thickness ribbons, as shown in Fig. 1b. We show that these deformation patterns are characterized by distinct size dependencies of the elastic costs of frustration, resulting in the morphologies undergoing “frustration escape” to bulk, unlimited states via distinct thermodynamic pathways. We predict a phase diagram of ground-state morphologies that determines the stable regions of both straight-polar and circular finite-thickness ribbons, as well as a uniformly strained bulk (unlimited) state.

We next consider the case where misfit polarity is misaligned with the underlying two-dimensional lattice (binding) directions, which we show leads to a complex evolution of finite-width strip morphologies that mix features of the straight-polar and circular ribbon states. Last, we show that while the elastic cost of frustration is generically smallest for ribbons growing along the misfit polarity direction (i.e. straight-polar ribbons), ground-state ribbons generically still favor growth along the binding directions due to the energy cost of having free boundaries.

The article is organized as follows. In Section 2, we introduce a microscopic model of WJ particle assembly and the ingredients of a discrete subunit computational model of the elastic ground-states. We then describe the pairwise inter-particle deformation modes and formulate a continuum elastic description of assembled domains of WJ particles. Then in Section 3, we study the case of misfit polarity aligned to one of the 2D binding directions, showing the emergence of two distinct, self-limiting morphologies, characterized by distinct patterns of strain gradients. We construct a phase diagram of the ground-state morphologies in terms of the cohesion between bound particles. In Section 4, we study when the misfit polarity is not aligned to either lattice direction. We first consider domains whose free edges are aligned along binding directions before generalizing to free edges that align with off-lattice directions to determine the preferred growth direction of long ribbons. We conclude with a discussion of the implications on the broader emergent effects of frustration on assembly, as well as for the engineering of self-limiting assemblies through design and fabrication of intentionally misfitting subunits.

2 Warped-jigsaw particle model

Here, we introduce the warped-jigsaw (WJ) particle model of Spivack and coworkers.²⁹ The shape of WJ particles is described by rigid quadrilaterals in the 2D plane, shown in Fig. 2. The geometry and interactions of the WJ particles simultaneously favor binding along rows and columns in a 2D square lattice as well as gradients in particle orientation along a certain in-plane direction, which frustrates the uniform lattice order. Each edge has a pair of elastic binding sites with specific interactions that bind together particles along two orthogonal directions: “left-right” (dark red-dark orange and light red-light orange sites) and “top-bottom” (dark blue-dark green and light blue-light green sites) edges of neighboring particles. Such specific interactions, of the type that can be generated e.g. via sequence-specific DNA linker binding,^42,43 imply that the left (bottom) edges only bind favorably to the right (top) ones and that bound neighbors favor a certain relative orientation. The relative orientations are described in a local frame of orthonormal unit directors {ê₁, ê₂} aligned to the respective horizontal and vertical directions of a particle, but can rotate in the xy-plane within an assembly. In the unfrustrated state, the particle is simply a square (Fig. 2a) and when bound to neighbors, the ground state arrangement is that of a square lattice (Fig. 2b), with the frames of all particles aligned uniformly.

Fig. 2 Discrete subunit model for warped-jigsaw particles. (a) Unfrustrated square subunit. Colored circles indicate elastic binding sites with specific binding between dark green and dark blue, light green and light blue, dark yellow and dark orange, and finally light yellow and light orange. The interaction between binding sites is described by ideal zero rest-length springs. The pairs of binding sites on each edge are separated by a distance v while the overall subunit size is a. (b) Unfrustrated square subunits bind to form a square lattice. (c) and (d) Frustrated quadrilateral subunit and the preferred binding horizontal and vertical binding of neighbors. The unit vector N̂ (red arrows) characterizes the direction of preferred curvatures κ₀i.e. the rotation of ê₁ and ê₂ along directions perpendicular to N̂, described by the unit vector T̂. (e) Deformation of an assembly (light red) away from its uniform 2D square lattice reference state (light blue), resulting in a displacement field u(x, y) and rotation field θ(x, y).

Frustration is introduced by shape asymmetry, as shown in Fig. 2c, where opposite edges of the unfrustrated square shape are skewed away from each other, forming taper angles α₁ and α₂ along the horizontal ê₁ and vertical ê₂ directions, respectively. This shape asymmetry due to the taper angles results in a preferred relative rotation between neighboring particles and consequently, nonzero preferred curvatures along columns and rows of these particles (Fig. 2d), a local arrangement that is incompatible with the positional order (square lattice) of the unfrustrated state. As outlined in the introduction, this frustration is polar or vectorial in nature in that it has a direction and magnitude. We characterize the polarity of shape-misfit in terms of the preferred curvature vector

κ₀ = κ₀N̂, (1)

where κ₀ is the magnitude of preferred curvature and N̂ characterizes the preferred bending direction of rows or columns of particles. The resulting direction of rotation of the particle frames along a row or column can be described by the unit vector T̂ = −ẑ × N̂, where ẑ is normal to the xy plane of assembly. The preferred curvatures along rows and columns are given by κ₁ ≡ κ₀·ê₁ = (2 tan α₂)/a and κ₂ ≡ κ₀·ê₂ = (2 tan α₁)/a, which describe the amount of preferred bending orthogonal to local ê₁ and ê₂ directions, respectively. Hence, quantifies the magnitude of the misfit while N̂ = −sin ψê₁ − cos ψê₂ describes its orientation, which we parameterize by the angle ψ in the local particle frame with sin ψ = κ₁/κ₀ and cos ψ = κ₂/κ₀.

Shown in Fig. 2e is an example of an assembly formed by trapezoidal WJ particles with ψ = 0, whose shapes are symmetric along their ê₁ (horizontal) directions and asymmetric along their ê₂ (vertical) directions, as characterized by the downward misfit polarity N̂ = −ê₂. This asymmetry favors curvature along only the “horizontal” rows causing the assembly to deviate from an unfrustrated square lattice reference state. To describe and model the elastic energy WJ assemblies, we take two approaches. In the first, we perform energy minimizations of domains of elastically bound rigid WJ particles. As we show, emergent ground states take the form of highly anisotropic rectangular domains with either finite numbers of columns M or rows N. We first assume that the domain boundaries align with binding directions so that domain edges are characterized by broken bonds along either the ê₁ or ê₂ directions, although we relax this assumption and consider the possibility of domain edges aligning along different directions in Section 4. Cohesive costs of an exposed (unbound) WJ particle edge in direction ê_α is parameterized by the binding energy Σ_α.‡ Deformations of interactions between bound particles (i.e. 2 per bound edge) are modeled by Hookean elastic springs with zero rest-length between specific binding sites, with stiffnesses parameterized by k_α for bonds in the ê_α direction. For a given M × N array of WJ particles with a specific frustration and cohesive stiffness, the elastic ground state is simply computed by energy minimization of the coupled network of springs and rigid particles (see Appendix A).

A second approach to model the elastic ground states is to study the continuum elastic theory of arbitrary rectangular domains of assembled WJ particles. This kind of continuum theory should hold in the limit of large assemblies (i.e. N and M ≫ 1) along with an assumption that the deformations of the ground state from uniform 2D square lattice order, which we take to be the reference state of our theory, are sufficiently small. In other words, we describe the positions and orientations of WJ particles relative to their corresponding configuration in an unfrustrated lattice assembly, where they have positions r = xx̂ + yŷ and particle frames ê₁ = x̂, ê₂ = ŷ. Upon introducing frustration, particles displace and re-orient as they relax the elastic costs of frustration. As illustrated in Fig. 2e, we describe the displacement of WJ centers by the field u(r) and their rotations relative to the fixed 2D Cartesian coordinates by the angle θ(r), i.e. the particle frames can be written as ê₁ = cos θ(r)x̂ − sin θ(r)ŷ. As we describe in the following section, our continuum elastic model is based on a small-θ expansion of the elastic energy of discrete WJ particle arrays (see Appendix B).

It should be noted in this work that we neglect excluded volume interactions in both the simulations and the continuum elastic theory, the latter of which is of course a must in order to make the theory analytically tractable. As can be seen in Fig. 2d and e, there is some degree of overlap due primarily to the shape of the corners of the WJ particles, a feature that will likely be present in experimental realizations of the WJ particles. However, as observed in the work of Spivack and coworkers,²⁹ which studied a discrete WJ particle model with excluded volume interactions, there is good agreement between the continuum elastic theory and simulations despite the inclusion of excluded volume interactions. In fact, in the limit of small frustration where the taper angles α₁, α₂ tend to zero, the degree of potential overlaps can be significantly reduced as the strains also tend to zero. It is also worth noting that the “shape” of the WJ particles that is responsible for generating polar frustration is more due to the positions of the elastic attractive sites on the edges and less the full body of the particles, meaning that one could redesign the particles (e.g. by truncating the corners) to reduce overlaps and still retain the same polar frustration. With this in mind, we neglect excluded volume and focus on the elastic effects of polarized frustration in a frustrated assembly.

2.1 Elasticity of intra-domain distortions

We briefly introduce the mechanical modes of intra-domain deformation, relate the corresponding moduli to microscopic features of the frustrated WJ particles, and finally formulate the continuum elastic energy of finite WJ domains.

As shown schematically in Fig. 3, deformations of inter-particle order relative to the preferred non-zero orientation of the WJ particles can be parameterized by three types of pure elastic modes: inter-particle stretching, shearing, and bending. The elastic cost of these modes is controlled by the geometry and stiffness of elastic bonds between bound particles. At the pairwise level, we can describe the elastic energy of a pair of bound edges along the ê_α direction as

where δ_‖, δ_⊥ and δθ parameterize the respective stretch, shear, and bend displacements (labeled in Fig. 3, left). The elasticity of these respective modes are determined by stiffnesses and placement of the inter-WJ bonds (see Appendix B) and are given bywhere v ≤ a is the distance between linkers along a particle edge (see Fig. 2a and c) and k_α is the elastic bond stiffness along the ê_α direction. Note that in all simulations from here on, we set the distance between the linkers to its maximum value v = a. This choice is not so important as it simply sets a scale for the bending moduli B_α, but it does allow us to maximize the bending moduli. As an example of the complex interplay between these deformation modes in frustrated WJ domains, we show in Fig. 3 the simulated elastic ground state of a M × M square domain, illustrating the spatial distributions of the three elastic deformation modes for a case ψ = 0 or polarity along the “vertical” ê₂ direction. We return to a detailed analysis of these patterns below, but upon initial inspection it is clear that shearing is largest along left-right free edges, compression/stretching is largest along the bottom-top free edges, and (un)bending of rows away from their preferred curvature is largest in the core of the domains.

Fig. 3 Three types of “pure” elastic deformations between pairs of WJ particles illustrated relative to their preferred non-zero orientation of parallel edges where there is no extension of the ideal springs between binding sites. Extension or compression occurs when bound edges of subunits are pulled apart or pushed together by a displacement δ_‖ perpendicular to their edges. Shear occurs when the bound edges of subunits slide against each other by a displacement δ_⊥ parallel to their edges. (Un)Bending occurs when bound subunits rotate away from their preferred orientation of aligned edges by a total angle δθ. Corresponding examples of 10 × 10 structures colored by the amount of extension/compression, shear, and bend are shown on the right for isotropic elastic constants k₁ = k₂.

To understand these emergent patterns of intra-domain distortion, their dependence on domain shape, and their impact on assembly thermodynamics, we study the ground states of a continuum elastic energy of WJ particle assemblies. This continuum theory was also derived in ref. 29, but only studied for a limited set of “tall” finite-width ribbon domains with isotropic binding constants. As shown in Appendix B, the elasticity of WJ domains in the continuum limit is described by a combination of positional and orientational strains. Positional strains, resulting from inter-particle stretching and shearing, are measured by the 2D tensor

γ_ij(r) = ∂_iu_j − [R_ij(−θ) − δ_ij], (4)

where the matrix R_ij(θ) describes 2D counterclockwise rotations by an angle θ and δ_ij is the Kronecker delta. The form of the strain is invariant under simultaneous rigid rotations of particle centers and orientations, ensuring global rotational invariance of the elastic energy. Orientational strains are measured by gradients of particle rotation relative to the preferred row bending

β(r) = ∇θ − κ₀T̂, (5)

where T̂ describes the preferred direction of 2D inter-particle rotations in the local particle frame. Note that the three pure elastic modes illustrated in Fig. 3 correspond to γ₁₂ = β₁ = 0 for extension and compression, γ₁₁ = β₁ = 0 for shear, and γ₁₁ = γ₁₂ = 0 for bend. In the analysis below, we linearize the theory for small frame rotations away from the uniform reference lattice aligned with x̂, ŷ, resulting in

γ_ij(r) ≃ ∂_iu_j − θ(r)ε_ij; β(r) ≃ ∇θ − κ₀T̂₀, (6)

where ε_ij is the 2D anti-symmetric Levi-Cevita tensor and T̂₀ = −cos ψx̂ + sin ψŷ is the (uniform) approximation of the preferred rotation rate in the small θ(r) limit.

The elastic energy is computed as an area integral over a rectangular domain (the undistorted reference state)

where in the small-θ approximation we neglect the difference between {ê₁, ê₂} and {x̂, ŷ}, thereby identifying the corresponding coordinate indices as (α = 1) → (i = x) and (α = 2) → (i = y), and use the fact that the shear and stretch moduli are equal for the WJ model, that is, Y_α‖ = Y_α⊥. For a given domain, the elastic energy ground state satisfies a set of Euler–Lagrange equations corresponding to force balance

Y_x∂_x²u_x + Y_y∂_y²u_x = −Y_y∂_yθ; Y_x∂_x²u_y + Y_y∂_y²u_y = Y_x∂_xθ, (8)

and torque balance

B_x∂_x²θ + B_y∂_y²θ = −Y_x(∂_xu_y − θ) + Y_y(∂_yu_x + θ). (9)

Note that the coupling between displacements and rotations in mechanical equilibrium reveals the mechanism of frustrated elasticity. Preferred nonzero row curvature, which is favored by the taper angles of the WJ particles, results in nonzero rotation gradients ∇θ ≠ 0 and consequently generates a source field for positional strain gradients in eqn (8). At the same time, rotations of the local WJ lattice directions away from the particle frames generates a source for orientational strain gradients in eqn (9). Boundary conditions at the free edge of finite domain ∂D require vanishing forcesand torqueswhere n̂ = (n_x, n_y) describes the normal to domain boundary.

The equilibrium equations (eqn (8) and (9)) and boundary conditions (eqn (10) and (11)) constitute a set of coupled, linear PDEs that can be solved by a series expansion (see Appendix C). Non-trivial and spatially-inhomogeneous solutions are required by the torque free boundary conditions in eqn (11) which produce nonzero θ gradients (i.e. row bending) at free edges where n̂·T̂ ≠ 0 and vanishing θ gradients where n̂·T̂ = 0 (i.e. no column bending).

2.2 Thermodynamics of finite-size domains

We consider a model of the assembly thermodynamics of WJ domains that is determined by the interplay between the elastic costs of frustration and the cohesive gains of binding new particles to the free boundaries of finite domains. Such an approximation obviously neglects the effects of other possible low-symmetry configurations (e.g. defective) as well as entropic effects that are present at finite temperature. Nevertheless, the purpose is to capture the predominant morphologies that emerge in regimes where nearly all subunits belong to assembled domains, i.e. super-saturated conditions. In this case, assembly favors morphologies that minimize the per-subunit free energy of assembly, which includes the (free-)energetic gains and costs of forming multiunit structures.³³ In our discrete WJ assemblies, the assembly free energy per subunit for a M × N domain takes the formThe first term is the constant cohesive gain for bound particles within the bulk, parameterized by attractive binding per edge −Σ_x and −Σ_y along respective horizontal rows and vertical columns. The second and third terms together account for the per-subunit cost of unbound free edges on the vertical and horizontal sides of the domain. The final term denotes the excess energy and represents the elastic costs of frustration per subunit, computed via the minimization of the spring network energy in our discrete model. We can similarly use our continuum theory to compute the elastic energy, in this case for a domain of horizontal and vertical dimensions, W = Ma and H = Na respectively. In this case, the excess energy cost is determined by the mechanical equilibrium of eqn (7) for the W × H domain,where the * denotes the solution of the Euler–Lagrange equations.

Inspection of eqn (12) reveals the basic physical mechanism of the size-dependent thermodynamics of frustrated WJ domains. On one hand, the cost of open boundaries thermodynamically favors growing all dimensions of an assembly to arbitrarily large sizes, in order to minimize the fraction of unbound edges at the free edges of the domains. On the other hand, the elastic cost of a frustrated WJ domain grows superextensively with domain size,³⁴ meaning that it grows faster than the increase in subunit number in the domain. This implies that in certain regimes, the competition between boundary and frustration costs selects at least one of the two assembly dimensions to be finite. In the following sections we first analyze the thermodynamic assembly landscape resulting from this interplay for the particular case of ψ = 0 polarization before then considering a broader range of misfit polarities.

3 Alignment of misfit polarization and lattice anisotropy: ψ = 0 case

Here, we analyze the case of ψ = 0 where misfit polarity is aligned to a 2D assembly direction such that κ₀ is parallel to ê₂ and the “horizontal” rows of the WJ domain favor curvature. We first study the mechanics and elastic energy of inhomogeneous ground states as a function of domain shape and then describe the thermodynamics optimal morphologies.

3.1 Elastic energy landscape: ψ = 0

The detailed solution of the continuum theory (eqn (8) and (9)) is presented in Appendix C. For this particular polarization, the free boundary conditions (eqn (10) and (11)) are such that the curvature of WJ rows achieves its preferred magnitude κ₀ on the vertically oriented sides of the domains, while its value in the interior of the domain may deviate from κ₀ depending on the domain size. The relaxation derives from a competition between preferred row bending and accumulated costs of either intra-row stretching or inter-row shears, and is thus characterized by “bend penetration” length scales determined by the ratios between elastic constants,Note that for our particular model design and in terms of the spring stiffnesses (eqn (3)), these elastic length scales become and λ_y ∼ v. In Appendix C.2, we introduce an approximation based on observations made on the structures (e.g.Fig. 4d) that allows us to determine the elastic ground state energy of W × H domains, which is given in closed form bywhere is the (per particle) flattening cost to completely unbend the rows from their preferred curvature, and the rescaled width and height arerespectively. Notably, it can be easily verified that the second term in eqn (15) vanishes in the double limit of w → ∞ and h → ∞, which shows that the bulk state is uniformly flattened, except at the free boundaries. In general, the excess energy approaches this upper limit with increasing size. It is interesting to note that in the limit of asymptotically small assemblies w ≪ 1 and h ≪ 1, the excess energy takes on the form (see Appendix C.2)which is symmetric about h = w, indicating that elastic anisotropy only enters beyond quadratic order in domain sizes.

Fig. 4 Energetics and shape of WJ particle assemblies. (a) Excess energy landscapes for M × N assemblies from simulations and continuum theory. Simulations have isotropic elastic constants with k_x = k_y. (b) Excess energy of tall assemblies with fixed height N = 30 and varying width M (left) and wide assemblies with fixed width M = 30 and varying height N (right). The inset compares the excess energy of tall and wide assemblies. (c) Centerline curvatures as a function of position s along an arc for tall M × 30 (left) and wide 30 × N (right) assemblies. The inset shows the height dependence of the curvature. (d) Example assemblies highlighting the degree of extension/compression, shear, and bend in tall and wide assemblies.

Fig. 4 shows a comparison of excess energy from simulations of M × N assemblies of discrete subunits with continuum theory of W × H (or Ma × Na) domains, showing qualitative and quantitative agreement notwithstanding the small-θ approximation made in the continuum formulation. The excess energy landscapes show two zero energy states corresponding to single rows (wide) and columns (tall). As the width of the tall assemblies and the height of the wide assemblies increases, there is an accumulation of strain due to the frustration or shape misfits, leading to a super-extensive growth in excess energy. The excess energy for tall (fixed height and varying width) and wide (fixed width and varying height) assemblies are shown in Fig. 4b, left and right, respectively. In general, these show that the excess energies increase monotonically with both the width and height of the domains, albeit with obvious differences in their rate of increase. In particular, the excess energy of tall assemblies accumulates more gradually than that of the wide assemblies.

In Fig. 4c and d, we show the mechanical equilibria of the two classes of anisotropic morphologies, tall (H ≫ W) and wide (W ≫ H) ribbon-like domains, from discrete particle simulations comparing their distinct intra-domain mechanics. First, we note that away from the boundaries of either type of domain, elastic ground states relax to a pattern that is effectively uniform along the longer dimension. Next, we observe that tall and wide morphologies both continuously approach the flattened state (i.e. uniformly unbent rows) as their narrow dimension grows, but they do so via distinct patterns. Tall ribbons exhibit large shears between rows at their free edges over a boundary zone of size ∼λ_x, with rows nearly fully unbending away from this boundary zone. In contrast, flattening in wide ribbons is mediated by a competition between bending and intra-row compression, with shears vanishing along most of their long dimension. Unlike shears in tall ribbons, the differential compression along rows in wide ribbons is not confined to a finite boundary layer, and instead extends through the thickness. Another key distinction is the fact that row bending at the center of the wide ribbons remains non-zero and continues to decrease as the ribbons grow in their narrow dimension, whereas row bending in tall ribbons falls to zero away from the shear boundary layer. Fig. 4c shows this contrasting behavior for row curvature behavior along the central horizontal row of tall and wide ribbons, with the former relaxing to uniformly unbent rows away from the boundary zone while the latter adopts a non-zero value that continuously unbends with increasing thickness.

These distinct mechanics are captured in the analytical forms of the elastic energy in the infinite width or height limits of the continuum solutions. The two morphologies of interest are those that grow tall with H ≫ W or wide with W ≫ H. In the limit of H/W → ∞, the excess energy of tall ribbons simplifies to

which was derived previously in ref. 29. Notably, the exponential form of the width dependence is a consequence of the finite shear boundary at the sides of the ribbon. That is, the boundary conditions required that the curvature at the free edges of the domain achieve their preferred value, but as one moves away from the free edges, the resistance to inter-row shear exponentially screens the curvature beyond a distance comparable to λ_x. Hence, local arrangements within the core of the domain become effectively insensitive to size for tall ribbons with W ≫ λ_x. This is consistent with the asymptotic form indicating that the entire width of the ribbon is unbent, with the exception of a finite boundary layer proportional to λ_x where the elastic penalty of inter-row shearing is consequently lower.

In the limit of wide ribbons W/H → ∞, the excess energy takes on the form

which corresponds to nearly uniformly-curved (i.e. circular) ribbons characterized by a mid-row curvatureA detailed calculation for this can be found in Appendix D. A key distinguishing feature of the wide ribbons is their mesoscale curvature which varies with the finite thickness of the domain itself. This mesoscale curvature results in a large size scale since κ_c(h)⁻¹ ≥ κ₀⁻¹ and we generally consider the limit of small curvature frustration where κ₀H ≪ 1. To clarify, due to the finite curvature “wide” ribbons, it is not strictly possible to extend W beyond 2πκ^(wide)_c(h)⁻¹, as the free ends of the ribbons will meet. However, we expect that in the small frustration limit (i.e. κ₀a ≪ 1) the elastic energy of the small-θ expansion accurately describes the circumferentially-uniform strain distributions in closed annular ribbons. Notably, the thickness dependence of row-curvature at the core of domains is power-law as opposed to exponential, with κ_c(h) ∼ h⁻² indicating that even for h ≫ 1, residual curvature propagates to the center of the domain and is not screened by a free boundary (Fig. 4c, inset). While for narrow assemblies (i.e. w ≪ 1 or h ≪ 1), both types of assemblies exhibit the same quadratic growth of excess energy with narrow dimension ( and ), the lack of boundary screening in wide-ribbon assemblies leads to a distinct and more rapid power-law approach to the asymptotically flattened limit .

The analytic forms of excess energy for tall and wide limits are shown as black lines in Fig. 4b and are compared in the inset, highlighting the respective similarity and distinctions in the narrow and thick ribbon limits. We next consider the thermodynamic consequences of those differences in the thickness dependence of excess energy for self-limitation.

3.2 Self-limiting domain size selection: ψ = 0 case

We now consider the consequences of anisotropy in the excess energy landscapes on the selection of self-limiting domain dimensions arising from the competition between the cohesive costs at the boundaries and the intra-domain elastic cost due to frustration. These self-limiting domains are described by minima in the per subunit assembly free energy, eqn (12). We first consider the case of isotropic binding strengths along orthogonal directions, Σ_x = Σ_y, and plot the per-subunit assembly free energies for both the discrete and continuum models in Fig. 5a. As can be seen, there are two low-energy channels through which the assembly free energy can decrease: one for fixed W and H → ∞ corresponding to tall ribbons and another for fixed H and W → ∞ corresponding to wide (annular) ribbons. Broader considerations of the ψ = 0 assembly free energy landscapes for arbitrary values of Σ_x and Σ_y reveal that these two ribbons states, i.e. tall (extending along ê₂) or wide/annular (extending along ê₁), are the only possible minima. In addition, these low energy states assemble unlimited along their long dimensions while the narrow dimensions may be finite depending on the value of the cohesive energy of their free edges. The thermodynamics of the narrow dimensions are determined by the (per subunit) assembly free energy,where σ_x and σ_y are rescaled edge binding energies given byThese rescaled binding energies can be interpreted as the ratio of the cohesive cost of having free boundaries to the cost of flattening the boundary of a ribbon over the bend-penetration lengths λ_x or λ_y for tall or wide ribbons, respectively. The free energies of tall (blue) and wide (orange) assemblies are plotted in Fig. 5b along with the local minima indicated by the empty circles. The self-limiting sizes can be determined by locating the local minimum of f^(tall)(w) and of f^(wide)(h) for given binding strengths σ_x and σ_y, plotted in Fig. 5c. In general, the self-limiting sizes increase with edge binding strength, which can be understood by a simple consideration of the small cohesion and consequently small self-limiting size limit of the continuum theory. In this limit, the balance between the quadratic growth of excess energy with size and the inverse dependence of edge cost with size leads to the power laws w_* ≃ (3σ_x/2)^1/3 and h_* ≃ (3σ_y/2)^1/3, which also shows that self-limiting domain sizes decrease with frustration as ∼κ₀^−2/3 in this weak cohesion regime.

Fig. 5 Assembly free energy and self-limiting sizes of WJ particle assemblies. (a) Free energy landscapes for with isotropic elasticity k_x = k_y for simulations and σ_x = σ_y = 0.8 for theory. (b) Free energies for tall H/W → ∞ (left) and wide W/H → ∞ (right) assemblies for varying cohesions 0 ≤ σ_x < 1.2 and 0 ≤ σ_y < 1.2. Solid dots are simulations of M × 100 and 100 × N assemblies with k_x/k_y = 10. Empty circles indicate positions of the local minima with dashed lines showing the flattening energy. (c) Dimensionless self-limiting sizes of tall (left) and wide (right) assemblies as functions of cohesion. The inset shows the centerline curvature of wide assemblies.

Beyond this weak cohesion and narrow thickness limit, tall and wide ribbons exhibit distinct thermodynamic dependencies of finite domain thickness on cohesion, most notably at the self-limiting to bulk transitions. As described previously,³³ the nature of the self-limiting to bulk transition depends on the specific form of the saturation of excess energy in the limit of large domain sizes. As described in the previous section, the elastic energies of tall and wide ribbons approach that of the bulk state with leading corrections −1/w and −1/h² respectively, differences that characterize a relatively more gradual vs. rapid expulsion of gradients in the ground states for increasingly larger domain thicknesses.

Tall ribbons exhibit a second-order transition between finite-width ribbons and bulk states in which w_* increases continuously with σ_x up to a critical value σ^(tall)_c = 1 where it diverges continuously. This is due to the fact that as the edge cohesive strength is increased, the local minimum of the free energy for tall assemblies remains below the energy of the bulk flattened (w → ∞) state and approaches that state continuously as σ_x → σ^(tall)_c (Fig. 5b, left). In Fig. 5c, we compare the analytical prediction of the self-limiting sizes from continuum theory to those obtained from the discrete WJ simulations for a range of values of elastic anisotropy k_x/k_y. Notably, this ratio controls the size of the bend-penetration λ_x relative to the discrete WJ particle size a. As k_x/k_y is increased, row bending becomes increasingly stiffer than inter-row shearing, leading to a much longer range of frustrated curvature propagation. The effect is that the size scale over which elastic energy accumulates and consequently the self-limiting sizes proportionately increase with the elastic anisotropy k_x/k_y or λ_x, as seen by the collapse of the data in Fig. 4b and 5c. This increase in size range results in the tall ribbons exhibiting more distinct, integer values of finite width prior to the divergent region for σ_x ≲ σ^(tall)_c. While in principle the equilibrium self-limiting width of tall ribbons is predicted to grow arbitrarily large, , as the strength of cohesion approaches σ^(tall)_c, it is likely that the combined effects of width-fluctuations²⁹ and the extreme sensitivity of size near σ^(tall)_c allow any practical control of width in physical experiments.

Fig. 5c shows the corresponding prediction for the growth of the self-limiting thickness h_* of wide assemblies with σ_y, which exhibits a first-order transition between self-limiting and bulk states. This means that above the transition value , the local minimum of f^(wide)(h) exceeds the free energy of the bulk state. We refer to this local minimum as being metastable, that is, while the wide assemblies can be locally stable about the energy minimum at a self-limiting thickness h_*, there is another state that is more energetically favorable, in this case, the bulk unlimited state. This metastable branch of self-limiting annular ribbons eventually loses stability to the bulk uniform state above an edge cohesion σ^(wide)_s = 9/8. This first-order transition implies that there is an upper limit to the equilibrium self-limiting thickness of annular ribbons given by . Notably, within our discrete WJ particle model studied here, the elastic length scale λ_y = v is independent of the elastic interaction stiffnesses and is therefore less than or equal to particle size a since v ≤ a, which implies that the size of equilibrium self-limiting annular assemblies is comparable to that of the particle dimension a and is therefore comparatively more limited than tall ribbons, whose size selection can vary with k_x/k_y. While not the goal of this study, it is possible to access a broader range of bend to stretch ratios through additional different microscopic features, for example, by introducing specific repulsive binding sites on the faces of WJ particles, as described previously.²⁹ Notwithstanding the limited range of equilibrium self-limiting thicknesses of wide ribbons in the present model, we compare the dependence of the self-limiting thicknesses with σ_y of discrete WJ simulations for annular ribbons to the continuum theory predictions in Fig. 5c, right.

Finally, we note that the predicted increase in self-limiting thickness with cohesion is accompanied by a corresponding decrease in the mesoscopic curvature of the ribbon itself, due to the increasing differential compression within the rows. The predicted variation of ribbon curvature κ_c(h_*) with σ_y, as well as that obtained from discrete WJ simulations, is plotted in the inset of Fig. 5c. This shows that the mesoscopic radius of curvature varies from the preferred value κ₀⁻¹ for σ_y → 0 up to twice that value at the transition point at σ_y = σ^(wide)_c. Hence, the mesoscopic dimensions of annular ribbon morphologies exhibit distinct thermodynamic sensitivity to inter-particle cohesion. In particular, changes in the domain thickness at the particle scale manifest in large changes in the gross morphology of wide ribbons at the size scale κ₀⁻¹ ≫ a.

3.3 Phase diagram: ψ = 0 case

As described above, among the rectangular domains for ψ = 0 misfit polarization, there exist two classes of minima describing tall and wide ribbons. These states differ qualitatively in terms of their gross morphology, internal strain distributions, and elastic costs of frustration. As a consequence, even for the case of equal scaled edge cohesions, σ_x = σ_y, their ground state free energies differ, with the free energy of tall ribbons generically lower than that of wide ribbons, which are therefore metastable for this isotropic binding case.

As shown in Fig. 5b, the ground state free energies and are governed by binding strengths along distinct edge directions σ_x and σ_y, respectively. Therefore, the thermodynamic selection between these distinct states is influenced by the ratio of binding along distinct directions, reflecting a preference to extend a domain along a low-edge energy face. The phase diagram of ground state morphologies is shown in Fig. 6, showing that WJ assembly exhibits polymorphism among three mesoscopically distinct states: tall, wide, and bulk assemblies.

Fig. 6 Phase diagram for tall (H > W) and wide (W > H) assemblies. Tall assemblies are energetically stable in blue regions, wide assembles are energetically stable in orange regions, and bulk assemblies are stable in gray regions. The patterned regions indicate where at least one other structure is metastable i.e. that structure, while having a local minimum in free energy, is not the lowest energy state.

The phase diagram is divided into six main regions, determined by the existence of a local minimum and whether that minimum is the ground state. If it is not the ground state, then it is metastable, meaning that while it is possible to stabilize that state, the assemblies generally prefer the lower energy ground state structures. Recall that the critical binding strengths at which the local minima of tall and wide assemblies disappear are σ_x = 1 and σ_y = 9/8, respectively. Therefore, when binding is too strong (i.e. σ_x > 1 and σ_y > 9/8), bulk assembly of unlimited flattened 2D structures occurs, as indicated by the solid gray region (vi). Decreasing the horizontal (ê₁) binding strength below its critical value (σ_x < 1), results in tall straight ribbons extending unlimited along the vertical (ê₂) binding directions, as indicated by the solid blue region (i). Lowering the vertical binding strength below its critical value (σ_y < 9/8), a local minimum develops in the free energy of the wide annular ribbons. For σ^(wide)_c < σ_y < 9/8, this minimum is metastable since the lowest energy state is still that of tall straight ribbons or bulk assembly, as indicated by the patterned blue (ii) or gray (v) regions. When the vertical binding strength is sufficiently weak (iii), finite thickness annular ribbons become the equilibrium state while the tall straight ones become metastable. The equilibrium transition between the two morphologies occurs when σ_y < σ^(wide)_c and indicating that for equal binding strengths σ_x = σ_y, tall ribbons are generally favored. That is, in the absence of cohesive anisotropy, the polarization of misfit and its elastic effects bias the thermodynamics towards mesoscopically tall straight ribbons, owing to the relatively “softer” accumulation of elastic energy of tall ribbons with width. Despite this, we note that only a relative small bias in binding energy along the horizontal and vertical directions is needed to stabilize the otherwise metastable wide and mesoscopically curved ribbons.

4 Misfit polarity along off-lattice directions: 0 < ψ < π/4

In Section 3, we found that the emergent morphologies of frustrated WJ assembly derive from two sources of anisotropy: (i) the binding strengths along the two binding directions ê₁ and ê₂ and (ii) the polarization of the WJ shape misfit itself. In this section we consider a more general class of ribbon morphologies where the misfit polarization κ₀ is not aligned to one of the two lattice directions (ê₁ or ê₂), i.e. the cases where 0 < ψ < π/4. Note that ψ = π/4 represents the largest possible deviation from a nearest-neighbor binding direction since ψ < 0 or ψ > π/4 simply correspond to mirror reflections of the WJ particles and hence do not change any results. For the ψ = 0 case, we found that the elastic cost of frustration is lowest for tall ribbons, that is, when assembly grows arbitrary long along the misfit polarization direction. In that case, the bias of elastic energy towards forming tall ribbons is compatible with the bias of low-edge energies towards forming 2D rectangular domains, which are lowest when the direction of free edges is aligned with a nearest-neighbor binding direction. More generally, assembly of finite-width ribbons could occur in off-lattice directions, a possibility which we characterize by introducing the angle ϕ between the long-axis of the ribbon and the ê₂ binding direction in the assembly, which we show schematically in Fig. 7. Notably for ϕ ≠ nπ/2 (where n denotes an integer), free edges are “terraced” leading to a greater cohesive cost per unit length of the ribbon compared to the flat edges of rectangular domains discussed in Section 3, a microscopic effect underlying facet formation more generally found in crystal assembly.⁴⁴

Fig. 7 Assembly along off lattice directions. (a) Definition of the angle ϕ between the long-axis of a ribbon and the ê₂ binding direction in the reference state. (b) Relaxed assembly for off-lattice direction ϕ = π/8 and shape asymmetry ψ = π/4, with curvatures κ_xa = κ_ya = 0.405.

In this section we analyze the thermodynamic competition between assembly of finite domains along low-edge energy directions of the 2D WJ “lattice” (i.e. ϕ = nπ/2) versus assembly along the direction favored by misfit polarization (i.e. ϕ = ψ). We first consider the assembly of free energy landscapes for domain edges aligned to the lattice directions ê₁ and ê₂, but for misfit polarities which are off those lattice directions (i.e. ϕ = 0 and 0 < ψ < π/4). We then consider possible ribbon ground states whose free edges are oriented away from the low-edge energy directions (i.e. ϕ ≠ nπ/2). For simplicity, we restrict our analysis in this section to the case of isotropic elastic parameters, meaning that k_x = k_y and correspondingly B_x = B_y, Y_x = Y_y and λ_x = λ_y.

4.1 Thermodynamic landscapes for ϕ = 0 and 0 < ψ < π/4 domains

In Fig. 8 we compare the per-subunit assembly free energy landscapes from discrete WJ simulations with isotropic binding Σ_x = Σ_y for three cases of misfit polarization, ψ = 0, π/8 and π/4. Per the design in Fig. 2, these changes in polarization are affected by appropriate changes in the taper angles, α_x and α_y, leading to kite-shaped quadrilaterals for ψ ≠ nπ/2. Upon inspection, it is clear that the gross features of the landscape are preserved with changes of polarization. Namely, the existence of two low-energy basins, oriented vertically and horizontally, represent two classes of finite-thickness ground states whose configurations have one dimension that extends unlimited. More careful inspection shows that the asymmetry between the vertical and horizontal basins is progressively reduced as the misfit polarization is increased from ψ = 0. This change is also reflected in the ground state shapes of the two ribbon classes. For ψ = 0, vertical ribbons are mesoscopically straight, while horizontal ribbons are highly curved. Upon increasing the misfit polarization to ψ = π/8, vertical ribbons adopt a slight bend, while horizontal ribbon curvature is slightly reduced. Finally, for maximal off-lattice polarization of ψ = π/4, the two ribbon morphologies adopt the same mesoscopic curvature. Indeed, in this case there is clear symmetry in the assembly free energy landscape about the M = N line. This emergent symmetry in the landscape implies that as ψ → π/4 the two ribbon morphologies become degenerate.

Fig. 8 Free energy landscapes for frustration directions ψ = 0, π/8, π/4 with fixed frustration κ₀a = 0.1 and isotropic elastic constants k_x = k_y. Example structures shown for the low energy channels in the energy landscapes and the saddle point.

This basic dependence of the ground state thermodynamics on the polarization direction can be understood simply in terms of the properties of the solutions of the continuum theory. From inspection of eqn (8) and (9), it is clear that positional and orientational strains, u and θ respectively, are linearly coupled to each other as well as to the WJ particle frustration κ₀via the boundary condition of eqn (11), which for isotropic elasticity reduces to (∇θ − κ₀T̂)·n̂ = 0 at the free edges of the domain. Hence, for ψ = 0 this requires preferred row bending at the vertically oriented free edges and vanishing column bending at the top and bottom of the domains, leading to a solution in which u and θ are linear functions of κ₀ and a corresponding elastic cost that is quadratic in preferred row curvature κ₀. Likewise, the case of ψ = π/2 corresponds to the same set of solutions just rotated by π/2, i.e. using the results of eqn (15), . For the case of ψ ≠ nπ/2, misfit polarity N̂ (and T̂) have projections along both row and column directions resulting in preferred curvatures κ₁ = κ₀ sin ψ and κ₂ = κ₀ cos ψ, respectively. In this case of off-lattice polarization, the solution for the positional strains is simply a linear superposition of the two solutions for ψ = 0 (vertical polarization) and ψ = π/2 (horizontal polarization). A more subtle, but valuable, feature of this superposition is that the elastic energy for misfit polarization along orthogonal directions in rectangular domains is uncoupled. That is, the excess energy for a general ψ value can be rewritten using the solution for the ψ = 0 elastic energy eqn (15) as

Formally, this result follows from the orthogonality of strain modes in the series solutions detailed in Appendix C.1. Intuitively, we can understand this as a result of a simple symmetry argument that elastic energy must be invariant under the reflection of the misfit polarization or shape through either the horizontal or vertical axis since mirroring the WJ particle shapes and domains will not change the (ground state) energy of the assemblies. This argument implies that there can be no terms in the elastic energy proportional to the product of column and row curvatures κ₁κ₂, and that the only terms allowed at quadratic order are proportional to κ₁² and κ₂², resulting in the form of eqn (23).

Taken together, these arguments imply that both the ground state shapes and elastic energy behave as simple mixtures of the “vertical” and “horizontal” solutions when misfit polarization is not aligned to 2D lattice directions of the WJ reference domain. Hence, this implies that the excess energy itself becomes symmetric with respect to the exchange of wide and tall directions (i.e. under (w,h) → (h,w) for ψ → π/4 where κ₁ = κ₂). It is most useful to consider the implications of this superposition in terms of the infinite ribbon domains given by

with corresponding mesoscopic curvaturesThese results show that the off-lattice misfit polarization leads to a mixing of the two distinct elastic modes for anisotropic ψ = 0 domains: the softer shear-bend mode in tall and mesoscopically-straight domains and the relatively stiffer stretch-bend mode in wide and mesoscopically-curved domains. As a consequence, for a fixed domain size, the ground-state energy and mesoscopically-curved shapes of tall and wide domains are most distinct when the polarization is aligned to one of the two binding directions, and are most similar and degenerate when the polarization is equally vertical and horizontal at ψ = π/4.

4.2 Off-lattice ribbon domains (ϕ ≠ 0)

The prior sections show that two competing mechanisms of anisotropy selection underlie self-limiting ribbon domains of the WJ assembly model. On one hand, anisotropy of binding tends to favor domains that grow longest along the relatively high binding energy directions, a generic effect of anisotropic edge energies in 2D assembly. On the other hand, elastic costs of frustration are generically lowest for domains that extend along the misfit polarization direction N̂. In Section 3, we showed that for ψ = 0 WJ particles the elastic bias for “tall” ribbons stabilized these ribbons over “wide” (annular) ribbons even for binding that is weakly favorable for the latter (i.e. σ_y ≲ σ_x). Here we consider and analyze the broader possibility that polar frustration may select ground state domains that extend along off-lattice directions, i.e. other than the ê₂ (ϕ = 0) or ê₁ (ϕ = π/2) binding directions.

To address this we consider the per-subunit assembly free energy of a ribbon whose long axis extends along a direction that makes angle ϕ relative to ŷ in the 2D reference lattice (see Fig. 7),

where [small script l] characterizes the thickness of the narrow dimension of the ribbon (perpendicular to the long dimension). For off-lattice directions 0 < ϕ < π/2, the edge energy per unit length of the ribbon varies with ϕ according to the cohesive cost breaking additional bonds when cutting a lattice along directions that are not aligned to nearest-neighbor binding directions.⁴⁴ Based on this generalized ribbon model, we plot in Fig. 9 the per-subunit free energy landscape as a function of finite-width [small script l] and orientation ϕ for several cases of misfit polarization ψ and binding anisotropy σ_x/σ_y.

Fig. 9 Energetics of off-lattice growth for ribbons with edges at an angle ϕ with respect to the vertical binding direction and thickness [small script l]. (a)–(c) Increasing misfit polarization angle ψ = 0, π/8, and π/8. Example assemblies are for the local minima ϕ = 0 and ±π/2 and saddle point near ϕ = ±π/4. (c)–(e) Increasing ratio of binding energies with . (f) Phase diagram of the transition between self-limiting and bulk growth for varying frustration polarization direction ψ and edge cohesive strength σ_⊥. The assembly free energies are plotted on a fixed range f/f_∞ ∈ [0.7,1.02] with white regions outside that range.

First, consider Fig. 9a–c which show a sequence of energy landscapes for ψ = 0, π/8, and π/4 for isotropic binding σ_x = σ_y = 0.5. All three cases show the existence of two local minima in the ([small script l], ϕ) plane, corresponding to tall, vertically-oriented (ϕ = 0) and wide, horizontally-oriented (ϕ = ±π/2) ribbons with finite thickness. Notably, even for 0 < ψ < π/4 these local minima do not shift away from the lattice directions, despite the fact that the elastic costs of frustration are lower (for a given thickness [small script l]) for ribbons extending along the polarization direction ϕ = ψ. Hence, all three landscapes are characterized by an intermediate unstable ribbon morphology at saddle points near ϕ = ±π/4, shown for each of the three cases.

In Fig. 9c–e, we consider the energy landscapes for fixed off-lattice polarization ψ = π/4 with increasing binding anisotropy, which show that the local minima corresponding to tall and wide ribbons can still exist, although for sufficient anisotropic binding the otherwise metastable wide ribbon morphology ultimately disappears, becoming thermodynamically unstable to bulk assembly in that direction.

Taken together, these assembly free energy landscapes demonstrate that changes in misfit polarization influence the ground state thermodynamics and gross morphological features of the self-limiting ribbon morphologies. While ribbon morphologies that extend along the misfit polarization direction (ϕ = ψ) tend to have lower elastic costs, these results show that the possible orientations of equilibrium and metastable morphologies appear to be strongly locked to the underlying 2D lattice directions due the anisotropic biases of edge energy. Mathematically, the relatively stronger effect of edge energy derives from its sharp and non-analytic dependence near to local minima (i.e. |ϕ − nπ/2|) in comparison to a much softer harmonic angular dependence (i.e. (ϕ − ψ)²) from elastic effects. In Fig. 9, we plot the generalized equation of state , the self-limiting ribbon width, where σ_⊥ is defined as the scaled edge energy for the free edge perpendicular to the ribbon's growth direction. This shows that for any ψ ≠ 0 or π/4 the transition between self-limiting and bulk states becomes at least weakly first order at critical cohesive strength that decreases continuously from σ_c = 1 at ψ = 0 to a lower value at ψ = π/4.

5 Discussion and conclusions

In this article, we studied the ground-state thermodynamics of the 2D planar frustrated self-assembly of polar misfitting warped-jigsaw particles via simulations of discrete particle arrays and continuum elasticity theory. Our analysis reveals three key results regarding how polar frustration effects the elasticity and morphologies of frustrated (WJ) assemblies, which we summarize as follows:

• Polar frustration set by the shape of WJ particles breaks the symmetry of the underlying frustration-free binding order, resulting in assemblies with distinct mesoscopic morphologies (in our case tall-straight and annular-like ribbons) and elastic responses to frustration (different elastic length scales and approach to bulk unlimited states) when grown parallel or perpendicular to the frustration polarization direction.

• The possible directions of growth are strongly locked to the underlying lattice binding directions due to the cost of free boundaries, and the polar frustration in our case is not sufficient to break away from those directions.

• For potential experimental realizations of WJ particles, particle designs that can control the ratios of bending to stretching, particularly and in our case, are key to assembling self-limiting structures with well-controlled sizes much larger than that of a single particle.

To be more precise, our results illuminate how anisotropic mesoscopic morphologies in the form of self-limiting multi-particle domains are selected by the distinct anisotropic features of the WJ particles themselves. We showed that polarity of the misfit direction is transmitted to the multi-particle scale such that the ground state shapes and thermodynamics of ribbons that grow “tall” versus “wide” relative to the misfit direction are mesoscopically distinguishable, each characterized by distinct mechanical modes of stress accumulation with the finite thickness of the domain. While there are those distinct morphologies oriented relative to the misfit direction, we find that much like with standard unfrustrated assembly systems, such as the Wulff construction of 2D crystals,⁴⁴ the anisotropic inter-particle binding biases assemblies to grow relatively larger in directions corresponding to low edge energies (in our case along the 2D binding lattice directions).

WJ assembly thermodynamics are distinct from other classes of misfit symmetry, such as scalar or nematic, in that the equilibrium self-limiting morphologies are imprinted with a distinctly polar character of the particle shape at the local scale. For example, tall ribbons of ψ = 0 WJ particles exhibit asymmetric shapes at their top and bottom ends, while wide ribbons exhibit mesoscopic “upwards” or “downwards” bending, which itself is sensitive to the ribbon thickness. Mesoscopic polarity arises because the microscopic asymmetry of WJ misfits breaks the otherwise 2-fold symmetry of the 2D domain anisotropy, unlike scalar or nematic misfit symmetries which preserve the underlying 2-fold symmetry. As a result of the incommensurate symmetries of polar misfit and the background 2D order of assembly, ground state thermodynamics of WJ assembly exhibits a rich polymorphism, with multiple mesoscopically distinct self-limiting morphologies. Shapes and thermodynamics of WJ assembly are controlled not only by the ratios of cohesion to elastic cost, as in standard GFA models, but also by the relative orientation of misfit polarity and 2D assembly order.

The WJ model and the dependence of the mesoscopic morphologies on misfit polarity provide new design parameters for experimental efforts to realize and control self-limiting assembly via microscopic misfit particle features such as shape and interactions. The necessary features to implement WJ assembly are directly implementable, for example, using DNA origami based particles^25,45,46 or even so-called magnetic hand-shake particles.^47,48 Both programmable subunit platforms combine the necessary ingredients to engineer and control the 2D misfitting shapes (e.g. trapezoids) of nanoscale subunits as well as their specific interactions needed to promote oriented edge bending, whose relative strengths could be dynamically tuned to stimulate transitions between the mesoscopically distinct self-limiting polymorphs.

Of course, the present theoretical study is limited in several respects that may be key to connecting with specific experimental systems. Here we note two primary limitations which require further study. Foremost, our study is restricted to ground states, essentially T = 0 thermodynamics, whereas equilibrium assembly necessarily requires finite-temperature effects as well as finite range interactions between subunits. Notably, the equilibrium self-limitation behavior is predicted to depend on the ratio of inter-particle cohesion to inter-particle elastic stiffness, both features of which are mediated by the interactions that bind neighbor particles together. As argued previously for simple models of inter-subunit binding,^29,30 equilibrium self-limitation places restrictions on the possible range of interactions. For example, the thermodynamics of tall ribbons with self-limiting widths that are multiple WJ particles wide requires k_x/k_y > 1, which can be achieved by introducing relatively short range binding in the ê₁ direction to make rows stiffer to bend while introducing longer range binding in the ê₂ direction to allow for softer shearing between rows. How much the equilibrium windows of tall ribbon assembly can be extended by controlling interactions with anisotropic range, not to mention the potential kinetic limitations of finite-range interactions on reaching those equilibrium states, remains to be studied. In addition, we have neglected excluded volume in our analysis of the ground states for our purposes of studying the elastic effects of polar frustration, although again such effects where included in the prior “discrete-particle” numerical study of WJ particles,²⁹ and led to no clear discrepancies with continuum theory. While at least for sufficiently low-frustration it is clear that excluded volume could be neglected, in more general contexts its effects may not be so fully eliminated in some realizations of WJ particles. For example, excluded volume effects may result in elastic moduli that are stiffer or more non-linear than in the present model as well as the presence of defects in the assemblies when certain particle arrangements are inhibited, and therefore deserves further study.

A second major open issue is the stability of the 2D planar domains to out-of-plane deformations and what kinds of complex 3D morphologies emerge. As an example, we show in Fig. 10 a simple extension of our planar WJ model that allows out-of-plane bending at a finite bending stiffness. As an example, we achieve out-of-plane stiffness by creating WJ particles with thickness along their ê₁ × ê₂ axis and adding an addition “layer” of springs to their top and bottom faces, resulting in an out-of-plane bending modulus, as shown in Fig. 10. Notably, we observe that both wide and tall ribbon domains exhibit complex 3D buckled shapes when the out-of-plane bending stiffness falls below an apparent threshold. What determines the threshold stiffness to suppress out-of-plane buckling in WJ domains, how the microscopic features such as the shape and interactions of WJ particles affect the out-of-plane morphologies, and in particular how misfit frustration is relaxed by deforming into the third dimension remain to be explored.

Fig. 10 Three-dimensional version of WJ particles (top) with a thickness t. Examples of out-of-plane deformations of tall and wide assemblies (bottom) when the out-of-plane bending modulus B_α⊥ ∼ k_αt² is decreased.

Notwithstanding these limitations on our current understanding of WJ assembly, this study points to a more general classification of GFAs in terms of the misfit symmetry of the subunits and its resulting effects at the mesoscale. Such a framework may be useful for analyzing the ingredients and outcomes of GFA for a much broader class of misfitting particles. For example, ref. 23 and 24, through a study of misfitting elastic polygons, have noted that for multi-protein assembly some measure of shape misfit is more likely the rule rather than the exception. It is reasonable to expect that misfitting shape in such models could be mapped onto a multi-pole expansion in the misfit symmetry, and further that mesoscopic assembly outcomes (e.g. self-limiting domain size and shape) may be systematically correlated to the corresponding moments of misfit symmetry at the particle scale. We should emphasize, however, that the resulting morphologies and ground-state thermodynamics derive from the combination of both the symmetry of the building blocks as well as the symmetries of the underlying assembly, which for small enough frustration derive from the symmetries of the unfrustrated order. The way in which (frustrated) particle symmetries can interact with elastic and cohesive responses of an assembly may in fact result in similar symmetry-breaking states for distinct classes of frustrations. Consider, for example, the cases of scalar and hexapolar frustration shown in Fig. 1a and d, both of which are characterized by underlying hexagonal 2D crystal states of unfrustrated order. The combination of the elastic response and orientational dependence of domain energies on line-energies of a hexagonal lattice already select three possible directions for narrow or wide domain growth, even for the case of scalar frustration. Hence, the 3-fold symmetry of a hexapolar misfit particle, such as in the case of a flexible hexagon of alternating short-long edges,^23,24 is not expected to lower the symmetry of self-limiting ground-state thermodynamcs below that observed already for scalar misfit, although it may lead to some additional polarity of the ribbon morphologies themselves.

This basic notion of misfit symmetry, and the complexity of potential outcomes, suggests the possibility of a generic, field-theoretic (e.g. Landau theory) approach to classifying and studying the interplay between inter-particle order and distinct symmetries of frustrated units. Notably, recent progress along these lines has been made in the context of nematic liquid crystalline phases, which can be frustrated by a combination of distinct “modes”:⁴⁹ splay (polar), bend (polar), twist (pseudoscalar), and bi-axial splay (tensorial). The general consequences for these distinct modes for finite-size domain selection in 3D nematically ordered assemblies remains to be understood,^50,51 notwithstanding the much broader possibilities offered by distinct phases of long-range order (e.g. nematic, columnar, crystalline) in the assembled state.

Data availability

Data for this manuscript include numerical energy minimization codes for discrete “warped-jigsaw” particle assemblies, as well as example parameter inputs for generating the data analyzed in the plots. Additionally, these include Mathematica notebooks and numerical data used to generate plots and schematics in the manuscript. These data are available at https://hdl.handle.net/20.500.14394/56645, Scholarworks repository, hosted by UMass Amherst.

Conflicts of interest

There are no conflicts to declare.

Appendices

A Simulation details and energy relaxation

We model our discrete WJ particles as rigid bodies composed of pairs of elastic bonding sites on each of their four edges. Let the position and frame of the particle be described by R and {ê₁, ê₂}, respectively, as shown in Fig. 2c. The attractive sites of the top, bottom, left, and right edges of the quadrilateral-shaped particle are respectively located atwhere α₁ and α₂ are the taper angles of the left-right (ê₁) and top-bottom (ê₂) edges, respectively. The ± represents the two interaction sites on each edge. To simulate sheets of such particles, we arrange them on a reference square lattice at positions R^(m,n) = max̂ + naŷ with the particle frames aligned with the lattice directions {ê₁, ê₂} = {x̂, ŷ}, and bind the specific binding sites together (s^(m,n+1)_bottom,± with s^(m,n)_top,± and s^(m+1,n)_left,± with s^(m,n)_right,±) with zero rest length Hookean springs with stiffnesses k₁ and k₂ between the left-right and top-bottom edges of neighboring particles, respectively.

Once a pre-assembled sheet of WJ particles with a given connectivity is created, we can relax the structure by computing the forces and torques acting on each particle and translating and rotating them according to overdamped equations of motion. To be more precise, the position R of a particle is updated via

R(t + δt) = R(t) + δtf_total (30)

where is the sum of the forces acting of the elastic binding sites of the particle and δt is the step size. The orientation θ of the particle, defined by cos θ = ê₂·ŷ, is updated via

θ(t + δt) = θ(t) + δtτ_total·ẑ (31)

where is the sum of torques due to the forces on the elastic binding sites and r_site is the position of an elastic binding site relative to the center of the particle.

B Continuum elasticity of WJ particle sheets

We here derive the continuum theory for describing M × N sheets of WJ particles. As described in Appendix A, we treat the WJ particles as rigid subunits with a pair of elastic binding sites on each edge. We start with an M × N sheet arranged on a reference square lattice such that particle (m,n) is located at R^(m,n) = max̂ + naŷ + u^(m,n), where u^(m,n) is a displacement relative to the reference configuration. Using the positions of the elastic binding sites (eqn (29a)–(29d)), we can compute the separations between the specific interaction sites between the top-bottom edges and left-right edges of neighboring subunits, which are given bywhere we have defined Δ_nu^(m,n) = u^(m,n+1) − u^(m,n), Δ_nê^(m,n)₁ = ê^(m,n+1)₁ − ê^(m,n)₁, and 〈ê^(m,n)₂〉_n = (ê^(m,n+1)₂ + ê^(m,n)₂)/2 (and similarly for Δ_m and 〈〉_m for the (m + 1,n) and (m,n) subunits). Assuming the interaction sites interact via Hookean springs with stiffnesses k₁ and k₂ for the left-right and top-bottom edges respectively, the energy of a sheet isThe continuum limit can be obtained by taking Δ_m/a → ∂_x, Δ_n/a → ∂_y, and resulting inwhere we identify the moduli and curvatures asand κ₁ = (2 tan α₂)/a and κ₂ = (2 tan α₁)/a, which represent the preferred bending orthogonal to the ê₁ and ê₂ directions, respectively. Defining rotations θ relative to the reference state, we have ê₁ = cos θx̂ + sin θŷ and ê₂ = −sin θx̂ + cos θŷ. In terms of the positional and orientational strains defined by

γ_ij = ∂_iu_j − [R_ij(−θ) − δ_ij], (36a)

β = ∇θ + κ₂x̂ − κ₁ŷ, (36b)

where R_ij is the 2D rotation matrix and δ_ij is the Kronecker delta, the energy becomesFor small rotations θ from the reference state of a square lattice aligned with x̂, ŷ, we may write ê₁ ≃ x̂ + θŷ and ê₂ ≃ −θx̂ + ŷ. The linearized energy functional becomesMinimizing with respect to the displacement u and rotations θ, we obtain the Euler–Lagrange equations for force balance

Y_x∂_x²u_x + Y_y∂_y²u_x = −Y_y∂_yθ (39a)

Y_x∂_x²u_y + Y_y∂_y²u_y = Y_x∂_xθ (39b)

and torque balance

B_x∂_x²θ + B_y∂_y²θ = −Y_x(∂_xu_y − θ) + Y_y(∂_yu_x + θ) (40)

The boundary conditions for the free edges along x = ±W/2 are

∂_xu_x|_x=±W/2 = 0 (41a)

(∂_xu_y − θ)|_x=±W/2 = 0 (41b)

(∂_xθ + κ₂)|_x=±W/2 = 0 (41c)

and along y = ±H/2 are

(∂_yu_x + θ)|_y±H/2 = 0 (42a)

∂_yu_y|_y=±H/2 = 0 (42b)

(∂_yθ − κ₁)|_y±H/2 = 0 (42c)

C Analytical solution and approximation

Note that due to linearity, we can solve the continuum theory by splitting it into two separate problems: one for κ₂ ≠ 0, κ₁ = 0 and the other for κ₂ = 0, κ₁ ≠ 0. We focus on κ₂ ≠ 0, κ₁ = 0 since the other case can simply be obtained by swapping indices x↔y and curvatures κ₂ ↔ −κ₁.

To start, it will be useful to briefly summarize the solution for infinitely tall ribbons with κ₂ ≠ 0, κ₁ = 0 in ref. 29, as we will use it to construct the general solution. Consider such an infinite ribbon bound between −W/2 ≤ x ≤ W/2. Since there will be no y dependence due to symmetry, the Euler-Lagrange equations become

Y_x∂_x²u_x = 0 (43a)

Y_x∂_x²u_y = ∂_xθ (43b)

B_x∂_x²θ = −Y_x∂_xu_y + (Y_x + Y_y)θ (43c)

The boundary conditions along x = ±W/2 are the same as eqn (41a)–(41c). This system of equations can easily be solved for the displacement and rotation fields, which are

u_x^(∞)(x) = 0 (44a)

Using these results, we construct the solution for a finite domain bound by −W/2 ≤ x ≤ W/2 and −H/2 ≤ y ≤ H/2 as follows. We write the displacement and rotation fields as

u_x(x, y) = u^(∞)_x(x) + U_x(x, y) (45a)

u_y(x, y) = u^(∞)_y(x) + U_y(x, y) (45b)

θ(x, y) = θ^(∞)(x) + Θ(x, y) (45c)

Substituting these into eqn (39a), (39b) and (40), we find that U_x(x, y), U_y(x, y), and Θ(x, y) satisfy the same force and torque balance equations

Y_x∂_x²U_x + Y_y∂_y²U_x = −Y_y∂_yΘ (46a)

Y_x∂_x²U_y + Y_y∂_y²U_y = Y_x∂_xΘ (46b)

B_x∂_x²Θ + B_y∂_y²Θ = −Y_x(∂_xU_y − Θ) + Y_y(∂_yU_x + Θ) (46c)

with the boundary conditions

∂_xU_x|_x=±W/2 = 0 (47a)

(∂_xU_y − Θ)|_x=±W/2 = 0 (47b)

∂_xΘ|_x=±W/2 = 0 (47c)

(∂_yU_x + Θ)|_y=±H/2 = −θ^(∞)(x) (48a)

∂_yU_y|_y=±H/2 = 0 (48b)

∂_yΘ|_y=±H/2 = 0 (48c)

We start with constructing Θ(x, y). Note that due to symmetry, Θ(x, y) must be odd in x. From the Θ boundary conditions (eqn (47c) and (48c)), we can writeSimilarly, U_x(x, y) must also be odd in x. From the boundary condition eqn (47a), we can writefor some function F_m(y). Substituting this into eqn (46a) and using the orthogonality of the set , we have that F_m(y) satisfiesthe general solution of which isThe boundary condition eqn (48a) tells us thatwhich givesFinally, for U_y(x, y), the boundary condition eqn (48b) allows us to writeSubstituting into eqn (46b), we find that G_n(x) satisfiesBy symmetry, U_y(x, y) must be even in x, and so the general solution for G_n(x) can be written asThe boundary condition eqn (47b) tells us thatwhich gives for n > 0Note that g₀ just gives a global translation, so we pick g₀ = 0 for simplicity. To determine an equation for A_m,n, we take the constructions for Θ, U_x, and U_y and substitute them into eqn (46c). Using the orthogonality of the sines and cosines, we have for n evenand for n odd

Note that for n odd, g_n and c_m are sums of A_m,n. By inspection of eqn (61), we can conclude that A_m,n = 0 for all odd n satisfies eqn (61). Consequently, this means that c_m = 0 for all m and g_n = 0 for odd n. Thus, we find that Θ(x, y) and U_y(x, y) are even in y and U_x(x, y) is odd in y. We summarize the series solution.

where s_m, g_n, and A_m,n for odd m and even n are given by eqn (54a), (59) and (60), respectively. Unfortunately, due to s_m and g_n depending on infinite sums of A_m,n, we effectively have an infinite system of equations for determining an infinite set of coefficients, which is impossible to solve. However, in Appendix C.2, we consider approximating the solution by truncating the series solution.

So far, we solved the case of κ₂ ≠ 0, κ₁ = 0. Fortunately, we do not need to repeat the same cumbersome calculation for κ₂ = 0, κ₁ ≠ 0. All that needs to be done is to take the κ₂ ≠ 0, κ₁ = 0 solution and simply swap all indices x ↔ y, lengths W ↔ H, and curvatures κ₂ ↔ −κ₁, and sum up the corresponding displacement and rotation fields. Putting this all together, we have that the general solution for arbitrary κ₂,κ₁ is

where , , and satisfy similar relationships as eqn (54a), (59) and (60) but with the swaps x ↔ y, W ↔ H, and κ₂ ↔ −κ₁. Note that these solutions are written so that the terms in the first line are proportional to κ₂ while those in the second line are proportional to κ₁.

C.1 Absence of κ₂κ₁ terms in the elastic energy. Naturally, we are interested in the elastic energy for general κ₂,κ₁, which should contain terms proportional to κ₂², κ₂κ₁, and κ₁². We show that the elastic energy contains no κ₂κ₁ terms i.e. the elastic energy can be written as the sum of the individual elastic energies for κ₂ ≠ 0, κ₁ = 0 and κ₂ = 0, κ₁ ≠ 0. To see this, let us write the solutions as

θ(x, y; κ₂, κ₁) = θ(x, y; κ₂, 0) + θ(x, y; 0, κ₁) (64a)

u_x(x, y; κ₂, κ₁) = u_x(x, y; κ₂, 0) + u_x(x, y; 0, κ₁) (64b)

u_y(x, y; κ₂, κ₁) = u_y(x, y; κ₂, 0) + u_y(x, y; 0, κ₁) (64c)

The parity of each piece is as follows: θ(x, y; κ₂, 0) is odd in x and even in y while θ(x, y; 0, κ₁) is even in x and odd in y; u_x(x, y; κ₂, 0) is odd in both x and y while u_x(x, y; 0, κ₁) is even in both x and y; and finally u_y(x, y; κ₂, 0) is even in both x and y while u_y(x, y; 0, κ₁) is odd in both x and y.

Let us examine each term in the elastic energy eqn (38). We suppress the x, y coordinates to keep things simple. The bending energy terms are

(∂_xθ + κ₂)² = [∂_xθ(κ₂, 0) + κ₂]² + [∂_xθ(0, κ₁)]² + 2[∂_xθ(κ₂, 0) + κ₂]∂_xθ(0, κ₁) (65)

(∂_yθ − κ₁)² = [∂_yθ(κ₂, 0)]² + [∂_yθ(0, κ₁) − κ₁]² + 2∂_yθ(κ₂, 0)[∂_yθ(0, κ₁) − κ₁] (66)

Since ∂_xθ(κ₂, 0) + κ₂ is even in both x and y while ∂_xθ(0, κ₁) is odd in both x and y, their product integrates to zero. Similarly, since ∂_yθ(κ₂, 0) is odd in both x and y while ∂_yθ(0, κ₁) − κ₁ is even in both x and y, their product also integrates to zero. For the stretching energy terms, we have

(∂_xu_x)² = [∂_xu_x(κ₂, 0)]² + [∂_xu_x(0, κ₁)]² + 2∂_xu_x(κ₂, 0)∂_xu_x(0, κ₁) (67)

(∂_xu_y − θ)² = [∂_xu_y(κ₂, 0) − θ(κ₂, 0)]² + [∂_xu_y(0, κ₁) − θ(0, κ₁)]² + 2[∂_xu_y(κ₂, 0) − θ(κ₂, 0)][∂_xu_y(0, κ₁) − θ(0, κ₁)] (68)

(∂_yu_x + θ)² = [∂_yu_x(κ₂, 0) +θ(κ₂, 0)]²+ [∂_yu_x(0, κ₁) +θ(0, κ₁)]² + 2[∂_yu_x(κ₂, 0) + θ(κ₂, 0)][∂_yu_x(0, κ₁) +θ(0, κ₁)] (69)

(∂_yu_y)² = [∂_yu_y(κ₂, 0)]² + [∂_yu_y(0, κ₁)]² + 2∂_yu_y(κ₂, 0)∂_yu_y(0, κ₁) (70)

Following the same procedure, note that ∂_xu_x(κ₂, 0) is even in x and odd in y while ∂_xu_x(0, κ₁) is odd in x and even in y; ∂_xu_y(κ₂, 0) − θ(κ₂, 0) is odd in x and even in y while ∂_xu_y(0, κ₁) − θ(0, κ₁) is even in x and odd in y; ∂_yu_x(κ₂, 0) + θ(κ₂, 0) is odd in x and even in y while ∂_yu_x(0, κ₁) + θ(0, κ₁) is even in x and odd in y; and finally ∂_yu_y(κ₂, 0) is even in x and odd in y while ∂_yu_y(0, κ₁) is odd in x and even in y. Therefore, the products of these terms integrate to zero, and so we see that the elastic energy can be decomposed into the sum of the elastic energies for κ₂ ≠ 0, κ₁ = 0 and κ₂ = 0, κ₁ ≠ 0, which are proportional to κ₂² and κ₁², respectively.

Note that this result could also have been obtained from a simple symmetry argument. Flipping the sign of κ₂ or κ₁ simply creates a subunit that is reflected across ê₁ or ê₂, respectively. All this does is create an assembly that is a reflection of the original assembly, which means that the energy must remain the same. Therefore, there cannot be any terms proportional to κ₂κ₁ in the elastic energy.

C.2 Truncated series approximation. The solution of the coupled PDEs eqn (8) and (9) requires solving an infinite system of equations eqn (60) for all the coefficients A_m,n for all n even, which therefore cannot be done in closed-form.

To acquire a closed-form expression for the elastic ground states, we perform a truncation on the series solution eqn (62). We here focus on the case of κ₂ ≠ 0, κ₁ = 0 to keep things simple since the general solution can be obtained by simply swapping indices as explained in Appendix C. To determine how to truncate the solution, we make the following observation on the behavior of u_x and u_y when going from an infinitely tall assembly N → ∞ to one that is finite in height. Note that for the infinite tall structures (eqn (44a)–(44c)), the row extension or compression is uniform all along y with u_x = 0 and vertical displacements have an exponential-like profile . When N is finite, two things occur:

• The extension and compression are no longer uniform all along y. Instead, as shown in Fig. 3, the top edge of an assembly experiences extension while the bottom edge, compression. In other words, u_x must have some nonzero y dependence, particularly concentrated near the top and bottom edges of the assembly.

• The vertical displacements u_y switch from a cosh profile where there can be a flattened core to a profile that resembles a circular arc throughout most of each row, as shown in Fig. 4d. In other words, the largest correction to u_y will be the change in its x dependence when the height is finite. Any variation along y is essentially a higher order i.e. the top rows curve similarly to the bottom rows.

These physical observations inform us about an approximation. We mostly need to account for the x dependence changes of u_y (from cosh to parabolic or circular) and the leading order dependence of u_x on y since u_x = 0 in the infinitely tall ribbons. This can be achieved by taking the leading n mode i.e. n = 0 and setting A_m,n = 0 for n > 0. The result is the following approximation

where we define A_m = A_m,0. This approximation combined with the conditions eqn (54a) and (60) for n = 0 allows us to write in closed-formand

Substituting this approximation in the elastic energy functional, we obtain

where we rescaled by the flattening energy _flat = B_xκ₂²/2.

We can check several limits. For tall assemblies with h → ∞, we find

For wide assemblies with w → ∞ (equivalently, an annulus), we findFinally for small assemblies with w, h ≪ 1, we have to leading orderWe briefly prove that this is symmetric about h = w. We use the following identity (see eqn (75))The excess energy can then be rewritten asUsingandwe finally write the excess energy in the small size limit as
This is symmetric about h = w.

D Mesoscopic curvature of wide ribbons

We here derive a discrete version of the centerline curvatures of wide or annular assemblies. Consider an annular assembly that is N particles thick. Since this annular assembly is radially symmetric, we can treat it as a ring of radial stacks of particles. Suppose that the inner ring has radius R and adjacent radial stacks are simply rotated relative to each other by an angle δθ around the center of the annulus. For simplicity (and from observation of simulations of wide assemblies), we assume that there is not significant radial compression of the rings i.e. u_r ≪ R. The elastic energy between two adjacent rows is simplyMinimizing this energy with respect to R and δθ, we haveandwhere R_c is the centerline radius. Note that for N = 1, we have that the centerline radius is R_c = a/(2 tan α) = κ₀⁻¹, which is exactly the preferred radius of a single ring of particles. For N ≫ 1, we have R_c ≃ N²a/(6tanα) = N²/(3κ₀). So, the centerline curvature scales goes as κ_c ≃ κ₀/(N²/3).

This scaling of the centerline curvature with thickness can also be obtained from a simple energy argument. The bending energy for an annulus with curvature κ is . The amount of extension or compression energy near the outer and inner rings can be estimated as . Minimizing the sum of these energies with respect to κ, we find

where and γ is some constant. Comparing with the discrete limit, we find γ = 1/3.

Acknowledgements

The authors are grateful to R. Mathew, I. Spivack, N. Hackney, J. Berengut, and L. Lee for valuable discussions about this work. This study was supported by the US National Science Foundation through awards NSF DMR-2028885, DMR-2349818 and the Brandeis Center for Bioinspired Soft Materials, an NSF MRSEC, DMR-2011846. Simulation studies of the monomer model were performed on the UMass UNITY Cluster at the Massachusetts Green High Performance Computing Center.

Notes and references

1. G. H. Wannier, Phys. Rev., 1950, 79, 357.
2. M. F. Collins and O. A. Petrenko, Can. J. Phys., 1997, 75, 605.
3. J. Vannimenus and G. Toulouse, J. Phys. C: Solid State Phys., 1997, 10, L537.
4. R. A. Reddy and C. Tschierske, J. Mater. Chem., 2006, 16, 907.
5. H. Takezoe and Y. Takanishi, Jpn. J. Appl. Phys., 2006, 45, 597.
6. C. Fernández-Rico, M. Chiappini, T. Yanagishima, H. de Sousa, D. G. A. L. Aarts, M. Dijkstra and R. P. A. Dullens, Science, 2020, 369, 950.
7. R. Ghafouri and R. Bruinsma, Phys. Rev. Lett., 2005, 94, 138101.
8. S. Armon, H. Aharoni, M. Moshe and E. Sharon, Soft Matter, 2014, 10, 2733.
9. M. Zhang, D. Grossman, D. Danino and E. Sharon, Nat. Commun., 2019, 10, 3535.
10. E. A. Matsumoto, G. P. Alexander and R. D. Kamien, Phys. Rev. Lett., 2009, 103, 257804.
11. T. Gibaud, E. Barry, M. J. Zakhary, M. Henglin, A. Ward, Y. Yang, C. Berciu, R. Oldenbourg, M. F. Hagan, D. Nicastro, R. B. Meyer and Z. Dogic, Nature, 2012, 481, 348.
12. P. Sharma, A. Ward, T. Gibaud, M. F. Hagan and Z. Dogic, Nature, 2014, 513, 77.
13. A. Aggeli, I. A. Nyrkova, M. Bell, R. Harding, A. N. S. L. Carrick, T. C. M. McLeish and N. Boden, Proc. Natl. Acad. Sci. U. S. A., 2001, 98, 11857.
14. M. S. Turner, R. W. Briehl, F. A. Ferrone and R. Josephs, Phys. Rev. Lett., 2003, 90, 128103.
15. G. M. Grason and R. F. Bruinsma, Phys. Rev. Lett., 2007, 99, 098101.
16. Y. Yang, R. B. Meyer and M. F. Hagan, Phys. Rev. Lett., 2010, 104, 258102.
17. A. I. Brown, L. Kreplak and A. D. Rutenberg, Soft Matter, 2014, 10, 8500.
18. D. M. Hall, I. R. Bruss, J. R. Barone and G. M. Grason, Nat. Mater., 2016, 15, 727.
19. G. Meng, J. Paulose, D. R. Nelson and V. N. Manoharan, Science, 2014, 343, 634.
20. W. T. M. Irvine, V. Vitelli and P. M. Chaikin, Nature, 2010, 468, 947.
21. R. E. Guerra, C. P. Kelleher, A. D. Hollingsworth and P. M. Chaikin, Nature, 2018, 554, 346.
22. S. Li, R. Zandi, A. Travasset and G. M. Grason, Phys. Rev. Lett., 2019, 123, 145501.
23. M. Lenz and T. A. Witten, Nat. Phys., 2017, 13, 1100.
24. H. Le Roy, M. M. Terzi and M. Lenz, Phys. Rev. X, 2025, 15, 011022.
25. J. F. Berengut, C. K. Wong, J. C. Berengut, J. P. K. Doye, T. E. Ouldridge and L. K. Lee, ACS Nano, 2020, 14, 17428.
26. T. Uchida, K. Abe, Y. Endo, S. Ichiseki, S. Akita, S. Liu, S. Aradachi, M. Saito, A. Fukuchi, T. Kikkawa, T. Dammaretz, I. Kawamata, Y. Suzuki, S. M. Nomura and S. Murata, Small, 2017, 13, 1702158.
27. B. Tyukodi, F. Mohajerani, D. M. Hall, G. M. Grason and M. F. Hagan, ACS Nano, 2022, 16, 9077.
28. N. Tanjeem, D. M. Hall, M. B. Minnis, R. C. Hayward and G. M. Grason, Phys. Rev. Res., 2022, 4, 033035.
29. I. R. Spivack, D. M. Hall and G. M. Grason, New. J. Phys., 2022, 24, 063023.
30. D. M. Hall, M. J. Stevens and G. M. Grason, Soft Matter, 2023, 19, 858.
31. F. Serafin, J. Lu, N. Kotov, K. Sun and X. Mao, Nat. Commun., 2021, 12, 4925.
32. S. Meiri and E. Efrati, Phys. Rev. E, 2021, 104, 024703.
33. M. F. Hagan and G. M. Grason, Rev. Mod. Phys., 2021, 93, 025008.
34. G. M. Grason, J. Chem. Phys., 2016, 145, 110901.
35. C. Sigl, E. M. Willner, W. Engelen, J. A. Kretzmann, K. Sachenbacher, A. Liedl, F. Kolbe, F. Wilsch, S. A. Aghvami, U. Prozter, M. F. Hagan, S. Fraden and H. Dietz, Nat. Mater., 2021, 20, 1281.
36. D. Hayakawa, T. E. Videbaek, D. M. Hall, H. Fang, C. Sigl, E. Feigl, H. Deitz, S. Fraden, M. F. Hagan, G. M. Grason and W. B. Rogers, Proc. Natl. Acad. Sci. U. S. A., 2022, 119, e2207902119.
37. Y. Suzuki, G. Cardone, D. Restrepo, P. D. Zavattieri, T. S. Baker and F. A. Tezcan, Nature, 2016, 369.
38. K. T. Sullivan, R. C. Hayward and G. M. Grason, Phys. Rev. E, 2024, 110, 024602.
39. M. Wang and G. M. Grason, Phys. Rev. E, 2024, 109, 014608.
40. N. W. Hackney, C. Amey and G. M. Grason, Phys. Rev. X, 2023, 13, 041010.
41. S. Schneider and G. Gompper, Europhys. Lett., 2005, 70, 136.
42. F. Cui, S. Marbach, J. A. Zhen, M. Holmes-Cerfon and D. J. Pine, Nat. Commun., 2022, 13, 2304.
43. W. B. Rogers and J. C. Crocker, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 15687.
44. C. Rottman and M. Wortis, Phys. Rep., 1984, 59.
45. D. Hayakawa, T. E. Videbæk, G. M. Grason and W. B. Rogers, ACS Nano, 2024, 18, 19169–19178.
46. J. Lee, J. Kim, G. Posnjak, A. Ershova, D. Hayakawa, W. M. Shih, W. B. Rogers, Y. Ke, T. Liedl and S. Lee, Nano Lett., 2025, 25, 16–27.
47. R. Niu, C. X. Du, E. Esposito, J. Ng, M. P. Brenner, P. L. McEuen and I. Cohen, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 24402–24407.
48. C. X. Du, H. A. Zhang, T. G. Pearson, J. Ng, P. L. McEuen, I. Cohen and M. P. Brenner, Soft Matter, 2022, 18, 6404–6410.
49. J. V. Selinger, Annu. Rev. Condens. Matter Phys., 2022, 13, 49–71.
50. J. Pollard and G. P. Alexander, New J. Phys., 2021, 23, 063006.
51. L. C. B. da Silva and E. Efrati, New J. Phys., 2021, 23, 063016.

---

Footnotes

† This classification in terms of the shape-misfit symmetry may very well be extended beyond 2D.

‡ Greek indices α, β = 1, 2 are used to denote vector components in the local particle frame, while Roman indices i, j = x, y refer to vector components in the fixed Cartesian frame.

---