DOI:
10.1039/D5SM00136F
(Paper)
Soft Matter, 2025,
21, 5423-5446
Self-limiting states of polar misfits: frustrated assembly of warped-jigsaw particles
Received
7th February 2025
, Accepted 24th April 2025
First published on 2nd May 2025
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 membranes7–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 energy32 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.29 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,29 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 morphologies19,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 frustration40 (which also exhibit similar branching morphologies as crystalline assembly on spherical surfaces), frustrated hyperbolic crystalline assemblies,30 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.29 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 {ê1, ê2} 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 (red arrows) characterizes the direction of preferred curvatures κ0i.e. the rotation of ê1 and ê2 along directions perpendicular to , described by the unit vector . (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 α1 and α2 along the horizontal ê1 and vertical ê2 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
| κ0 = κ0 , | (1) |
where
κ0 is the magnitude of preferred curvature and
![[N with combining circumflex]](https://www.rsc.org/images/entities/b_char_004e_0302.gif)
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 with combining circumflex]](https://www.rsc.org/images/entities/b_char_0054_0302.gif)
= −
ẑ ×
![[N with combining circumflex]](https://www.rsc.org/images/entities/b_char_004e_0302.gif)
, where
ẑ is normal to the
xy plane of assembly. The preferred curvatures along rows and columns are given by
κ1 ≡
κ0·
ê1 = (2
![[thin space (1/6-em)]](https://www.rsc.org/images/entities/char_2009.gif)
tan
α2)/
a and
κ2 ≡
κ0·
ê2 = (2
![[thin space (1/6-em)]](https://www.rsc.org/images/entities/char_2009.gif)
tan
α1)/
a, which describe the amount of preferred bending orthogonal to local
ê1 and
ê2 directions, respectively. Hence,

quantifies the magnitude of the misfit while
![[N with combining circumflex]](https://www.rsc.org/images/entities/b_char_004e_0302.gif)
= −sin
ψê1 − cos
ψê2 describes its orientation, which we parameterize by the angle
ψ in the local particle frame with sin
ψ =
κ1/
κ0 and cos
ψ =
κ2/
κ0.
Shown in Fig. 2e is an example of an assembly formed by trapezoidal WJ particles with ψ = 0, whose shapes are symmetric along their ê1 (horizontal) directions and asymmetric along their ê2 (vertical) directions, as characterized by the downward misfit polarity
= −ê2. 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 ê1 or ê2 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 = x
+ yŷ and particle frames ê1 =
, ê2 = ŷ. 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 ê1 = cos
θ(r)
− 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,29 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 α1, α2 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
|  | (2) |
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 by
|  | (3) |
where
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”
ê2 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 k1 = k2. | |
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) = ∂iuj − [Rij(−θ) − δij], | (4) |
where the matrix
Rij(
θ) 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) = ∇θ − κ0 , | (5) |
where
![[T with combining circumflex]](https://www.rsc.org/images/entities/b_char_0054_0302.gif)
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
γ12 =
β1 = 0 for extension and compression,
γ11 =
β1 = 0 for shear, and
γ11 =
γ12 = 0 for bend. In the analysis below, we linearize the theory for small frame rotations away from the uniform reference lattice aligned with
![[x with combining circumflex]](https://www.rsc.org/images/entities/b_char_0078_0302.gif)
,
ŷ, resulting in
| γij(r) ≃ ∂iuj − θ(r)εij; β(r) ≃ ∇θ − κ0 0, | (6) |
where
εij is the 2D anti-symmetric Levi-Cevita tensor and
0 = −cos
ψ![[x with combining circumflex]](https://www.rsc.org/images/entities/b_char_0078_0302.gif)
+ 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)
|  | (7) |
where in the small-
θ approximation we neglect the difference between {
ê1,
ê2} and {
![[x with combining circumflex]](https://www.rsc.org/images/entities/b_char_0078_0302.gif)
,
ŷ}, 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
| Yx∂x2ux + Yy∂y2ux = −Yy∂yθ; Yx∂x2uy + Yy∂y2uy = Yx∂xθ, | (8) |
and torque balance
| Bx∂x2θ + By∂y2θ = −Yx(∂xuy − θ) + Yy(∂yux + θ). | (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 forces
|  | (10) |
and torques
|  | (11) |
where
![[n with combining circumflex]](https://www.rsc.org/images/entities/b_char_006e_0302.gif)
= (
nx,
ny) 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
·
≠ 0 and vanishing θ gradients where
·
= 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.33 In our discrete WJ assemblies, the assembly free energy per subunit for a M × N domain takes the form |  | (12) |
The 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, |  | (13) |
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,34 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 κ0 is parallel to ê2 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 κ0 on the vertically oriented sides of the domains, while its value in the interior of the domain may deviate from κ0 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, |  | (14) |
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 by |  | (15) |
where
is the (per particle) flattening cost to completely unbend the rows from their preferred curvature, and the rescaled width and height are |  | (16a) |
|  | (16b) |
respectively. 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) |  | (17) |
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 kx = ky. (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
|  | (18) |
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
|  | (19) |
which corresponds to nearly uniformly-curved (
i.e. circular) ribbons characterized by a mid-row curvature
|  | (20) |
A 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)
−1 ≥
κ0−1 and we generally consider the limit of small curvature frustration where
κ0H ≪ 1. To clarify, due to the finite curvature “wide” ribbons, it is not strictly possible to extend
W beyond 2π
κ(wide)c(
h)
−1, as the free ends of the ribbons will meet. However, we expect that in the small frustration limit (
i.e. κ0a ≪ 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−2 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 ê2) or wide/annular (extending along ê1), 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, |  | (21a) |
|  | (21b) |
where σx and σy are rescaled edge binding energies given by |  | (22a) |
|  | (22b) |
These 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 ∼κ0−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 kx = ky 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 kx/ky = 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,33 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/h2 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 kx/ky. Notably, this ratio controls the size of the bend-penetration λx relative to the discrete WJ particle size a. As kx/ky 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 kx/ky 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-fluctuations29 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 kx/ky. 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.29 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 κ0−1 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 κ0−1 ≫ 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 (ê1) binding strength below its critical value (σx < 1), results in tall straight ribbons extending unlimited along the vertical (ê2) 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 ê1 and ê2 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 κ0 is not aligned to one of the two lattice directions (ê1 or ê2), 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 ê2 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.44
 |
| Fig. 7 Assembly along off lattice directions. (a) Definition of the angle ϕ between the long-axis of a ribbon and the ê2 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 ê1 and ê2, 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 kx = ky and correspondingly Bx = By, Yx = Yy 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 κ0a = 0.1 and isotropic elastic constants kx = ky. 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 κ0via the boundary condition of eqn (11), which for isotropic elasticity reduces to (∇θ − κ0
)·
= 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 κ0 and a corresponding elastic cost that is quadratic in preferred row curvature κ0. 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
(and
) have projections along both row and column directions resulting in preferred curvatures κ1 = κ0
sin
ψ and κ2 = κ0
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
|  | (23) |
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
κ1κ2, and that the only terms allowed at quadratic order are proportional to
κ12 and
κ22, 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 κ1 = κ2). It is most useful to consider the implications of this superposition in terms of the infinite ribbon domains given by
|  | (24) |
|  | (25) |
with corresponding mesoscopic curvatures
|  | (26) |
|  | (27) |
These 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
. 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 ê2 (ϕ = 0) or ê1 (ϕ = π/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),
|  | (28) |
where
![[small script l]](https://www.rsc.org/images/entities/char_e146.gif)
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.
44 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]](https://www.rsc.org/images/entities/char_e146.gif)
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 . (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 (
, ϕ) 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
) 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. (ϕ − ψ)2) 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,44 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 particles25,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 kx/ky > 1, which can be achieved by introducing relatively short range binding in the ê1 direction to make rows stiffer to bend while introducing longer range binding in the ê2 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,29 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 ê1 × ê2 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αt2 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”:49 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 {ê1, ê2}, 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 at |  | (29a) |
|  | (29b) |
|  | (29c) |
|  | (29d) |
where α1 and α2 are the taper angles of the left-right (ê1) and top-bottom (ê2) 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) = ma
+ naŷ with the particle frames aligned with the lattice directions {ê1, ê2} = {
, ŷ}, 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 k1 and k2 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) + δtftotal | (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
θ =
ê2·
ŷ, 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
rsite 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) = ma
+ 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 by |  | (32a) |
|  | (32b) |
where we have defined Δnu(m,n) = u(m,n+1) − u(m,n), Δnê(m,n)1 = ê(m,n+1)1 − ê(m,n)1, and 〈ê(m,n)2〉n = (ê(m,n+1)2 + ê(m,n)2)/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 k1 and k2 for the left-right and top-bottom edges respectively, the energy of a sheet is |  | (33) |
The continuum limit can be obtained by taking Δm/a → ∂x, Δn/a → ∂y, and
resulting in |  | (34) |
where we identify the moduli and curvatures as |  | (35) |
and κ1 = (2
tan
α2)/a and κ2 = (2
tan
α1)/a, which represent the preferred bending orthogonal to the ê1 and ê2 directions, respectively. Defining rotations θ relative to the reference state, we have ê1 = cos
θ
+ sin
θŷ and ê2 = −sin
θ
+ cos
θŷ. In terms of the positional and orientational strains defined by | γij = ∂iuj − [Rij(−θ) − δij], | (36a) |
| β = ∇θ + κ2 − κ1ŷ, | (36b) |
where Rij is the 2D rotation matrix and δij is the Kronecker delta, the energy becomes |  | (37) |
For small rotations θ from the reference state of a square lattice aligned with
, ŷ, we may write ê1 ≃
+ θŷ and ê2 ≃ −θ
+ ŷ. The linearized energy functional becomes |  | (38) |
Minimizing with respect to the displacement u and rotations θ, we obtain the Euler–Lagrange equations for force balance | Yx∂x2ux + Yy∂y2ux = −Yy∂yθ | (39a) |
| Yx∂x2uy + Yy∂y2uy = Yx∂xθ | (39b) |
and torque balance | Bx∂x2θ + By∂y2θ = −Yx(∂xuy − θ) + Yy(∂yux + θ) | (40) |
The boundary conditions for the free edges along x = ±W/2 are | (∂xuy − θ)|x=±W/2 = 0 | (41b) |
| (∂xθ + κ2)|x=±W/2 = 0 | (41c) |
and along y = ±H/2 are | (∂yux + θ)|y±H/2 = 0 | (42a) |
| (∂yθ − κ1)|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 κ2 ≠ 0, κ1 = 0 and the other for κ2 = 0, κ1 ≠ 0. We focus on κ2 ≠ 0, κ1 = 0 since the other case can simply be obtained by swapping indices x↔y and curvatures κ2 ↔ −κ1.
To start, it will be useful to briefly summarize the solution for infinitely tall ribbons with κ2 ≠ 0, κ1 = 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
| Bx∂x2θ = −Yx∂xuy + (Yx + Yy)θ | (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
|  | (44b) |
|  | (44c) |
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
| ux(x, y) = u(∞)x(x) + Ux(x, y) | (45a) |
| uy(x, y) = u(∞)y(x) + Uy(x, y) | (45b) |
| θ(x, y) = θ(∞)(x) + Θ(x, y) | (45c) |
Substituting these into
eqn (39a), (39b) and (40), we find that
Ux(
x,
y),
Uy(
x,
y), and
Θ(
x,
y) satisfy the same force and torque balance equations
| Yx∂x2Ux + Yy∂y2Ux = −Yy∂yΘ | (46a) |
| Yx∂x2Uy + Yy∂y2Uy = Yx∂xΘ | (46b) |
| Bx∂x2Θ + By∂y2Θ = −Yx(∂xUy − Θ) + Yy(∂yUx + Θ) | (46c) |
with the boundary conditions
| (∂xUy − Θ)|x=±W/2 = 0 | (47b) |
| (∂yUx + Θ)|y=±H/2 = −θ(∞)(x) | (48a) |
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 write
|  | (49) |
Similarly,
Ux(
x,
y) must also be odd in
x. From the boundary condition
eqn (47a), we can write
|  | (50) |
for some function
Fm(
y). Substituting this into
eqn (46a) and using the orthogonality of the set

, we have that
Fm(
y) satisfies
|  | (51) |
the general solution of which is
|  | (52) |
The boundary condition
eqn (48a) tells us that
|  | (53) |
which gives
|  | (54a) |
|  | (54b) |
Finally, for
Uy(
x,
y), the boundary condition
eqn (48b) allows us to write
|  | (55) |
Substituting into
eqn (46b), we find that
Gn(
x) satisfies
|  | (56) |
By symmetry,
Uy(
x,
y) must be even in
x, and so the general solution for
Gn(
x) can be written as
|  | (57) |
The boundary condition
eqn (47b) tells us that
|  | (58) |
which gives for
n > 0
|  | (59) |
Note that
g0 just gives a global translation, so we pick
g0 = 0 for simplicity. To determine an equation for
Am,n, we take the constructions for
Θ,
Ux, and
Uy and substitute them into
eqn (46c). Using the orthogonality of the sines and cosines, we have for
n even
|  | (60) |
and for
n odd
|  | (61) |
Note that for n odd, gn and cm are sums of Am,n. By inspection of eqn (61), we can conclude that Am,n = 0 for all odd n satisfies eqn (61). Consequently, this means that cm = 0 for all m and gn = 0 for odd n. Thus, we find that Θ(x, y) and Uy(x, y) are even in y and Ux(x, y) is odd in y. We summarize the series solution.
|  | (62a) |
|  | (62b) |
|  | (62c) |
where
sm,
gn, and
Am,n for odd
m and even
n are given by
eqn (54a), (59) and (60), respectively. Unfortunately, due to
sm and
gn depending on infinite sums of
Am,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 κ2 ≠ 0, κ1 = 0. Fortunately, we do not need to repeat the same cumbersome calculation for κ2 = 0, κ1 ≠ 0. All that needs to be done is to take the κ2 ≠ 0, κ1 = 0 solution and simply swap all indices x ↔ y, lengths W ↔ H, and curvatures κ2 ↔ −κ1, and sum up the corresponding displacement and rotation fields. Putting this all together, we have that the general solution for arbitrary κ2,κ1 is
|  | (63a) |
|  | (63b) |
|  | (63c) |
where

,

, and

satisfy similar relationships as
eqn (54a), (59) and (60) but with the swaps
x ↔
y,
W ↔
H, and
κ2 ↔ −
κ1. Note that these solutions are written so that the terms in the first line are proportional to
κ2 while those in the second line are proportional to
κ1.
C.1 Absence of κ2κ1 terms in the elastic energy.
Naturally, we are interested in the elastic energy for general κ2,κ1, which should contain terms proportional to κ22, κ2κ1, and κ12. We show that the elastic energy contains no κ2κ1 terms i.e. the elastic energy can be written as the sum of the individual elastic energies for κ2 ≠ 0, κ1 = 0 and κ2 = 0, κ1 ≠ 0. To see this, let us write the solutions as | θ(x, y; κ2, κ1) = θ(x, y; κ2, 0) + θ(x, y; 0, κ1) | (64a) |
| ux(x, y; κ2, κ1) = ux(x, y; κ2, 0) + ux(x, y; 0, κ1) | (64b) |
| uy(x, y; κ2, κ1) = uy(x, y; κ2, 0) + uy(x, y; 0, κ1) | (64c) |
The parity of each piece is as follows: θ(x, y; κ2, 0) is odd in x and even in y while θ(x, y; 0, κ1) is even in x and odd in y; ux(x, y; κ2, 0) is odd in both x and y while ux(x, y; 0, κ1) is even in both x and y; and finally uy(x, y; κ2, 0) is even in both x and y while uy(x, y; 0, κ1) 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θ + κ2)2 = [∂xθ(κ2, 0) + κ2]2 + [∂xθ(0, κ1)]2 + 2[∂xθ(κ2, 0) + κ2]∂xθ(0, κ1) | (65) |
| (∂yθ − κ1)2 = [∂yθ(κ2, 0)]2 + [∂yθ(0, κ1) − κ1]2 + 2∂yθ(κ2, 0)[∂yθ(0, κ1) − κ1] | (66) |
Since ∂
xθ(
κ2, 0) +
κ2 is even in both
x and
y while ∂
xθ(0,
κ1) is odd in both
x and
y, their product integrates to zero. Similarly, since ∂
yθ(
κ2, 0) is odd in both
x and
y while ∂
yθ(0,
κ1) −
κ1 is even in both
x and
y, their product also integrates to zero. For the stretching energy terms, we have
| (∂xux)2 = [∂xux(κ2, 0)]2 + [∂xux(0, κ1)]2 + 2∂xux(κ2, 0)∂xux(0, κ1) | (67) |
| (∂xuy − θ)2 = [∂xuy(κ2, 0) − θ(κ2, 0)]2 + [∂xuy(0, κ1) − θ(0, κ1)]2 + 2[∂xuy(κ2, 0) − θ(κ2, 0)][∂xuy(0, κ1) − θ(0, κ1)] | (68) |
| (∂yux + θ)2 = [∂yux(κ2, 0) +θ(κ2, 0)]2+ [∂yux(0, κ1) +θ(0, κ1)]2 + 2[∂yux(κ2, 0) + θ(κ2, 0)][∂yux(0, κ1) +θ(0, κ1)] | (69) |
| (∂yuy)2 = [∂yuy(κ2, 0)]2 + [∂yuy(0, κ1)]2 + 2∂yuy(κ2, 0)∂yuy(0, κ1) | (70) |
Following the same procedure, note that ∂
xux(
κ2, 0) is even in
x and odd in
y while ∂
xux(0,
κ1) is odd in
x and even in
y; ∂
xuy(
κ2, 0) −
θ(
κ2, 0) is odd in
x and even in
y while ∂
xuy(0,
κ1) −
θ(0,
κ1) is even in
x and odd in
y; ∂
yux(
κ2, 0) +
θ(
κ2, 0) is odd in
x and even in
y while ∂
yux(0,
κ1) +
θ(0,
κ1) is even in
x and odd in
y; and finally ∂
yuy(
κ2, 0) is even in
x and odd in
y while ∂
yuy(0,
κ1) 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
κ2 ≠ 0,
κ1 = 0 and
κ2 = 0,
κ1 ≠ 0, which are proportional to
κ22 and
κ12, respectively.
Note that this result could also have been obtained from a simple symmetry argument. Flipping the sign of κ2 or κ1 simply creates a subunit that is reflected across ê1 or ê2, 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 κ2κ1 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 Am,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 κ2 ≠ 0, κ1 = 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 ux and uy 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 ux = 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, ux must have some nonzero y dependence, particularly concentrated near the top and bottom edges of the assembly.
• The vertical displacements uy 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 uy 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 uy (from cosh to parabolic or circular) and the leading order dependence of ux on y since ux = 0 in the infinitely tall ribbons. This can be achieved by taking the leading n mode i.e. n = 0 and setting Am,n = 0 for n > 0. The result is the following approximation
|  | (71a) |
|  | (71b) |
|  | (71c) |
where we define
Am =
Am,0. This approximation combined with the conditions
eqn (54a) and (60) for
n = 0 allows us to write in closed-form
|  | (72) |
and
|  | (73) |
Substituting this approximation in the elastic energy functional, we obtain
|  | (74) |
where we rescaled by the flattening energy
flat =
Bxκ22/2.
We can check several limits. For tall assemblies with h → ∞, we find
|  | (75) |
For wide assemblies with
w → ∞ (equivalently, an annulus), we find
|  | (76) |
Finally for small assemblies with
w,
h ≪ 1, we have to leading order
|  | (77) |
We briefly prove that this is symmetric about
h =
w. We use the following identity (see
eqn (75))
|  | (78) |
The excess energy can then be rewritten as
|  | (79) |
Using
|  | (80) |
and
|  | (81) |
we 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. ur ≪ R. The elastic energy between two adjacent rows is simply |  | (82) |
Minimizing this energy with respect to R and δθ, we have |  | (83) |
and |  | (84) |
where Rc is the centerline radius. Note that for N = 1, we have that the centerline radius is Rc = a/(2
tan
α) = κ0−1, which is exactly the preferred radius of a single ring of particles. For N ≫ 1, we have Rc ≃ N2a/(6tanα) = N2/(3κ0). So, the centerline curvature scales goes as κc ≃ κ0/(N2/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
|  | (85) |
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
- G. H. Wannier, Phys. Rev., 1950, 79, 357 CrossRef.
- M. F. Collins and O. A. Petrenko, Can. J. Phys., 1997, 75, 605 CrossRef CAS.
- J. Vannimenus and G. Toulouse, J. Phys. C: Solid State Phys., 1997, 10, L537 CrossRef.
- R. A. Reddy and C. Tschierske, J. Mater. Chem., 2006, 16, 907 RSC.
- H. Takezoe and Y. Takanishi, Jpn. J. Appl. Phys., 2006, 45, 597 CrossRef CAS.
- 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 CrossRef PubMed.
- R. Ghafouri and R. Bruinsma, Phys. Rev. Lett., 2005, 94, 138101 CrossRef PubMed.
- S. Armon, H. Aharoni, M. Moshe and E. Sharon, Soft Matter, 2014, 10, 2733 RSC.
- M. Zhang, D. Grossman, D. Danino and E. Sharon, Nat. Commun., 2019, 10, 3535 CrossRef.
- E. A. Matsumoto, G. P. Alexander and R. D. Kamien, Phys. Rev. Lett., 2009, 103, 257804 CrossRef.
- 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 CrossRef CAS.
- P. Sharma, A. Ward, T. Gibaud, M. F. Hagan and Z. Dogic, Nature, 2014, 513, 77 CrossRef CAS PubMed.
- 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 CrossRef CAS.
- M. S. Turner, R. W. Briehl, F. A. Ferrone and R. Josephs, Phys. Rev. Lett., 2003, 90, 128103 CrossRef CAS PubMed.
- G. M. Grason and R. F. Bruinsma, Phys. Rev. Lett., 2007, 99, 098101 CrossRef PubMed.
- Y. Yang, R. B. Meyer and M. F. Hagan, Phys. Rev. Lett., 2010, 104, 258102 CrossRef PubMed.
- A. I. Brown, L. Kreplak and A. D. Rutenberg, Soft Matter, 2014, 10, 8500 RSC.
- D. M. Hall, I. R. Bruss, J. R. Barone and G. M. Grason, Nat. Mater., 2016, 15, 727 CrossRef CAS PubMed.
- G. Meng, J. Paulose, D. R. Nelson and V. N. Manoharan, Science, 2014, 343, 634 CrossRef CAS.
- W. T. M. Irvine, V. Vitelli and P. M. Chaikin, Nature, 2010, 468, 947 CrossRef CAS PubMed.
- R. E. Guerra, C. P. Kelleher, A. D. Hollingsworth and P. M. Chaikin, Nature, 2018, 554, 346 CrossRef CAS PubMed.
- S. Li, R. Zandi, A. Travasset and G. M. Grason, Phys. Rev. Lett., 2019, 123, 145501 CrossRef CAS.
- M. Lenz and T. A. Witten, Nat. Phys., 2017, 13, 1100 Search PubMed.
- H. Le Roy, M. M. Terzi and M. Lenz, Phys. Rev. X, 2025, 15, 011022 Search PubMed.
- 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 CrossRef CAS.
- 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 CrossRef.
- B. Tyukodi, F. Mohajerani, D. M. Hall, G. M. Grason and M. F. Hagan, ACS Nano, 2022, 16, 9077 CrossRef CAS PubMed.
- N. Tanjeem, D. M. Hall, M. B. Minnis, R. C. Hayward and G. M. Grason, Phys. Rev. Res., 2022, 4, 033035 CrossRef CAS.
- I. R. Spivack, D. M. Hall and G. M. Grason, New. J. Phys., 2022, 24, 063023 CrossRef.
- D. M. Hall, M. J. Stevens and G. M. Grason, Soft Matter, 2023, 19, 858 RSC.
- F. Serafin, J. Lu, N. Kotov, K. Sun and X. Mao, Nat. Commun., 2021, 12, 4925 CrossRef CAS PubMed.
- S. Meiri and E. Efrati, Phys. Rev. E, 2021, 104, 024703 CrossRef.
- M. F. Hagan and G. M. Grason, Rev. Mod. Phys., 2021, 93, 025008 CrossRef CAS PubMed.
- G. M. Grason, J. Chem. Phys., 2016, 145, 110901 CrossRef.
- 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 CrossRef CAS PubMed.
- 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 CrossRef CAS PubMed.
- Y. Suzuki, G. Cardone, D. Restrepo, P. D. Zavattieri, T. S. Baker and F. A. Tezcan, Nature, 2016, 369 CrossRef CAS.
- K. T. Sullivan, R. C. Hayward and G. M. Grason, Phys. Rev. E, 2024, 110, 024602 CrossRef CAS.
- M. Wang and G. M. Grason, Phys. Rev. E, 2024, 109, 014608 CrossRef CAS.
- N. W. Hackney, C. Amey and G. M. Grason, Phys. Rev. X, 2023, 13, 041010 CAS.
- S. Schneider and G. Gompper, Europhys. Lett., 2005, 70, 136 CrossRef CAS.
- F. Cui, S. Marbach, J. A. Zhen, M. Holmes-Cerfon and D. J. Pine, Nat. Commun., 2022, 13, 2304 CrossRef CAS.
- W. B. Rogers and J. C. Crocker, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 15687 CrossRef CAS PubMed.
- C. Rottman and M. Wortis, Phys. Rep., 1984, 59 CrossRef CAS.
- D. Hayakawa, T. E. Videbæk, G. M. Grason and W. B. Rogers, ACS Nano, 2024, 18, 19169–19178 CrossRef CAS PubMed.
- 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 CrossRef CAS.
- 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 CrossRef CAS.
- 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 RSC.
- J. V. Selinger, Annu. Rev. Condens. Matter Phys., 2022, 13, 49–71 CrossRef CAS.
- J. Pollard and G. P. Alexander, New J. Phys., 2021, 23, 063006 CrossRef.
- L. C. B. da Silva and E. Efrati, New J. Phys., 2021, 23, 063016 CrossRef.
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. |
|
This journal is © The Royal Society of Chemistry 2025 |
Click here to see how this site uses Cookies. View our privacy policy here.