{'fraction_non_alphanumeric': 0.0379250004088508, 'fraction_numerical': 0.012739790995469933, 'mean_word_length': 4.754564276303407, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
How geometric frustration shapes twisted fibers, inside and out: Competing morphologies of chiral filament assembly

Chirality frustrates and shapes the assembly of flexible filaments in rope-like, twisted bundles and fibers by introducing gradients of both filament shape (i.e. curvature) and packing throughout the structure. Previous models of chiral filament bundle formation have shown that this frustration gives rise to several distinct morphological responses, including self-limiting bundle widths, anisotropic domain (tape-like) formation and topological defects in the lateral inter-filament order. In this paper, we employ a combination of continuum elasticity theory and discrete filament bundle simulations to explore how these distinct morphological responses compete in the broader phase diagram of chiral filament assembly. We show that the most generic model of bundle formation exhibits at least four classes of equilibrium structure -- finite-width, twisted bundles with isotropic and anisotropic shapes, with and without topological defects, as well as bulk phases of untwisted, columnar assembly (i.e."frustration escape"). These competing equilibrium morphologies are selected by only a relatively small number of parameters describing filament assembly: bundle surface energy, preferred chiral twist and stiffness of chiral filament interactions, and mechanical stiffness of filaments and their lateral interactions. Discrete filament bundle simulations test and verify continuum theory predictions for dependence of bundle structure (shape, size and packing defects of 2D cross section) on these key parameters.


I. INTRODUCTION
Assemblies of chiral building blocks often inherit the lack of mirror symmetry at the mesoscale from the microscopic scale of those subunits [1]. This effect is well known in the context of chiral mesophases of liquid crystals [2,3], in which the lack of inversion symmetry of rodlike molecules leads to self-organized, twisted textures. One such state is the cholesteric phase, where the molecular axis adopts a handed twist throughout the sample. In a second class of systems, membranes formed by chiral molecules (e.g. surfactants), the local 2D packing of chiral elements ultimately leads to mesoscopically handed morphologies, including helicoidal and spiral ribbons [4][5][6]. In biological assemblies, where protein building blocks are universally chiral at the molecular scale, the effect of chirality to template morphology at a super-molecular length scale is well known [7][8][9].
In this paper, we consider the influence of chirality to shape rope-like assemblies of 1D filamentous elements, a common materials architecture. Chiral bundles and fibers form in a diverse range of systems, from supramolecular polymers, organogels and columnar fibers [10][11][12][13][14], to protein fibers, both functional (collagen [15,16], cellulose [17], fibrin [18]) and pathological (amyloids [19], sickle hemoglobin [20]). While these systems are clearly diverse in terms of both their structure and interactions at the molecular scale, they share a common geometrical template, with filaments or columnar molecular stacks winding helically along the long axis of the bundle or fiber (see Fig. 1A). This basic geometry leads to generic, yet non-trivial, consequences for the resulting structures and thermodynamics of assembly, owing to the geometrical frustration of long-range order by chiral tex-tures in self-organized bundles.
Geometric frustration refers to the incompatibility of preferred patterns of local order and global constraints, which prevents the propagation of distortion-free order in assemblies [21][22][23]. In rope-like twisted bundles, frustration arises from two geometric mechanisms. The first, orientational frustration, results from the tendency of chirality to promote non-parallel textures of filament backbone orientation throughout the assembly [24]. Lateral gradients of filament bending introduced by chiral ordering give rise to important size-dependent assembly costs that can establish the equilibrium (self-limiting) finite size of bundles [25][26][27][28]. The second mechanism is metric frustration, where the non-parallel, twisted textures of filaments favored by chirality make uniform inter-filament spacing impossible [29]. It has been shown that this metric frustration in twisted bundles can be mapped onto the metric frustration of packing elements on a positively-curved 2D surface [30,31]. Specifically, the distance of closest approach between filaments in a twisted bundle corresponds one-to-one to the geodesic distance between equivalent points on the "dual surface", whose curvature radius is proportional to the helical pitch of bundle twist (see Fig. 1B). This "dual problem" has received broad interest in diverse contexts ranging from the long-standing Thomson problem [32] to the structure of particle-assemblies at liquid droplet interfaces [33]. Previous studies have shown that metric frustration in both bundles and curved crystals leads to intrinsic gradients of stress in the assembly that in turn may trigger a rich set of morphological responses [23], including self-limiting domains [25,26], defective bulk order [32,[34][35][36][37] and anisotropic surface shapes [38][39][40][41].
In this paper, we explore the principles that se- Competing morphologies of self-assembled chiral filaments. Compression at the perimeter is due to geometric frustration of twisted textures. We show equivalent bundle cross section (A) and curved crystal structure (B) using the previously derived mapping. At the top, elements are shaded according to local strain of inter filament distance ∆, relative to ideal separation d. In (C), unwinding preferred chiral twist (or equivalently decreasing spherical curvature) allows a larger assembly with reduced compression. In (D), excess free surface of anistropic, tape-like bundles leads to compression that grows only with thinner cross-section dimension. In (E), the introduction of a 5-fold defect in the inter-filament order relaxes hoop compression.
lect among the multiple competing morphologies stabilized by geometric frustration in chiral filament bundles. These distinct morphological responses are illustrated in Fig. 1, associated with bundle geometry along the bundle (Fig. 1C), bundle surface geometry (Fig.  1D) and symmetries of inter-filament order in the bulk (Fig. 1E). While chiral filament interactions select a preferred amount of bundle twist, the cost of both intraand inter-filament strains may be relaxed via an untwisting of skewed filament packing (Fig. 1C), which corresponds to a geometric "flattening" of the dual metric surface as lateral bundle area grows and allows bundles to reach a larger diameter than possible for fixed twist. Indeed, when chiral interactions are sufficiently "soft', bundle assemblies may "escape" frustration by untwisting completely, and consequently grow to unbounded domains. Alternatively, a second mechanism to relax frustration cost is possible if surface energies (associated with fewer cohesive bonds at the bundle exterior) are sufficiently small compared to the cost to deform interfilament bonds, so that the bundle cross section becomes highly anisotropic as shown in the tape-like morphology in Fig. 1D. Analogous to the elastic instability for growing crystalline domains on spherical surfaces described below, this morphology is driven by a tendency to lower the overall growth of inter-filament strains with increasing filament number in the bundle. Finally, Fig. 1E shows that for sufficiently high overall twist in the assembly, it becomes advantageous to adopt (5-fold) disclination defects in the 2D bundle cross section to lower the cost of geometric frustration of inter-filament spacing, a phenomenon well-studied for the assembly of crystalline arrays of positively curved surfaces.
Several previous studies [24,26,30,31,36,37,41] have investigated how each of these "morphological responses" may effectively lower the cost of geometric frustration imposed by chiral order in bundles. Here we address how these distinct responses compete in the broader phase space of chiral filament assembly. Specifically, we use a combination of continuum elasticity theory and discrete-filament bundle simulations to show that equilibrium morphology is selected by a combination of key thermodynamic (bundle surface energy), geometric (preferred chiral twist) and mechanical (stiffness of filaments and their interactions) properties of filament assemblies. In Sec. II, we present a generic continuum model of defect-free twisted bundle assemblies, which allows us to compare the structure and thermodynamics of finite width bundle assemblies (possessing both isotropic and anisotropic cross sections) to the cost of unwinding chiral interactions to allow bulk order to grow without frustration. We show that the possibility of "escaping frustration" is governed by the ratio of stiffness of chiral filament interactions to the combined surface energy and mechanical costs of forming finite width twisted domains. Above a critical chiral interaction stiffness, bundles exhibit both distinct isotropic (cylinder) and anisotropic (tape-like) morphologies, while below this critical stiffness, finitediameter bundles are favored only within a narrow range of weakly-cohesive surface energies, above which equilibrium bundle growth causes assemblies to untwist to the parallel, unfrustrated state. In Sec. III, we present discrete filament bundle simulations that compare the relative stability of bundles with anisotropic lateral dimensions to bundles which possess topological defects in their interior order. We find that the optimal morphological response (boundary shape or bulk defects) is determined by the filament number N , the bundle twist and the ductility of inter filament bonds. For small-N and ductile inter-filament interactions, bundles pass directly from defect-free cylinders (i.e. isotropic boundary shape) to defective cylinders with increasing twist, while for large-N and brittle interactions, the defective cylinder phase is preempted by an intermediate phase of defect-free anisotropic (tape-like) bundles. We argue that the stability window of anisotropic tapes can be modeled via continuum theory estimates for the relative costs of straining inter-filament bonds (by geometric frustration) compared to the surface energy cost of introducing excess (anisotropic) boundary. We conclude with a discussion of these results in the context of outstanding questions regarding the role of geometric frustration in shaping the assembly of chiral filament bundles.

II. ESCAPING FRUSTRATION TO BULK ASSEMBLY: CONTINUUM MODEL OF MORPHOLOGY SELECTION
In this section, we present a continuum elastic theory for the formation of twisted bundles of chiral filaments that relates the free energy of assembly to bundle morphology, specifically the equilibrium twist and lateral dimensions of the cross section. We begin by introducing the model and its application to the thermodynamics of assembly with isotropic versus anisotropic bundle cross sections, and then analyze the possibility of bundle assemblies to "escape frustration" by relaxing interfilament twist with increasing size.

A. Model
We describe the assembly of twisted bundles as finite domains of filaments with 2D columnar order and preference for chiral textures of filament backbones [37]. A columnar bundle is described by the shape of its 2D cross section normal to the mean-orientation of filaments 2π/Ω r (ẑ axis), a local 2D displacement u(x) (⊥ toẑ) of local filament positions away from a parallel and equallyspaced (e.g. hexagonal) reference configuration, and the unit vector orientation field of filament backbones t(x) = (ẑ + ∂ z u)/ 1 + |∂ z u| 2 . Given the shape and geometry of the bundle (described by u(x) and t(x)) the following free energy describes bundle formation Here, represents the bulk and surface terms for a bundle of length L, cross sectional area and perimeter, A and P , respectively, and areal filament density ρ 0 . The net per filament gain of inter-filament cohesion (minus entropic loss) to add filaments to a bundle is −µ, while fewer cohesive interactions at the lateral surface of the bundle leads to a positive surface energy Σ [48]. The elastic energy has the form, where the first term represents the cost for intra-filament, bending deformation, with B and κ(x) = |∂ z t| the filament stiffness and curvature, respectively. The second term describes the continuum elastic cost of deforming lateral inter-filament order (spacing and symmetry) where σ ij and u ij are the respective in-plane stress and strain tensors (with components in the plane of lattice order, ⊥ toẑ). We assume a isotropic, linear elasticity appropriate for hexagonal inter-filament order u ij = Y −1 (1 + ν)σ ij − νδ ij σ kk where Y and ν are the respective Young's modulus and Poisson's ratio of the 2D filament array. The geometric coupling between positional and orientational degrees of freedom is built into the form of the non-linear columnar strain [42], The non-linear contribution from tilt is necessary to preserve rotational invariance of the theory, and more intuitively, represents the fact the distance of closest approach between filaments will be reduced when filaments tilt relative toẑ (i.e. t i = 0) without changing their positions in that plane (i.e. ∂ i u j = 0). The final contribution to F derives from orientation dependent interactions between chiral, rod-like elements, and hence, has the form of the twist-elastic contribution to the Frank energy of nematic phases [43], where K is the twist elastic constant (often written as K 22 ) and q 0 = 0 parameterizes the drive of chiral filaments to adopt a twisted texture with t · (∇ × t) = −q 0 [49].
In simple chiral nematics (which lack positional order) q 0 = 0 leads to a uniaxial cholesteric ground-state texture, with nematic director rotating along a direction orthogonal to the nematic axis with a pitch 2π/q 0 . Due to the 2D positional order, uniaxial cholesteric textures require constant longitudinal gradients of in-plane shears [44], and hence, are suppressed as L → ∞. Hence, for finite chiral columnar domains such as bundles, chirality favors an alternative "double twist" texture [45], in which filament orientations wind helically around a central axis, and filament positions are simply rigidlyrotated (shear-free to linear order) from one transverse section to the next along the bundle. Parameterizing this texture by the twist Ω of filaments around the axis, the backbone texture is where (r, φ) are polar coordinates with respect to the twist axis. In this geometry, 1) each filament adopts a helical shape and possesses a single pitch P = 2π/Ω; and 2) neighbor filaments spaced δr along the radial direction adopt a nearly constant inter-filament skew angle δθ Ωδr (see Fig. 2). Because t · (∇ × t) 2Ω ∝ δθ for this texture, the preference for double-twist in the bundle can therefore be linked to microscopic interactions between chiral filaments which favor a fixed interfilament skew between neighbors δθ 0 ∝ q 0 and generate an inter-backbone torque ∝ K(δθ − δθ 0 ) for deviations from this preferred local packing for which F chiral /V 2K(Ω − Ω 0 ) 2 , where Ω 0 ≡ −q 0 /2 is the rate of bundle twist preferred by chiral interactions. [24,26] Given this double twist texture, these elastic costs for intra-filament and inter-filament deformations may be derived for a given 2D cross-section. For example, since filament curvature grows linearly with distance from the bundle axis κ(x) = Ω 2 r, the bending energy (per unit volume) takes the simple form F bend /V = Bρ 0 Ω 4 r 2 /2, where r 2 is the moment of interia of the cross section. Calculation of the inter-filament elastic energy requires solution of the mechanical equilibrium equations, which are the Euler-Lagrange equations for u(x), ∂ i σ ij = 0 combined with a compatibility condition that relates inplane gradients of inter-filament stresses to in-plane gradients of filament tilt t ⊥ (in-plane components of t) where and plays a role equivalent to the Gaussian curvature in 2D elastic membranes, connecting gradients of out-of-plane tilt to gradients of in-plane stress. For the double-twist texture in eq. (6) K ef f = +3Ω 2 , where the small-tilt approximation assumes sufficiently large helical pitch such that (Ωr) 1, the frustration of inter-filament spacing is equivalent to crystalline order on a sphere of radius |Ω −1 |/ √ 3 [29]. In this model, consider two limiting classes of bundle cross-sectional shape: cylindrical and tape-like twisted bundles. Cylindrical bundles possess an isotropic (circular) cross-section, which would be optimal from the point of view of surface energy. For a cylindrical bundle of radius R the inter-filament elastic pressure has the form, σ ii = 3Y Ω 2 (R 2 − 2r 2 )/8, illustrating the tendency of compression to increase with distance from the bundle center, analogous to the geometric compression of a (flat) 2D membrane confined to the spherical surface. From this stress distribution and the second-moment of the cross section, r 2 = R 2 /2, we have In the second morphology, the helical tape, we assume the cross-section has a rectangular shape with one lateral dimension, the width w, much larger than its transverse thickness t. We approximate the stress distribution by taking the w t limit, where introducing Cartesian coordinates (x, y) along the respective width and thickness axes, stresses will be independent of x (away from distant ends at x = ±w/2). In this case, mechanical equilibrium and vanishing normal stress along the wide edges require σ yy σ xy 0. Combining compatibility with the condition of zero average stresses alongx direction, we find σ xx Y Ω 2 (t 2 − 12y 2 )/16. Notably, this solution shows that compression increases away from the line y = 0, leading to stress growth only along the thin dimension of the cross-section (as opposed to w). This illustrates the effect of introducing additional free surface to the frustrated assembly. Compatibility requires a nonzero stress drop across that bundle cross section, and the introduction of free surface allows that stress drop to take place along the narrow dimension of the domain through unconstrained displacement of the boundary, thereby lowering the elastic cost of frustration (at the expense of breaking cohesive bonds and raising the surface energy). Integrating the elastic stress through the section and using r 2 = (w 2 + t 2 )/12 we have have This shows the key result that elastic costs (per unit volume) of extrinsic (i.e. bending) geometry of anisotropic bundles are most sensitive to the wide dimension (growing as ∼ w 2 ) as compared to costs of intrinsic (i.e. metric) geometry which are controlled by the narrow dimension (growing as ∼ t 4 ).
To analyze the free energy density of competing morphologies it is useful to introduce two material dependent length scales: where the first length scale, the elasto-cohesive length, effectively characterizes the ratio of depth and elastic stiffness of inter filament bonds, and the second, called the bend penetration depth in the context of columnar order [43], characterizes the relative cost of intra-filament bending in comparison to inter-filament bond stretching. Further by introducing the dimensionless length scales, and energy scales, we arrive at dimensionless free energies for cylindrical and tape-like bundles where we neglect the τ 2 contribution to bending in the w t limit and we have dropped the constant contribution from the bulk chemical potential.

B. Cylindrical Bundles: Finite Diameter vs. Frustration Escape
First, we analyze the equilibrium conditions for bundles assumed to form cylindrical morphologies. Specifically, we analyze the role played by the strength of chiral inter-filament interactions, as parameterized by k, to control the ability to unwind the helical twist of the bundle (i.e. ω < 1) and allow the bundle to reach larger equilibrium sizes than possible for fixed twist, ω = 1. Here, we consider the canonical ensemble and assume conditions (e.g. saturated filament solutions) where all but a negligible fraction of filaments are condensed into bundles. In this case, assuming bundles adopt a mean radius ρ (with negligible size fluctuations), the preferred size and twist of bundles follows from minimization of f cyl with respect to ρ and ω yielding the equations of state, respectively and Solving the latter cubic equation for equilibrium twist we find, where Fig . 3A shows the radius dependence of equilibrium twist. Small bundles adopt preferred twist, ω cyl 1, while large bundles untwist as ω cyl ∼ (k/ρ 4 ) 1/3 as ρ → ∞ due to the size dependent elastic costs of twist frustration. The crossover between "twist-locked" and "mechanically-untwisted" bundles can be associated by the "half-twist size" ρ 1/2 defined by ω cyl (r 1/2 ) = 1/2 in eq. (17). For weak chiral interactions (k 1) ρ 1/2 ∼ k 1/2 indicating that bundles prefer to untwist upon reaching only a small radius, while for stiff chiral interactions (k 1) ρ 1/2 ∼ k 1/4 indicating that in the limit of clamped chirality (k → ∞) bundles retain preferred twist at all sizes.
Inserting the size-dependent twist, eq. (17), into eq. (15) yields the relation between surface energy and equilibrium bundle size ρ cyl . In the limit of s 1, bundles are driven by relatively weak cohesive interactions and retain preferred twist (ω cyl 1), reaching only a narrow size ρ cyl ∼ s 1/3 1 limited by intra-filament bending. The growth of bundles with increasing s is highly dependent on the stiffness of chiral interactions, k (see Fig. 3B). In the k → ∞ limit of "twist-locked" interactions, bundles crossover from s 1 to size limited by geometrically-induced inter-filament stress at large surface energy, ρ cyl ∼ s 1/5 1. For finite k, the reduction of equilibrium twist (shown in Fig. 3A) leads to larger bundle sizes for a given s, relative to the "twistlocked" limit, due to the untwisting of bundles with increasing size. This leads to a maximal s beyond which surface energy drives the bundle to completely untwist and reach infinite size. Thus, assembly can escape the costs of geometric frustration provided that the cohesive drive (as quantified by s) is sufficiently strong to untwist , plotting equilibrium twist of cylindrical bundles Ω (relative to value preferred by chirality, Ω0) as a function of outer bundle radius, for a range of reduced chiral stiffness, shows untwisting radius (defined by Ω = Ω0/2) increases with k. In (B), we show the selection of equilibrium cylinder size, R cyl , as a function of reduced surface energy , for a range of reduced chiral stiffness k. In (C), the outer filament twist angle tan θ = ΩR cyl of equilibrium cylinders, shows a generic non-monotonic dependence on surface energy. The inset in (C), shows the radius dependence of twisted cylinders near to the transition from finite-diameter bundles and bulk, untwisted assembly at s = s * for k = 10. In (B) and (C), dashed curves show metastable branches (for s * < s < smax) where finite cylinders are separated from (equilibrium) bulk assembly by a finite free energy barrier.
assembly to the parallel state. This follows from the fact that beyond r 1/2 bundle twist falls rapidly, and consequently the size-dependent mechanical costs fall (proportional to ω 4 cyl ) allowing bundles to grow to much larger sizes, peeling away from the k → ∞ behavior and ultimately exhibiting a singular divergence at a critical surface value s max . The k-dependence of this maximal s max can be estimated by using half-twist condition (ρ 1/2 and ω cyl = 1/2) in the equation of state, eq. (15), yielding s max ∼ k 3/2 for k 1 and s max ∼ k 5/4 for k 1. Thus, the stable range of self-limiting chiral filament bundle assembly vanishes as k → 0, and conversely, extends at all s for infinitely stiff chiral interactions.
Before considering the effect of "frustration escape" on the competition between cylindrical and tape-like bundles, we point out two features of the untwisting transition. Fig. 3C shows the non-monotonic evolution of the twist angle θ = arctan(ΩR) as a function of s for a range of chiral stiffness, k. For weakly cohesive assembly, θ generically increases with s due to the increasing radius and roughly constant pitch, while for sufficiently large s, bundles untwist substantially (say for s near the half-twist condition) leading to a maximum θ beyond which twist angle begins to decrease before a singular drop at s max . A second feature is illustrated in the inset of Fig. 3C, showing plots of the free energy vs. bundle radius near to s max , where f cyl is plotted relative to the bulk, untwisted (ω → 0, ρ → ∞) state. As surface energy approaches this value from below, the depth of the finite-ρ minimum increases until it reaches a critical value s * cyl < s max where finite cylinder assembly is in equilibrium with bulk, untwisted assembly. For s * cyl < s ≤ s max finite diameter bundles remain metastable in the regime of equilibrium bulk order (depicted as dashed curves in Fig. 3), while for s ≥ s max the barrier associated with untwisting to the parallel, unfrustrated state is removed. While bundle metastability is a generic feature for finite k assembly, the range of metastable s is relatively narrow, because s * cyl exhibits the same scaling dependence on stiffness of chiral interactions as s max .

C. Twisted Cylinders vs. Helical Tapes vs. Bulk, Untwisted Assembly
We now revisit a previous prediction for morphology transition between twisted cylindrical bundles and twisted tapes, which assumed fixed twist (i.e. fixed Ω) as a opposed to fixed chiral interactions (i.e variable Ω for fixed K and Ω 0 ). Previously, we showed that for fixed Ω, as Σ increases, the rapidly growing cost of interfilament frustration drives a morphology transition to highly anisotropic tapes. In this case, frustration by chirality is fixed, so that while tapes reach lateral dimensions (widths) that are larger than would be possible for cylindrical morphologies, they remain finite for all Σ due to the non-vanishing cost of finite twist. Hence, in the present model (fixed chiral interactions) chirality acts as a "soft frustration" mechanism, which can be overcome for sufficiently strong cohesive assembly drives by untwisting from the preferred chiral twist to parallel assembly.
Based on the continuum model, we consider the competition between three phases: cylindrical bundles (finite R and Ω), tape-like bundles (finite t, w and Ω), and untwisted bulk assembly (R → ∞ and Ω → 0). Along with the equations of state for cylindrical bundles, the free energy of tape-like bundles is determined via the solution for the equations of state of lateral bundle dimensions and twist (minimization of f tape with respect to τ , υ and ω), and The first two conditions are easily solved for a given twist υ = (24s/ω 4 ) 1/3 and τ = (80s/ω 4 ) 1/5 . When s ω 4 , due to the weaker growth of the bending energy (∼ υ 2 ) compared to the cost of inter-filament frustration (∼ τ 4 ), large cohesive energies cause tapes to grow increasingly more anisotropic, leading to w/t 1. Solving the additional condition for twist equilibrium, we find a qualitatively similar trend as described above for isotropic cylinders: tape dimensions increase with approximately constant ω tape 1 for small s, and at a critical size range equilibrium tapes rapidly unwind leading to a singular increase in tape size at critical s, the value of which increases with k. Indeed, because the thickness of tapes scales in the same way as cylinder radius for large s, the escape to bulk untwisted assembly (ω tape → 0) occurs at comparable, if somewhat larger, values of surface energy for tapes relative to cylindrical bundles.
The boundary between equilibrium cylinders and tapes is found by solving the condition f cyl = f tape for equilibrium structures, while the transition to bulk (untwisted) assembly is determined by comparing free energies of bundle morphologies to f bulk = 2k. For the k → ∞ limit of stiff chirality we find that reduced inter-filament frustration, at the expense of excess surface energy, in tapes leads to a transition from cylinders to tapes at s c/t (k → ∞) 3, 860, as previously reported in ref. [41]. For finite k, soft chiral interactions allow bundles, both tapes and cylinders, to untwist to bulk parallel assembly above a critical surface energy, as estimated by s max above. Whether this untwisting transition occurs from a cylinder phase or tape phase of bundles can crudely be determined by comparing the k dependence of s max to s c/t . This crude estimate suggest a critical range of k ≈ 100, above which stiff chiral interactions allow bundles to grow to sufficient size that inter-filament frustration drives a morphology transition to equilibrium tapes at large s. For "soft" chiral interactions (k < ∼ 100), frustration escape proceeds at sufficiently low s to prevent a broad range of stable equilibrium tapes before the bulk phase is reached at high s This crude picture is consistent with the 3-phase assembly phase diagram computed from exact solution of the equations of state above, as shown in Fig. 4. We note, however, that even in the regime of "soft" chiral interactions (k < ∼ 100) a narrow window of equilibrium persists between cylindrical bundles and bulk assembly due to the fact that bundle sizes grow dramatically with s near to the untwisting transition, generically accompanied by a sufficient rise in inter-filament frustration costs to drive a transition to tapes, albeit tapes which are also very close to the point of frustration escape via untwisting. A narrowing window of tapes persists down to a triple point at s tri = 6.65 and k tri = 2.4, below which very soft chiral interactions permit weak surface energy to drive finite cylinders directly to bulk, untwisted filament arrays.

III. RESHAPING THE BOUNDARY VS. REORGANIZING THE INTER-FILAMENT ORDER: DISCRETE-FILAMENT MODEL OF MORPHOLOGY
In this section, we investigate how frustrated chiral bundle assemblies select between two fundamentally distinct morphologies: defect-free bundles with anisotropic shapes (tapes) and bundles possessing topological defects in the interior, but with otherwise isotropic cross sections. Previous studies have investigated these distinct morphologies as independent structural responses to frustration in chiral bundles.
Disclinations are point-like locations in a 2D quasihexagonal packing where 6-fold inter-filament bond order breaks down. As shown schematically in Fig. 1E, 5-fold disclinations allow partial relaxation of the frustrationinduced compression along the azimuthal direction in twisted bundles [29,36], akin to how they relax hoop compression in spherically-curved crystals [46]. Defect configurations are characterized by their net topological charge Q, which is the sum of total deficit in coordination (e.g., assigning value +1/-1 to 5-fold/7-fold disclinations). Previous theory [36] and simulation studies [30,31] of twisted bundles with quasi-cylindrical boundary shapes show that the stability of excess disclinations in the ground state packing is controlled by the product ΩR ≡ tan θ, which corresponds to the area integral of Gaussian curvature of packings mapped onto the dual surface [29]. For bundles θ ≤ θ c = arctan( 2/9) 0.44, cylindrical bundles favor Q = 0 packings (largely defect free), while for θ > θ c ground states possess excess 5fold disclinations, Q > 0 with a number that increases from +1 to a maximum of +6 for highly twisted bundles (θ → π/2).
In a parallel study, Hall and coworkers [41] investigated the possibility for bundles to relax frustration without defects, through a reorganization of boundary shape. As described above, surface energy drives bundles to form larger diameters, which in turn increases the costs of frustration-induced elastic energy. Because the cost of metric frustration grows rapidly with the area of isotropic bundle cross sections A, growing as ∼ A 2 , relatively stiff inter-filament bonds ultimately lead to an elastic instability of boundary shape that drives cross sections to adopt anisotropic shapes (i.e. tape-like morphologies). This instability, which parallels the behavior of defect free crystal growth on spherical surfaces [38,40,47], lowers the elastic cost of inter-filament frustration by limiting the narrow dimension of the structure, at the expansion of increased surface energy and intra-filament bending cost. In this previous study, which considered the (fixed twist) canonical ensemble, we showed that this bundle to tape transition occurred at the critical value of surface energy, s c 3, 860 (simulation results showed a lower threshold, s c ≈ 20 − 60).
Here, we employ a discrete filament bundle simulation model to compare the relative stability of these competing morphological responses. For simplicity, we focus this study on the fixed twist and fixed filament number (microcanonical) ensembles. In contrast to previous studies, we compare numerically simulated bundle cross sections that alternately optimize boundary shape and interfilament order, and investigate the conditions that select between boundary shape vs. defect-mediated morphology response.

A. Discrete filament bundle model and simulation methods
We consider a discrete filament model of chiral filament bundles. Bundles are described by a fixed number of filaments N and fixed (post-equilibrium) value of twist Ω. The energy of bundle assembly is modeled by the energy of a cross-sectional slice of height ∆z, whose energy ∆E is described by The first term describes bending energy of filaments with κ i = Ω 2 ρ i /[1 + (Ωρ i ) 2 ] with ρ i the radial distance of the ith filament from the twist axis. The second term describes the cohesive interaction of filament i with filament j (by definition at an inner radial distance, ρ j ≤ ρ i ), which can be captured by the distance of closest approach ∆ ij between center lines of the filaments. Due to the 3D geometry of filament backbones, the closest contact between curves of i and j is generically out of the plane of the cross section, depending not only on distance between points in the cross section, but also on their distance from the central axis and the bundle twist Ω (see eqs. B8 and B9 in appendix B 1). Here, we choose the same class of cohesive potentials as in ref. [41] u(∆) = 6 5σ 11 where is the cohesive energy per unit length of filaments separated at an equilibrium distance d (e.g., their diameter). We choose a "5-11" potential since the form follows directly from the integration of a pair-wise Lennard-Jones potential between a point on i and a (locally straight) section of j, ∆ ij away, where the factor of (1 + κ j ∆ ij · N j ) −1/2 corrects for the curvature of j (with N j , its normal). Here, σ sets the ductility of pairwise bonds: σ d results in brittle inter-filament bonds and σ d results in ductile bonds. In eq. (21) ∆ i = ∆z 1 + (Ωρ i ) 2 is the length element of i in the cross section. From the energy per slice, we define the (intensive) energy per filament length, where ∆L = N i=1 ∆ i is the total filament length per cross-sectional slice. We focus simulations on the case of negligible bending energy content, and choose only a minimal non-zero value of B = 10 d 2 , so as to provide a slight inward force on filaments towards the bundle center.
To determine equilibrium morphologies for a given N , Ω, bending stiffness and inter-filament potential, we minimize E with respect to in-plane filament coordinates x i . As the energy landscape of this optimization problem is rough due to the underlying frustration, resulting equilibria are highly dependent of both initial configurations and sampling protocol. Rather than seek a rigorous ground state (say via advanced sampling), here we aim to identify candidate ground state structures which are assumed to closely parallel the true ground state in terms of energy and structure. As the structural competitors of ground states fall into at least two morphologically very distinct states (presumably separated by multiple large energy barriers) for each set of simulated bundle parameters we use a combination of the following methods.
Method A -This method considers filament positions which are initially defect-free and templated from an hexagonal lattice of spacing d, but fall within a crosssectional "stencil" of variable anisotropic structures. The stencil is described by a stadium curve with lateral inner dimensions, X 1 and X 2 , such that narrower ends of the stencil are capped by hemicircular arcs of radius |X 1 − X 2 |/2. The stencil dimensions and boundary filaments are adjusted such that the total number of enclosed filaments is N (see Appendix B 2), producing a distinct initial configuration that varies from nearly isotropic, cylinders (X 1 ≈ X 2 ) to highly anisotropic, tape-like bundles (X 1 X 2 ). From each of these distinct initial configurations, filament positions are first relaxed to avoid overlaps with ∆ ij < d (see Appendix B 2), and following this, relaxed according to conjugate gradient minimization of E, and the lowest energy (per length) is selected as the defect-free ground state candidate.
Method B -This method makes no assumption of initial triangular order and instead allows bundles to form with uncontrolled defect types, numbers and locations in the bundle cross sections. Bundle initial positions are generated by random placement of filament positions with a distance of 0.5d √ N from the bundle axis, followed by a post-initialization relaxation to remove overlaps (see Appendix B 2). From random initial configurations, filament positions are relaxed according to conjugate gradient minimization of E. For each set of bundle parameters, at least 3000 random initial configurations are generated, and the lowest final energy is selected as the ground state candidate.
Following application of both relaxation methods, ground state competitors are characterized according to lateral aspect ratio, α, defined by the dimensions of the rectangle with minimal area that encloses the cross section, including filament volumes. Additionally, the interfilament packing is triangulated, making use of the conformal projection from the dual curved surface to the 2D plane [31], bulk disclination defects are identified as positions of non 6-fold coordination not lying along the outer edge of the bundle, and the net topological charge Q is computed as the sum of total deficit in coordination (relative to 6-fold) from that set of defects.

B. Simulation Results and Discussion
Results of minimal energy morphologies for N = 181 filaments are shown in Fig. 5 for two different interaction ranges and a sequence of increasing twist, between Ωd = 0.015 − 0.134. For the relatively ductile potential σ = d, we find that bundle cross sections remain roughly circular, but adopt excess 5-fold disclinations above a critical twist (Q ≥ +1 for Ωd ≥ 0.075), and increasing topological charge with further twist, consistent with the results for the same potential in ref. [30]. In contrast, for the shorter-range, or relatively brittle, potential σ = 0.4d, at relatively low twist, before the onset of defects, optimal bundles become anisotropic. At Ωd = 0.03 the cross section becomes weakly anisotropic  N = 59 (B), where a threshold value of aspect ration α = 1.5 is used to distinguish cylinders from tapes and apparent aspect ratios of bars highlight optimal energy aspect ratios. Dashed line is drawn from predicted threshold for critical twist to stabilize 5-folds in cylindrical bundles and solid blue line is determined by continuum theory prediction for fixed N (and λB = 0) transition between defect-free cylinders and tapes (as discussed in text). α 1.49 (compared to the zero twist value of α = 1.04) and becomes increasing tapered as twist is raised to Ωd = 0.06 where α 2.94. For larger twist, optimal bundles return to relatively circular boundary shapes (defined according to α < 1.5), and punctuated with the same excess of 5-fold disclinations as the ductile potential states. Thus, for sufficiently brittle inter-filament potentials, anisotropic, defect-free bundles appear intermediate to defect-free cylinders at low-twist and defective cylinders at high twist.
In Fig. 6A, we present the morphology phase diagram of N = 181 bundles spanned by ductility, σ/d (ranging from 0.253 − 1.318) and reduced twist Ω √ A where A ≡ N ρ −1 0 with ρ −1 0 = √ 3d 2 /2 the area per filament in a uniform hexagonal packing. Consistent with the above, we observe that stable defect-free, anisotropic tapes are generally confined to a range of twist intermediate to defect-free and defective cylinders below a critical value of ductility, (σ/d) c 1.0, beyond which the boundary between cylinders and tapes shifts to lower twist with increasingly brittle interactions. In this range of potentials, the transition from defect-free bundles (either cylinders or tapes) to defective cylinders appears largely independent of ductility, occurring at a critical twist of roughly (Ω √ A) c ≈ 0.84. Assuming a circular symmetry, A = πR 2 , this critical twist is consistent with the previous prediction and simulation result for threshold twist at which the relaxation of geometric frustration per 5fold defect is sufficient to overcome the elastic self-energy for defect formation in isotropic bundles [36,37]. We note that transitions to higher topological charge packings (Q = +2, +3, +4) also appear to be relatively independent of ductility, occurring at roughly the same reduced twist for all values of σ/d.
In Fig. 6B, we show the morphology phase diagram for a smaller bundle, with N = 59, which shows the same generic trends as the larger bundle discussed above. While the threshold reduced twists for stabilizing defec- as shown in Fig. 7.
To understand the observed morphology transitions, in Fig. 8 we analyze the energetics of ground-state bundles for a sequence of increasing twist at fixed σ/d = 0.3 and N = 181 (i.e., sufficiently brittle to form tapes) and compare to continuum theory predictions at fixed N and B = λ B = 0 (see Appendix A). Specifically, we analyze the cohesive energy (relative to the Ω = 0 state), decomposed in the two pieces, E bulk and E surf , the energy of interior and surface filaments, respectively (see Appendix B 3). The energy of defect-free cylinders increases with small twist (both in the bulk and at boundary), in accordance with the elastic penalty due to frustration, ≈ Y A 2 Ω 4 , from eq. (8). For the fixed-N case, the transition to tapes can be estimated by the point at which the the inter-filament elastic energy of cylindrical bundles grows in excess of the surface energy cost, or Σ √ A > ∼ Y A 2 Ω 4 , indicating that it is advantageous to create additional surface in order to relax bulk frustration. A more careful analysis (see Appendix A) gives a reduced twist for the cylinder-tape transition ( These values of reduced twist are shown as blue curves in the phase diagrams in Fig. 6, and agree well with the defectfree cylinder-to-tape transition twist, and more significantly, its tendency to decrease with N and increase with σ/d, depending on the relative cost of inter-filament strain in comparison to bundle surface. As shown in Fig.  8, at an intermediate range of reduced twist (Ω √ A) cyl/tape < Ω √ A < (Ω √ A) c , anisotropic bundle shapes become favorable. From the continuum model (see Appendix A), the equilibrium aspect ratio (for α 1) is determined by a balance of excess surface energy, growing with narrow dimension Optimal aspect ratio, α, for simulation bundle ground states plotted versus a dimensionless ratio of predicted bulk-frustration to surface energy penalty for cylinders, A 5/2 Ω 4 /λS, for values of N = 59, 101 and 181. The dashed line shows large-α power law predicted from fixed-N continuum model, with morphology indicated as in Fig. 5 legend. Note that defective morphologies fall below α 1.5, the threshold we use to discriminate tapes from cylinders. w = √ Aα as ≈ Σ √ Aα, and inter-filament strain growing with thin dimension t = A/α as ≈ Y Ω 4 (A/α) 2 , yielding an equilibrium α * ≈ A(Ω 4 /λ S ) 2/5 that increases with filament number, twist and brittleness of inter-filament bonds. This power-law trend for tape aspect ratio is compared to simulated bundles in Fig. 9. Because equilibrium anisotropy increases with twist Ω, The surface energy E surf of tapes increases rapidly with twist, while E bulk remains relatively constant with Ω due to the mitigation of inter-filament strain by the narrow bundle thickness, as shown for stable tapes at intermediate-twist 0.46 < Ω √ A < 0.84 in Fig. 8. At high twist (above (Ω √ A) c 0.84), defective cylinders become stable, leading to a drop of surface energy due to the fewer number of boundary filaments in the circular cross section. For high twist, E bulk > E surf due to the straining of inter-filament bonds by the combined effect of frustration and topological defects. As reported previously in ref. [31] the steady increase of Q with Ω leads to a nominal plateau of energy density at high twist as interior defects effectively mediate the cost of increased frustration by twist. This generic flattening of the bulk frustration energy due to one or more excess 5-fold disclinations is responsible for the generic crossover from stable intermediate twist tapes, with energies that grow for modestly with twist (f tape ∼ Ω 4/5 as compared to f cyl ∼ Ω 4 for cylinders), to defective cylinders, which approach a relatively constant, if somewhat cuspy energy density due to the fixed circular boundary and neutralizing effect of 5-fold disclinations.
We use the predictions of the continuum model to present an overview of the fixed-N morphology selection in Fig. 10, here as a phase diagram spanned by the axes Ω 2 A and Ωλ S . The intersection of the cylindertape boundary Ω 2 A 13.2Ωλ S with the critical condi- tion for stabilizing 5-fold defects in cylinders Ω 2 A 0.71 yields the triple point where defect-free cylinder, defectfree tape and defective cylinder phases meet, giving also the prediction for the dependence of the critical ductility (σ/d) c on N , shown in Fig. 7. Here, we also show the effect of non-zero bending stiffness or Ωλ B = 0. The effect of bending stiffness in this morphology phase diagram is to penalize tape morphologies (which possess a larger second moment of cross-section area r 2 ) relative to cylinders, shifting the cylinder-tape boundary to larger values of Ω 2 A as shown in Fig. 10, particularly at low values of Ωλ S where brittle interactions otherwise favor very-wide tapes. This analysis shows that above a critical stiffness value, for (Ωλ B ) max ≥ 0.155, the window of intermediate twist tapes is shifted above the threshold for stable defects (which are insensitive to bending stiffness) and we therefore expect tape-like bundles to be preempted entirely by defective cylinders when filaments exceed this critical stiffness B max ≈ Ω −2 (d/σ) 2 .

IV. CONCLUSIONS
In this paper we investigate two models of morphology competition for self-assembled chiral filament bundles. The first considers the ability to "escape" frustration by untwisting chiral inter-filament skew and its effect on the otherwise self-limiting assembly of bundles. We find that finite stiffness of inter-filament chiral interactions, K, puts an upper limit on the (cohesive) surface energy Σ max ∼ K ν at which finite, frustration-limiting bundles form where ν = 3/2 for K K c and ν = 5/4 for K K c where K c = Ω 2 0 (ρ 0 B) 2 /Y . When chiral interactions are sufficiently stiff (K > ∼ 100K c ) equilibrium assembly includes windows of both cylindrical (small Σ) and tape-like bundles (intermediate Σ) as well as bulk, untwisted assembly (large Σ), while in the opposite limit of "soft" chiral interactions, the window of tapes narrows dramatically and ultimately vanishes, leading to direct transition from twisted cylinders to bulk assembly for even weakly cohesive bundles. The second model investigates the interplay between competing morphological responses to inter-filament frustration: bulk defects vs. boundary anisotropy. We found that when inter-filament bonds are sufficiently brittle, tape-like bundles are stable as an intermediate-twist phase between defect-free cylinder (low twist) and defective cylinder (high twist) morphologies, and that the selection of the energetically favorable morphology was controlled by the relative cost of straining inter-filament bonds in the bulk (as required to form large isotropic domains and topological defects) compared to the cohesive cost of creating excess free bundle surface, which has the critical effect of lowering frustration costs through low-energy boundary deformation.
The results demonstrate a surprising and generally overlooked properties of geometrically-frustrated assemblies [23]. Preferred chiral ordering makes uniform order, both intra-and inter-filament, impossible, leading to a competing set of heterogeneous, and morphological distinct, local free energy minimal. Thus, a proper understanding of even the equilibrium behavior of these systems is a far greater challenge than canonical, unfrustrated systems, requiring an analysis of all classes of metastable equilibrium, each of which is governed by a distinct set of thermodynamic or microscopic parameters describing assembly. In the context of the present work, an outstanding challenge still remains to compare and understand the competition of all four classes of morphological responses at once: self-limiting width, untwisting of chiral interactions, formation anisotropic boundary shape or topological defect in the interior.
Furthermore, a much broader challenge will remain to understand how the principles of morphology selection analyzed here also controls non-equilibrium assembly. For example, in [41] we analyzed the tape-like and cylindrical morphologies of twisted amyloid protein bundles in terms of the equilibrium (fixed twist) model of domain shape selection, and we find that tapes form at a degree of bundle twist, where based on the present study, we would expect equilibrium morphologies to relax frustration via one or more 5-fold interior defects while retaining a cylindrical shape. The observation of tapes in this regime suggests a kinetically limited assembly pathway, where the bundles are kinetically free to relax frustration upon reaching a critical size through anisotropic domain growth, but may be subject to very large barriers that inhibit the global reorganization of defect-free tapes into defective cylinders. A systematic study of the thermodynamically and kinetically accessible pathways between distinct morphologies of chiral filament bundles, or geometrically frustrated assemblies more broadly, remains a critical and unresolved aspect of the anomalous behavior of these systems.

Acknowledgments
The authors are grateful to acknowledge E. Matsumoto, G. Schröder-Turk and I. Bruss  Here, we describe the continuum model prediction for coexistence between cylindrical and tape morphologies, restricted to fixed twist rate Ω and fixed area A, equivalent to fixed filament number, N Aρ 0 . In this case, cylinder radius R and tape width w and thickness t satistfy A = πR 2 = wt. Defining dimensionless areaĀ ≡ Ω 2 A, elasto-cohesive lengthλ S ≡ Ωλ S , and bend penetration depthλ B ≡ Ωλ B , we rewrite eqs. 8 and 9 (neglecting the contribution from twist elasticity for fixed-twist ensemble and re-casting energy density f ≡ F/Y V ), The equilibrium value of aspect ratio α * is determined by the condition where we consider only the stable solutions where aspect ratio α * > 1, since the inter-filament elastic energy fails to hold for α < ∼ 1. For example, if we consider the limit of α 1 and λ B = 0 (as discussed in the text), we may neglect the bending terms and the surface contribution from the narrow ends of the tape section (∝ α −1/2 1), relative to the large contribution from the wide dimension (∝ α 1/2 1). From this limit of A3 we find equilibrium tape dimensions and energy, α * =Ā (80λ S ) 2/5 , f tape (α * ) = 5 4/5 2 9/5λ Coexistence is determined by the solution to f cyl (α * = f tape . Exact numerical solution of equilibrium conditions for λ B = 0 gives coexistence between cylinders and tapes along the line . With increasing bending cost, tapes become less stable due to the distribution of filaments to larger radial distances from the twist axis, and therefore, higher curvature. Notably, for Ωλ B > 0.155, tapes are unstable for all Ω 2 A below the threshold 2π/9 which defines coexistence of the defect-free cylindrical morphology with defective cylindrical morphology with (+1) topological charge (equivalent to (ΩR) c = 2/9 derived in [36]). The triple point where defective structures coexist with both defect-free tapes and cylinders, is found according to the above conditions for coexistence, subject to the constraint (Ω 2 A) c = 2π/9 0.7. Defining and for outer filament length i : Thus, the total energy depends on the length of the outer filament and the second-order deflection of j near to close contact with i, accounting for the factor of (1 + κ j ∆ ij · N j ) −1/2 . Taking ∆ i as the discrete approximation for ds i we arrive at the formula for ∆E ij used in eq. (21). When the interaction has the form of a Lennard-Jones potential, V (r) = V 0 [(σ 0 /r) 12 − 2(σ 0 /r) 6 ], then u(∆) = 6 [5(σ/∆] 11 − 11(σ/∆) 5 ] with preferred filament spacing d = σ = 0.947σ 0 and minimum energy = 1.686V 0 σ 0 [31]. Hence the range (and curvature) of the cohesive potential is locked to the equilibrium spacing for this particular model of cohesive "threads". Here, we use a broader class of inter-filament potentials, u(∆), that allow for independent adjustment of equilibrium spacing, d, and range, σ of attractive well: At the minimum ∆ = d, the curvature, or stiffness, of the potential is a function of σ only, The above form of this potential at the minumum controls the coarse-grained modulus of the inter-filament array (see [31,41]), Y 110 /( √ 3σ 2 ), while the filament diameter controls the surface energy Σ /d, allowing for a broadly tunable range of elasto-cohesive length λ S = Σ/Y The inter-filament spacing ∆ is computed according to an approximation previously established [30]. In short, for filaments with centerlines intersecting the cross section at positions x i and x j , the exact expression for interfilament distance comes from minimization with respect to the vertical component of separation z of the following: To approximate the contact separation z * , we use the following interpolation, This form interpolates between the exact solutions to d dz |∆ ij (z)| 2 = 0, in the limits of small and large Ω 2 ρ i ρ j and has been shown to be sufficiently accurate for potentials whose range is smaller that their thickness, a limit well satisfied by simulations here. 3 tapes are favored, while for the relative ductile case, σ/d = 1.0, roughly cylindrical bundles (α 1) are optimal.

Simulation Details
Method A, Lattice Sampling -Initial configurations were generated with filament positions on a triangular lattice with spacing d, at points within a stadium perimeter (connecting two half-circles with straight lines as in Figure 1A). The stadium aspect ratios were taken ranging from 1.0 to 7.5 in increments of 0.2, and taken with center on a lattice point, lattice edge and the center of of a triangular plaquette. For aspect ratio larger than 1, the longer dimension (tape width) was oriented along the close-packed lattice direction. Stadium dimensions were chosen to be the smallest so that each choice of as-pect ratio and stadium center enclosed the desired number of filaments. When this produced a structure with excess filaments (i.e. close to, but more than predetermined N ), additional points were omitted from the last close-packed row, starting from the end. This filament "pruning" is illustrated in Fig. 11A, where two outer points are omitted to produce the structure with 181 filaments. Fig. 11B shows the dimensions, and two limiting example structures, of lattice initialized configurations sampled for N = 181 filaments. Corresponding relaxed energies for each lattice structure sampled are shown in Fig. 11c for a twisted bundles at two sample values of ductility σ/d (brittle/ductile).
Method B, Random Sampling -Random initial configurations were generated by uniformly distributing filament positions within a circular boundary of radius 0.5 √ N d, such that the starting configuration always had slight compression. In Fig. 12, we show how these randomly seeded configurations lead to competing energy minima with a range of Q for a given bundle parameter values. For the purposes of our present study, we chose to exhaustively sample the defective configurations for N = 59 and N = 181, and sample defect configurations in N = 101 and N = 440 bundles more sparsely, only to confirm the transition from defect-free cylinders/tapes to defective packings at the critical twist. Accordingly, the number of randomly-initial configurations, n init per set of bundle parameters (Ω and σ/d) was varied with N : n init = 6000 for N = 59; n init = 3000 for N = 101; n init = 60000 for N = 181; and n init = 3000 for N = 440. For N = 1540, no attempt was made to sample random initial structures, due to the prohibitive computational cost for bundles at this size.
For each value of filament number N , simulations were run with twist Ω √ N ranging from 0 to 2.0 with 30 equally incremented values. The range of ductility σ/d varied with filament number N , and values of σ/d were spaced logarithmically with successive ratio 1.096: Structure Relaxation -From each set of initial filament positions, a local equilibrium was found by steepest descent according to the defined energy density, stopping when the largest component of the gradient (parameterized by filament coordinates) had magnitude less than 10 −4 /d. In areas of low twist, the best structures were generated by slight relaxation of configurations generated from a triangular lattice. Larger rearrangements were required in regions of high twist, and randomly generated starting configurations consistently produced more favorable structures. In order to allow larger rearrange- ments even with brittle systems and to remove overlaps with filament "cores" (∆ < d − σ), initial structures were first (before steepest descent) preconditioned by relaxation via steepest descent according to a soft repulsive potential,

Simulation analysis
Morphology Characterization -The lowest-energy structure at each point on the phase diagram was analyzed, taking both the aspect ratio and net topological charge of the interior. The aspect ratio was taken from the lateral dimension of the minimum-area rectangle that encloses the cross section. The minimal-area rectangle is found from considering all rectangles enclosing the structure's convex hull where one edge of the rectangle must coincide with an edge of the convex hull.
The net topological charge was computed by Delaunay triangulation of the transformed cross-sectional positions. As the true distances between filaments are properly encoded in the dual metric surface, we triangulate based on a conformal map of filament positions from the dual surface to the 2D plane . This map is performed by radial dilation of filament position, ρ = ρ 1+ √ 1+(Ωρ) 2 e √ 1+(Ωρ) 2 for a filament with crosssectional radial position ρ. This mapping results in the angles between filament separations reflecting the true inter-filament bond angles when measuring filament distances according to the geodesic distance on the dual surface [30]. The conformal map allows regular Delaunay triangulation, subdividing the convex hull of the packing into triangles connecting filament positions, so that the minimum inter bond-angle of the triangulation is maximized. The perimeter of the bundle was constructed by starting with the convex hull, then successively removing exterior bonds with large distance of closest approach, |(∆ ij /d)| > 1.4 and adding the interior point that shared a triangular plaquette according to the Delaunay triangulation. The net topological charge of the interior was then taken by totaling the excess coordination for each filament not belonging to the set of boundary filaments.
Surface and Bulk Cohesive Energies -Here we decompose the cohesive contribution from pair-wise interactions, into energies associated with surface formation and bulk inter-filament deformation. The surface energy is taken to be the deficit of exterior filament contact from a reference state with total cohesive energy 3 . The surface energy is then The bulk energy is taken as the difference, ∆E bulk = ∆E coh − ∆E surf .