Twist renormalization in molecular crystals driven by geometric frustration

Abstract

Symmetry considerations preclude the possibility of twist or continuous helical symmetry in bulk crystalline structures. However, as has been shown nearly a century ago, twisted molecular crystals are ubiquitous and can be formed by about 1/4 of organic substances. Despite its ubiquity, this phenomenon has so far not been satisfactorily explained. In this work we study twisted molecular crystals as geometrically frustrated assemblies. We model the molecular constituents as uniaxially twisted cubes and examine their crystalline assembly. We exploit a renormalization group (RG) approach to follow the growth of the rod-like twisted crystals these constituents produce, inquiring in every step into the evolution of their morphology, response functions and residual energy. The gradual untwisting of the rod-like frustrated crystals predicted by the RG approach is verified experimentally using silicone rubber models of similar geometry. Our theory provides a mechanism for the conveyance of twist across length-scales observed experimentally and reconciles the apparent paradox of a twisted single crystal as a finite size effect.

---

1 Introduction

Nearly a century ago Ferdinand Bernauer made the extraordinary claim that more than 25% of organic molecular crystals can be crystallized into helical forms.¹ Recently these results have been confirmed and elaborated using analytical techniques.^2,3 Moreover, it was demonstrated that twist is observed also for single crystals that display no discernible density of defects.⁴ Our definition of a crystal is thus challenged, as the existence of twist inevitably breaks translation invariance, and therefore cannot be supported by any crystallographic structure. Despite the ubiquity of twisted crystals, the precise mechanism leading to their formation is not fully resolved yet.

N-Benzoylglycine (hippuric acid, HA, shown in Fig. 1) is unique among twisted-crystal-forming componds in that it can form twisted morphologies crystallized from the melt, solution, or vapor phase. When grown by sublimation it crystallizes as thin rectangular rods with twist period (pitch) that increases with the rod thickness. The pitch ranges from a few microns to several hundred microns. Existing models aiming to explain the size dependent shape of HA rely on large concentrations of crystalline defects that cooperate in order to twist and deform an otherwise straight crystalline lattice.³ However, selected area electron diffraction images show sharp spots, SEM images demonstrate well defined facets and AFM images show clear uninterrupted atomic steps, all contributing to the observation that twisted HA is indeed a single crystal and showing no evidence of a high defect density.⁴ Moreover, quantitatively similar twisting behavior is observed when HA is crystalized from the melt⁵ or solution⁶ suggesting that the twisting is a result of material properties rather than due to specific crystallization dynamics.

Fig. 1 SEM image of N-benzoylglycine (hippuric acid, HA) crystals exhibiting size dependent pitch length.

Size dependent behavior, such as that described above, is characteristic of geometrically frustrated assemblies.⁷ Geometrically frustrated systems are associated with two or more mutually contradicting geometric tendencies often giving rise to a local order that cannot be globally extended. Such system have no stress free configuration even in the absence of external constraints, and can be realized, mathematically or physically in many different ways.^8–10 In the spirit of crystallization processes, a frustrated assembly could be formed by the ordered aggregation of incompatible units that must distort in order to fit next to one another in a periodic manner. This distortion in turn leads to residual stresses that may increase as the system grows and can gives rise to a variety of phenomena, such as defect creation,¹¹ selection of high aspect ratio and filamentation,^10,12,13 as well as growth arrest and size limitation.^7,9,14

Intrinsic twist tendency of the constituents of an assembly is in general incompatible with long range order and thus leads to geometric frustration. This incompatibility takes different forms in different systems. In the liquid crystalline blue phase it arises from the geometric inability to sustain a biaxial twist of the director field in a three dimensional volume leading to gradual deviations from the biaxial twist and to finite width domains.¹⁵ A distinct type of frustration arises in twisted bundles of filaments that are allowed to slide past one another.¹⁶ In this case frustration arises from the inability to sustain constant filament spacing given the uniform twist. Here too the frustration may change the morphology of the bundle and lead to growth arrest.

In this work, inspired by the twisted molecular crystals of HA, we study the role geometric frustration assumes in the formation of twisted molecular crystals of macroscopic pitch lengths that continuously vary with crystal size. To achieve this we model the crystalline assembly of twisted deformable cubes. The crystalline order is enforced by prescribing the local topology of the constituents and in particular precluding their relative sliding. We observe that these constraints necessitate some deviations from the desired intrinsic geometric structure and that the degree of deviation depends on system size. As the crystalline assembly grows it gradually untwists, approaching an untwisted conformation in the limit of infinite size. We thus interpret the observed twist in molecular crystals as finite size effect and provide a mechanism for the size dependent conveyance of twist from the molecular to the macro-scale.

2 Results

We construct the model of twisted molecular crystals by assembling chiral continuous building blocks in a simple cubic lattice. The building blocks are cubes that possess a right handed twist along a single axis in their rest configuration, as depicted in Fig. 2A. We align the twist axis of all cubes along the z-direction, and identify adjacent faces of cubes to form a simply connected solid, i.e. the internal faces of the cubes are in full non-slip contact with their immediate neighbors. Physically, this assembly rule is supposed to model adjacent molecules that interact along multiple points such that both their relative positioning and relative orientations have a single preferred value, and deviations from this value increase the elastic energy. In particular, it precludes crystalline defects.

Fig. 2 (A) The basic building block. (B) The lateral assembly of two cubes (2 × 1) necessitating some stress. (C) The longitudinal assembly of two cubes (1 × 2), requiring no stress. (D) 2 × 2 × 2 assembly of cubes. (E) 2 × 2 assembly of rods, which is equivalent to a 2 × 2 × 2 assembly of cubes. In D and E cubes and rods are separated from the main body for demonstration purposes. In all panels the color code corresponds to the elastic energy density obtained by minimizing the global elastic energy of the aggregate with respect to rod-like degrees of freedom of the constituents according to eqn (1).

An assembly grows by the successive addition of new building blocks. Any new building block can connect either across the straight faces which are normal to the twist axis or transversely across the curved surfaces. These two possibilities are not energetically equivalent. As depicted Fig. 2B and C, connecting longitudinally across the twist direction incurs no mechanical cost, while connecting the curved surfaces necessitates some deformation and stress.

Assuming the mechanical energy landscape is convex (i.e. there are no local energy minima, at least not in the vicinity of the global minimum), implies that the conformation of the assembled object is independent of the order in which the constituents came together. Often geometric frustration may lead to non-convex elastic energies. However, for the case considered here, the simple geometry of the assembly and the proximity of the frustrated equilibrium to the unconstrained elastic equilibrium of the disconnected constituents suggest that the resulting elastic energy may indeed be convex. This is also verified experimentally by the robust stability of the unique equilibrium of the frustrated assembly.

We are thus free to recast the assembly in the most convenient form; we choose to consider each solid as composed of columns of cubes stacked along the twist axis and then connected to each other laterally. We term each of these columns the elementary twisted rod. The main advantage of this approach is that each of the elementary twisted rods is stress free before attaching to its neighbors and shares the same twist as the twisted cubes from which it is comprised. For example we can view the cubic eight blocks assembly, as an assembly of four rods each composed of two blocks stacked along their twist axis (see Fig. 2D and E). Following this approach we will model a solid of dimensions w × w × L as comprised of N² twisted rods of square cross-section with side length w₀, length L and preferred twist t₀, such that w = Nw₀.

HA crystals grow as twisted rods that unwind as they grow thicker.^4,5 Very thick rods unwind to the extent that their twist becomes indiscernible. Similarly, in our model one can show that the bulk assembly (obtained by taking the limit N = w/w₀ → ∞) is morphologically trivial and displays no twist (see Appendix E). Thus to explore the origin of twist in molecular crystals one needs to consider them as assemblies of finitely many rods, rather than as a continua. This path to coarse graining is realized using a renormalization approach in which rods hierarchically assemble into thicker rods, as depicted in Fig. 3.

Fig. 3 Illustration of the RG process. At each step 4 identical rods are glued together in full contact, and the resulting structure is relaxed to its elastic energy minimum. In this example the twist renormalizes by a factor of ¼ at every step, as predicted theoretically for the aggregation of rod made from an isotropic and uniform material.

We model each of the elementary twisted rods mechanically as a Cosserat-rod, i.e. a rod that can twist about its mid-line, bend, stretch, shear its cross-section with respect to its mid-line and uniformly dilate its cross-section^17–19 (see Appendix A for more details). Each of these deformation modes is associated with a preferred reference value, denoted with an over-bar and elastic moduli prescribing the energetic cost of deviating from these preferred values. The elastic energy, to leading order, may be brought to the form

where the configuration variables vector is given by Ψ = (v₁, v₂, v₃, u₁, u₂, u₃, σ). v₁ and v₂ describe shear deformations of the rod, v₃ describes elongation in the direction normal to the cross section, u₁ and u₂ describe bending of the rod, u₃ measures twisting of the cross-section, and σ accounts for cross sectional dilation. For the elementary twisted rod parameterized by arc length the reference values of these quantities read . The entries in the coupling matrix G are the elastic moduli of the rod. We assume the length of rods is effectively infinite, and therefore we only consider deformations that are invariant along the rods’ long axes (i.e. we disregard finite length effects).

We start by assembling together four copies of the elementary twisted rod together to form a 2 × 2 aggregate where each of the rods is in full contact with it neighboring rods. The condition of full contact cannot be realized without some deviation from the reference parameters of the individual rods, i.e. cannot hold for all rods in the aggregate simultaneously. Finding the residually stressed configuration with least elastic energy gives the equilibrium configuration the 2 × 2 aggregate rod will adopt. The full contact condition and the symmetry of the aggregate require that the cross-section of the individual rods must remain normal to the z-axis, implying that the constituent rods wind around each other not by bending but rather by shearing. One can show that to leading order longitudinal stretching and cross-sectional dilations are not required, and that there is no linear coupling between shear and twist deformations, (see Appendix A). These lead to a simplified elastic energy per unit length of the aggregate that reads

where G⁰_s, G⁰_t and t₀ are the shear and twist moduli of the elementary twisted rod and its reference twist, respectively. For a square rod of width w₀ made from an isotropic and uniform material of Young's modulus Y and Poisson ratio ν these elastic moduli readHowever, if we consider the elementary rod as a vertical stack of the molecular building blocks, while it may be well described as a Cosserat rod, the above relation between the shear and twist moduli need not hold. The ratio between the twist and shear moduli can be used to define an effective mechanical width , which may differ from the physical width of the rod (noting again that for the isotropic and uniform case w_eff = w).

In the aggregate the full contact condition constrains the values the deformation may assume on the boundary connecting adjacent rods. Denoting the discontinuity of the rod's shear across the connecting face by Δv we obtain (see Appendix)

|Δv₁| = |Δv₂| = w₀u₃. (4)

Minimizing the elastic energy (2), these constraints yield v₁² = v₂² = u₃²w₀²/4 for each of the rods in the assembly. When substituted into (2), we obtain a one dimensional effective elastic energy as a function of u₃ alone, whose minimum is at
where α₀ ≡ (w₀/w_eff,0)² and reads 1 for an isotropic and uniform rod.

Thus, in the absence of further external constraints, the aggregate will adopt the configuration of a twisted straight rod with twist t₁. The resulting aggregate, provided that it still could be considered slender, may again be treated as a Cosserat rod. It is associated with the reference values , and its new elastic moduli matrix G₁ can be calculated directly from the response properties of its constituents.

This provides the basis for a renormalization group approach in which at every step four rods from generation n − 1 are joined to form the nth-generation rod of width w_n = 2ⁿw₀. The first two such generations are illustrated in Fig. 3. Starting with the reference values for the elementary twisted rod cited above we can show that all the reference values remain trivial except for the twist. Knowing the shear and twist moduli at every step we may repeat the minimization of the elastic energy (2) only replacing t₀ by t_n−1. The renormalization flow of the shear and twist moduli reads

This result can be intuitively understood by observing that a pure shear of the four rod aggregate corresponds to a pure and identical shear for each of its constituents. However, a pure twist of the four rod aggregate corresponds to an equal twist of the constituent rods accompanied by an unavoidable shear according to eqn (4). Moreover, eqn (5) suggests , which converges to the physical width w_n as n grows.

At every step the reference twist rescales according to

Note that for an isotropic and uniform material (i.e. α₀ = 1) we have for any step n

For a non isotropic material, initially the untwisting rate can be faster (for α₀ > 1) or slower (for α₀ < 1). However, for large n, α tends to 1, and thus the untwisting rate approaches that of an isotropic material, .

A rod theory is derived by taking a series expansion in the system's width while assuming slenderness, i.e. assuming that the width is the smallest length scale in the system. In the present case of twisted rod assembly, the length of the rods is assumed to be infinite, and the only remaining length scale to compare the width to is the pitch (inversely related to the twist, t). Therefore a rod is slender if wt ≪ 1. When the width of the system is increased the assumption of slenderness may break and thus invalidate the use of rod theory. However, in the case considered here, whenever then doubling the width results in increasing the pitch by more than two for all steps. Moreover, even if and wt will increase in a few initial steps, if the initial w₀t₀ is sufficiently small then the validity of the rod theory may hold. Therefore, the system becomes effectively more slender when it thickens, and the use of rod theory throughout is justified.

In each renormalization step the aggregate rod untwists further, and residual energy per rod, , accumulates according to

As n → ∞ the aggregate rod fully unwinds and the total residual energy per rod converges to the energy required to completely unwind the elementary twisted rod e_∞ ≡ G⁰_tt₀², as can be observed in Fig. 4.

Fig. 4 Residual energy per rod divided by e_∞, the energy of a fully untwisted rod (also the asymptotic residual energy), for different values of α₀. The dashed line represents the constant ‘trivial’ energy required to assemble untwisted rods.

Throughout the above treatment we assumed that the assembly forms a simply connected solid without cuts or holes. More specifically, if identical two dimensional coordinate grids were drawn on each of the faces of the rods, then the assembly rule implies that these coordinate grids must coincide on the common face between two adjacent rods. In general the direction of the third material coordinate (pointing away from the common surface) need not agree between adjacent rods. However, we observe that in the elastic equilibrium these directions also coincide. Framed more mathematically, our assembly rule requires only continuity of the embedding, yet at equilibrium we obtain also continuity of derivatives. Consequently, the configuration adopted by the aggregate rod coincides with the linear profiles of a uniform rod of similar cross-section despite the non-uniform profile of residual stresses in it. This motivates a top-down approach complementary to the renormalization approach employed above. In this top-down approach a large rod of cross-section A is decomposed to N × N rods of smaller cross sectional area A₀ but with similar aspect ratio. This top-down approach allows us to generalize the result obtained from the renormalization approach, namely that t ∼ t₀A₀/A to arbitrary N (not only 2ⁿ), and to non-square elementary rod cross-sections, and even extend to rectangular assemblies of square rods to account for anisotropy in the assembly of the aggregate. See Appendix D.

The model presented here aims to describe how molecular constituents come together to produce twisted crystals of pitch length measuring tens to hundreds of microns. However, as the main results appearing in eqn (5) and (6) are not scale dependent, we may test them using a table-top centimeter sized model. We cast silicone rubber elastic twisted rods and glue them to one another according to the full contact assembly rule described above (see Appendix F for methods used). Fig. 5 shows 1 × 1, 2 × 2, 3 × 3 and 4 × 4 aggregate rods, along with their measured twist, which are shown to agree well with their values predicted by eqn (6).

Fig. 5 Experimentally measured and theoretically predicted twist for 1 × 1, 2 × 2, 3 × 3 and 4 × 4 aggregate rods. The images at the bottom are photographs of the corresponding silicone rubber models. The red dashed line is the theoretical prediction for an isotropic material, interpolated between integer values. Inset shows the rotation measuring technique.

3 Discussion

We consider the problem of assembly of molecular constituents that possess some internal tendency to twist relative to one another. By reformulating the problem as the equivalent lateral assembly of square twisted rods we are able to follow the growth process and predict the gradual straightening of the twisted molecular crystal, as well as the increasing residual elastic energy associated with frustration in the assembly. This allows us to directly relate the moderate twist observed in macroscopic molecular crystals with the potentially nanometric pitch of the individual molecular building blocks, and to the elastic properties of the molecular assembly.

It is natural to seek a continuum approach to the problem. Starting from a field theoretic approach one could write a theory that includes both elastic contributions, and liquid crystalline contributions (which include higher order spatial derivatives) that together reproduce the long wavelength behavior of the system.¹⁴ However, the origin of twist in this case is less transparent. Obtaining these field equations from first principles for molecular crystals is a challenge.

A step towards obtaining such a continuum theory from first principle is to identify the geometric charge associated with the frustrated assembly of twisted rods. However, the homogenization of the assembly presented here, obtained by considering a domain of constant size W × W × L composed of twisted square rods of dimensions w × w × L and of twist t₀, and taking the limit w → 0 with the parameters W, L and t₀ kept constant, results in a trivial Euclidean bulk with no residual elastic energy (see Appendix E for proof). Thus we cannot associate the frustration caused by a local twist with a simple Riemannian charge as one could do for the case of uniformly bent constituents²⁰ or for incompatible surfaces.²¹

Geometrically frustrated systems often display a super-extensive elastic energy, i.e. that the energy per unit volume grows as the assembly becomes larger. For twisted rod aggregates we observe super extensive elastic energy for small aggregates (as observed in Fig. 4). However, for large assemblies the elastic energy per unit volume converges to a constant, rendering the elastic energy asymptotically extensive. Geometrically frustrated systems that possess a strongly super-extensive elastic energy exhibit a preference to form filamentous structures. In the present case, such a tendency may only exist in the early stages of the growth. Thus the filamentous profiles observed experimentally are likely not due to frustration alone but rather a signature of the early dynamics of the crystallization or of some other growth limiting factors.

Considering approximately isotropic rods in which w₀ ∼ w_eff and α₀ ∼ 1, the twist decays inversely proportional to the cross-sectional area of the aggregate from the very first renormalization step. However, if w₀ ≪ w_eff then α₀ ≪ 1 and the twist decays slower in the initial stages. As the aggregate grows, α_n → 1 and the twist will again decay with the inverse of the cross-sectional area of the aggregate. As the molecular crystals observed in experiments are of micron size, orders of magnitude larger than the size of their molecular constituents they should be interpreted as large aggregates in which the value of α_n is close to unity, and at every width doubling one expects the twist to diminish by a factor of four. It is important to note that in this late regime it is difficult to measure t₀ and α₀ separately. Naively extracting t₀ while assuming α_n = 1 for all n may lead to strong overestimation of the initial twist t₀ for cases where α₀ is small.

While the initial unwinding rate depends on the material parameter through α₀, in all cases we considered, geometric frustration provides an effective means for conveying twist from the molecular scale to the scale of the aggregate. This results in size dependent pitch unwinding the crystal as it grows, and reconciles the simultaneous appearance of continuous twist and crystalline order as a finite size effect. Translational symmetry may not lead to the optimal arrangement of constituents for small crystals, but rather arise as an imperfect compromise in view of stringent geometric constraints as the crystal grows in size.

Conflicts of interest

There are no conflicts to declare.

Appendices

A Cosserat rod energy for the isotropic case

A Cosserat rod consists of a parametrized curve r⃑(s) and a moving Lagrangian orthonormal material frame of vectors d̂₁(s), d̂₂(s), d̂₃(s), such that the three dimensional embedding of the rod is given by

R⃑(s,x,y) = r⃑(s) + σ(xd̂₁(s) + yd̂₂(s)). (7)

Here (x,y) are the Lagrangian cross-section coordinates, not be confused with the absolute directions in space, and σ is the isotropic cross-sectional dilation. The evolution of the frame vectors with respect to the curve parameter s is given by the set of ODE's (using the Einstein summation rule over three Euclidean coordinates indexed by i, j, k):

∂_sd̂_i = −ε_ijku_jd̂_k, (8)

where ε_ijk is the Levi-Civita symbol and {u_j} are the components of the Darboux vector u⃑ = {u₁, u₂, u₃}. The elements u_j of the Darboux vector can in general be functions of s. Considering long rods we will limit our theory to only constant Darboux entries. The curve r⃑(s) is defined by the projections v_i of its tangent on the frame vectors:

∂_sr⃑(s) = v_id̂_i, (9)

where we similarly limit ourselves to consider only constants for the projections v_i. If then s is the arc length parameter for the curve r⃑(s).

The elements u_i and v_i have clear geometrical interpretations; u₁ and u₂ correspond to bending while u₃ corresponds to twisting. v₁ and v₂ measure shearing of the cross section with respect to the tangent curve, while v₃ accounts for elongation in the direction normal to the cross section, (see Fig. 6).

Fig. 6 The fundamental deformations of a Cosserat rod with uniform parameters.

We proceed to derive the energy of a Cosserat rod. The energy is a function of the strain, which is proportional to the difference between the metric of the rod's current configuration and its reference configuration metric. We use an overbar, ^–, to mark the parameters of the reference configuration. The reference dilation is set to 1 (without loss of generality).

We calculate the metric tensors g_ij, ḡ_ij using (7):

We then use the metric description of elasticity described in ref. 8 and 12 to calculate the energy, per unit length:where the lowering and raising of the indices is performed using the reference metric ḡ. Y is Young's modulus, and ν is Poisson's ratio. Ω is the domain of the cross section, , thus defining w as the (Lagrangian) width of the cross section.

We assume that the rods’ reference state is close to that of a straight simple rod (with vanishing shear, bend, twist and stretch), and that strains are small, implying that the parameters σ, u_i and v_j are close to their reference values. We do not restrict the relative magnitude of the different strain components, and simply expand to second order in the deviations from trivial geometry and in the strains.

We then take a 2nd order Taylor expansion of the energy in all of the reference variables, and assign a scaling exponent to each of the reference variables (e.g.), and retain the lowest significant orders.

The second assumption is that of small strains. To do this we expand up to second order the configuration variables about the references (e.g. we expand u₃ about ū₃). We assign to each strain variable (i.e. terms like u₃–ū₃) a different scaling exponent, and any remaining reference configuration variables are assigned the same scaling exponents as in the first assumption. We retain only the lowest significant orders in each of the independent scaling exponents and derive the energy of a single rod:

Alternatively, the energy can be cast in the quadratic form as in eqn (1), with the elastic moduli matrix G given by:

Furthermore, as shown in Appendix D, in an aggregate rod with an elementary unit that has only non-zero reference twist (i.e. ū₃ ≠ 0, ū₁ = ū₂ = 0), the only degrees of freedom (d.o.f) that evolve with system size or composition are the shear and the twist. We can therefore use a simplified energy, as used in the main body (2):

B Full contact conditions on adjacent rods

In this section we express the ‘full contact’ assembly rule through constraints on the parameters that adjacent rods may assume. We consider a pair of adjacent rods of width w, and find the constraints by requiring continuity of the embeddings of the two rods on their common boundary (that the parametric curves from both rods coincide).

We first consider an x border, i.e. a left–right border, corresponding to the coordinate value for the left rod, and to for the right rod. The continuity condition is then:

We have so far not restricted the rod parameters σ, v_i, u_i. They are, at this stage, allowed to all differ in both rods. The above equalities need to hold for all s and all y. We can therefore examine them element-wise as a polynomials in y and match the different coefficients. Comparing the y¹ terms we obtain

σ_Ld̂_2L(s) = σ_Rd̂_2R(s).

As the d̂_i frame vectors are of unit length, we must have that σ_L = σ_R ≡ σ. We are thus left with

d̂_2L(s,u_1L,u_2L,u_3L) = d̂_2R(s,u_1R,u_2R,u_3R).

i.e. the d̂₂ frame vector is shared by two adjacent rods. The two remaining frame vectors d̂₁ and d̂₃ have, in principle, a degree of freedom of a constant angle of rotation θ about the shared d̂₂, such that d̂_1L·d̂_1R ≡ cos(θ). In the following passages we will show that when two rods with shear and bend free reference configurations are joined in full contact, the equilibrium configuration has θ = 0, such that all the frame vectors between the two rods are shared. We examine the first derivative of d̂₂ in the left and right rods:Taking the dot product with d̂_1L and d̂_3L, respectively we obtain:Taking the second derivative of d̂₂, and examining the part perpendicular to d̂₂, we obtain:

u_3Lu_2Ld̂_3L + u_1Lu_2Ld̂_1L = u_3Ru_2Rd̂_3R + u_1Ru_2Rd̂_1R. (17)

Taking the dot product with d̂_1L we findwhere the last equality is obtained using (16). Thus we express the u_i parameters of the left rod in term of that of the right. We proceed to do the same for the v_i parameters. We examine the y⁰ term of (14). Taking its s derivative we obtain:Taking the dot product with d̂₂, d̂_1L and d̂_3L respectively, we obtain:We use an energy functional for a rod of width w, which is a generalization of the isotropic energy (12) found in Appendix A:Note that here, for ease of calculation, we specifically retain the dimensional part of the elastic moduli, in contrast to the energy functional used in the main text. We set v̄₁ = v̄₂ = 0, ū₁ = ū₂ = 0 and v̄₃ = 1 to account for rods that have a shear and bend free reference configuration. We express the combined energy of both rods in terms of the parameters of the right rod, using (16), (18) and (20). However, due to some non-linear couplings in (16) and (20), the combined resulting energy of the two rods has terms of higher order than the original quadratic form of each of the single rods. Retaining terms up to second order in all variables, we attempt to find local minima in the energy landscape by equating the gradient of the combined energy to zero. For the ∂_uiR and ∂_viR terms in the gradient, an exact solution exists in terms of θ and σ. The ∂_θ term can however only be solved numerically, and it can be shown that for a wide applicable range of system parameters the solution is θ = 0 (θ = π is also a solution but has higher global energy). After setting θ = 0 an exact solution for the ∂_σ term is also obtained.

We have thus found that when the θ degree of freedom is available to an assembly of two rods, the minimal energy configuration has θ = 0. Assuming the energy landscape is convex about its minima, we set θ = 0 permanently for all further calculations. This implies that the frame vectors d̂_i are shared and that all u_i are also shared. Moreover the shear and elongation parameters along an x border are related by:

Repeating the calculation, for a y border (using labels S for ‘south’ and and N for ‘north’, instead of left and right) we find that

C Adjacent rods in full contact can be described as a single rod

Two linear functions on adjacent domains that are continuous across the domains’ boundary will not necessarily form a linear function on the union of the domains. For example the functionis continuous at x = 1, linear in each of the segments but is not linear in 0 ≤ x ≤ 2. If first derivatives are shown to be also continuous across the domains, the resulting function is indeed linear over the union of domains.

Within the Cosserat rod theory the deformations in the cross-section of the rod are assumed linear. Thus, when joining different Cosserat rods, one may in general obtain a deformation profile that is not linear in the cross-section of the union of the rods. In this appendix we show that while the full contact condition may seem to require only continuity of the embedding, it in fact enforces continuity of derivatives as well.

Using the result of Appendix B we can show that

∂_xR⃑_L(s, x, y)|_x=w/2 = d⃑₂ = ∂_xR⃑_R(s, x, y)|_x=−w/2.

Consequently, when adjacent rods are joined, the deformation profile along their joint cross-section is linear, and one may describe them as a single rod without needing to invoke any further assumptions or approximations.

Thus each of the constituents in a rod aggregate share the bend, twist and dilation with those of the aggregate rod and inherit their shear and elongation from the rod aggregate by appropriately translating the mid-curve by a cross-sectional vector (α₁,α₂) ≡ [small alpha, Greek, vector], and examining the behavior of the mid-curve's tangent vector.

Explicitly we may write:where ṽ_i are the shear and elongation parameters of a constituent rod whose center is at .

D Top-down approach

Based on the finding that rods in a full-contact can be equivalently described as a single rod (see Appendix C), we provide an alternative method to calculate the full configuration of a rod aggregate (reference parameters, elastic moduli and residual energy) in a single step. Moreover, with this method we can generalize the description to rods of a general rectangular cross section, subdivided into any N × M square rods (compared to the 2^N× 2^N square rods assembly accessible through the renormalization).

We consider an N × M rod assembly, of overall width W and thickness T. The width of the constituent square rods that make up the generally rectangular aggregate is . Therefore, we can express . The displacement of each constituent rod mid-curve from the center of the aggregate rod is given by

We use the energy functional from (21), for a rod of width w, which is a generalization of the isotropic energy (12) found in Appendix A. Note again that for ease of calculation we retain the dimensional part of the moduli; for example we write G_tw⁴ instead of simply noting the twist modulus G_t as we do in (2). We proceed to express the rod energy of each of the constituent rods as a function of the parameters of the aggregate rod. We do this by replacing the v_i parameters in each rod according to (25), using the displacements defined in (27). Additionally, we replace the width of each constituent rod , N being the number of times the aggregate is split across the width of the aggregate. Thus the energy of a single constituent rod is E_rod[n,m], and the overall aggregate energy is the sum of all the constituent rods:

We use (26) to express the strain variables on the constituent rods in terms of the strain variables of the aggregate. This relation couples the cross-sectional dilation σ to bend and twist deformations non-linearly, leading to high order terms that are outside of the rod-like approximation. Omitting these we obtainwhere the residual energy E_residual depends on reference parameters only and is given by:

(29) had the same quadratic form as that of a single rod, (as in (21) and (12)), and as such its minimizers (or new references) are found trivially. Note that only the reference bend and reference twist evolve non-trivially with N. For the particular case of vanishing reference bends of the constituent rods, we obtain that the rod assembly has only a non-trivial twist. This twist, however, contributes to both twist and shear of the constituent rods through (26), justifying the form of the simplified rod-energy (2).

Thus the top-down approach is an alternative and more general way to calculate the energy of a rod aggregate for any rectangular assembly. Setting we recover the results of the bottom-up approach presented in the main text (recalling that N² = 4ⁿ, where N is the number of top-down splits and n is the number of RG steps).

Taking N → ∞ while keeping W and T constant is akin to a homogenization procedure, whereby a fixed volume is subdivided into elementary rods with vanishing size. In this limit we find that the residual energy (30) vanishes, as do the twist and bend references. These results can be seen also as consequences of the triviality of the Riemannian limit of the system, discussed in the next appendix.

E Trivial homogenization

In this appendix we prove that the Riemannian homogenization of the twisted rod assembly is trivial. This is in contrast with other frustrated assembly processes where the homogenization is not trivial and its properties dominate the system's behavior. One of the simplest cases is that of an assembly of rods of constant radius of curvature, ρ and width δ connected to one another.²² The discontinuity in geodesic curvature across the curve connecting two adjacent stripes renders the metrics of the assembly at any finite δ non-differentiable. However, the limit of δ → 0 while keeping ρ constant results in a smooth geometry of constant Gaussian curvature K = −1/ρ², as has been demonstrated in ref. 22 and can be inferred from the vanishing splay and constant bend compatibility condition for two dimensional liquid crystals.²⁰

Examining the case of twisted rod assembly leads us to consider joining N × N twisted straight rods of twist t and width w, such that the total width of the assembly's cross section is W = Nw. The homogenization of the rod assembly corresponds to the limit w → 0 while keeping W and t constant.

The metric of a single twisted straight rod with respect to the cross sectional Lagrangean coordinates x and y is given by

We define ḡ_w as the reference metric obtained by joining N × N such rods. The x and y coordinates are continuous between rods. The resulting metric satisfies

ḡ_w(x,y) = ḡ_w(x + w,y) = ḡ_w(x,y + w).

The variation in the diagonal elements of ḡ_w is of order w² whereas the variation of the off diagonal elements is of order w. The diagonal elements of ḡ_w are continuous across the connecting faces but their first derivatives are discontinuous, much like the case of ribbons of constant bend. However, some of the off diagonal terms are themselves discontinuous across connecting faces.

For simplicity we consider the case of a vanishing Poisson ratio for which, up to constants, the elastic energy of a metric g with respect to ḡ_w is given by

where the L² norm with respect to the reference metric ḡ_w reads

The inverse metric for a straight twisted rod reads

Note also that |g| = 1 and the volume element . We set the limiting energy functional to read
Noting that ḡ_w = I + [scr O, script letter O](w) and ḡ_w⁻¹ = I + [scr O, script letter O](w) we obtain

E_w = E₀ + [scr O, script letter O](w). (31)

The difference E_w–E₀ is a second order polynomial in the components of g with coefficients that depend on ḡ_w which in turn is close to the identity matrix. It is thus straightforward to prove that E₀ is indeed the Γ-limit of E_w as w tends to 0. The existence of a recovery sequence is immediate by examining the constant sequence
It remains to show that for every sequence of metrics g_n → g and w_n → 0 as n → ∞ we have
For every metric g such that E₀ < ∞ it is straightforward to show that |E_w(g) − E₀(g)|/w < C for all w, for some 0 < C that may depend on E₀ but not on w. Therefore .

Thus for bounded Euclidean metrics g

Note that not only the minimizer of the limiting functional is the identity, but also that the elastic energy in the limit vanishes. Thus all residual stress is eliminated by considering the homogenization, and finite twist pitches are a finite size effect.

F Methods

The silicone rubber twisted rod assemblies were created by casting liquid silicone rubber (Smooth-On, Inc. Mold-Max™ NV14) in 3D printed molds (Formlabs, Inc, Form 2 printer, using Grey resin). The individual elementary rods were cast as twisted straight rods of pitch length of 16 cm and square cross-section of width 0.8 cm. White and grey pigmented cured rods were clamped into a straight conformation to align their faces and then were glued together using red-pigmented uncured resin (of the same material as the rod) resulting in continuous rod aggregates with no discernible variation in material properties. After curing, the aggregate rods were photographed from below in their unconstrained state. Image edges were detected via a Roberts cross operator (MATLAB ‘edge’ function) and the twist was inferred by identifying the feature orientation using a Radon transform.

Acknowledgements

This work was supported by the National Science Foundation under Award DMR-1608374. It was supported partially by the Materials Research Science and Engineering Center (MRSEC) program of the National Science Foundation (NSF) under Award DMR-1420073. H. A. was supported by NSF grant no. DMR-1262047. E. E. is the incumbent of the Ernst and Kaethe Ascher career development chair and thanks the Ascher foundation for their support, and also acknowledges the support by the Alon fellowship. This work was funded by NSF/DMR-BSF grant 2015670, the Minerva foundation grant number 712273, and by the ISF grant number 1479/16.

Notes and references

1. F. Bernauer, Gedrillte Kristalle, Borntraeger, 1929.
2. A. Shtukenberg, E. Gunn, M. Gazzano, J. Freudenthal, E. Camp, R. Sours, E. Rosseeva and B. Kahr, ChemPhysChem, 2011, 12, 1558–1571.
3. A. G. Shtukenberg, Y. O. Punin, A. Gujral and B. Kahr, Angew. Chem., Int. Ed., 2014, 53, 672–699.
4. A. G. Shtukenberg, A. Gujral, E. Rosseeva, X. Cui and B. Kahr, CrystEngComm, 2015, 17, 8817–8824.
5. A. G. Shtukenberg, J. Freudenthal and B. Kahr, J. Am. Chem. Soc., 2010, 132, 9341–9349.
6. P. Calvert and D. Uhlmann, J. Polym. Sci., Part B: Polym. Phys., 1973, 11, 457–465.
7. G. M. Grason, J. Chem. Phys., 2016, 145, 110901.
8. E. Efrati, E. Sharon and R. Kupferman, J. Mech. Phys. Solids, 2009, 57, 762–775.
9. S. Armon, H. Aharoni, M. Moshe and E. Sharon, Soft Matter, 2014, 10, 2705–2908.
10. G. Meng, J. Paulose, D. R. Nelson and V. N. Manoharan, Science, 2014, 343, 634–637.
11. W. T. M. Irvine, V. Vitelli and P. M. Chaikin, Nature, 2010, 468, 947–951.
12. E. Efrati, E. Sharon and R. Kupferman, Soft Matter, 2013, 9, 8187–8197.
13. S. Schneider and G. Gompper, Europhys. Lett., 2005, 70, 136–142.
14. G. M. Grason, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 041919.
15. P. G. de Gennes and J. Prost, The physics of liquid crystals, Clarendon Press, 1993, p. 324.
16. D. M. Hall, I. R. Bruss, J. R. Barone and G. M. Grason, Nat. Mater., 2016, 15, 727–732.
17. E. Cosserat and F. Cosserat, Bull. Am. Math. Soc., 1913, 19, 242–246.
18. N. Chouaieb and J. H. Maddocks, J. Elasticity, 2004, 77, 221–247.
19. S. Antman, Nonlinear Problems of Elasticity, Springer, 2nd edn, 2005, p. 845.
20. I. Niv and E. Efrati, Soft Matter, 2018, 424–431.
21. S. Armon, E. Efrati, R. Kupferman and E. Sharon, Science, 2011, 333, 1726–1730.
22. D. W. Henderson and D. Taimina, Math. Intell., 2001, 23, 17–28.

---