Defect-driven shape instabilities of bundles

Topological defects are crucial to the thermodynamics and structure of condensed matter systems. For instance, when incorporated into crystalline membranes like graphene, disclinations with positive and negative topological charge elastically buckle the material into conical and saddle-like shapes respectively. A recently uncovered mapping between the inter-element spacing in 2D columnar structures and the metric properties of curved surfaces motivates basic questions about the interplay between defects in the cross section of a columnar bundle and its 3D shape. Such questions are critical to the structure of a broad class of filamentous materials, from biological assemblies like protein fibers to nano- or micro-structured synthetic materials like carbon nanotube bundles. Here, we explore the buckling behavior for elementary disclinations in hexagonal bundles using a combination of continuum elasticity theory and numerical simulations of discrete-filaments. We show that shape instabilities are controlled by a single material-dependent parameter that characterizes the ratio of inter-filament to intra-filament elastic energies. Along with a host of previously unknown shape equilibria---the filamentous analogs to the conical and saddle-like shapes of defective membranes---we find a profoundly asymmetric response to positive and negative topologically charged defects in the infinite length limit that is without parallel to the membrane analog. The highly non-linear dependence on the sign of the disclination charge is shown to have a purely geometric origin, stemming the from the distinct compatibility (or incompatibility) of effectively positive- (or negative-) curvature geometries with lengthwise-constant filament spacing.


I. INTRODUCTION
Topological defects are fundamental to the properties of ordered materials, from their structure and thermodynamics to their dynamics and mechanical response. There is long history, dating back to some of the earliest mathematical models of defects [1][2][3], of understanding the non-linear influence of topological defects on material structure through the lens of differential geometry. In such a description, topological defects are understood as sources for metric deformation in solid media [4], specifically curvature multipoles, leading to intrinsic stresses that reshape the material and its stress response. Far more than descriptive, the relationship between intrinsic (Gaussian) curvature and topological defects, becomes even more profound for ordered systems that are free to "reshape their metrics", such as 2D ordered membranes of both crystalline and liquid crystalline (e.g. 2D in-plane polar, nematic or smectic order) varieties [5][6][7][8][9]. The flexibility of out-of-plane deformations in such systems gives rise to an instability, in which the in-plane metric of a thin membrane may adapt to its 3D environment to accommodate the non-Euclidean geometry favored by disclinations (and multi-pole combinations thereof, such as dislocations), leading to a spontaneous buckling of sufficiently thin membranes [6]. Beyond the relevance to self-organized matter, the mechanical and geometric principles of defects in elastic sheets is now a cornerstone of the current approaches for designing and engineering 2D origami and kirigami materials [10][11][12][13].
In this paper, we analyze the defect-induced geometric instability of a parallel class of 2D ordered matter: columnar or filamentous bundles. This is done using a structural model that describes cohesive assemblies of many quasi-1D elements (e.g. filaments or columns) possessing 2D order transverse to their backbone. This model applies to a broad class of materials, from assemblies of flexible filaments (e.g. protein filaments [14][15][16][17][18] or synthetic nanotubes [19][20][21][22]) that self-assembled into cable-like structures via attractive interactions between fibers of self-stacking molecules [23,24], to finite domains of columnar liquid crystals [25][26][27], or even condensed phases of vortices on multicomponent superconductors [28]. This class of material is two-dimensionally ordered in the sense that it retains translational symmetry along the filaments or columns [29]. Recent works [30][31][32] demonstrate that columnar systems, like membranes, are also capable of altering their geometry by modifying their 3D embedding. But unlike membranes which deflect outof-plane, this deformation occurs through the geometrically non-linear coupling between the columns' orientations and transverse spacings (i.e. their metric). A well known example of such a coupling in bulk columnar media is the Helfrich-Hurault instability [33], in which a sufficiently large transverse tension drives a non-uniform tilt of the columns, thus maintaining a more uniform intercolumn spacing at the expense of bending [34]. Driven by the mechanics of metric deformation, this instability is the columnar analogy to the Euler buckling of a 2D elastic sheet under compression, which underlies the 3D buckling behavior of defective crystalline membranes.
In contrast, for columnar media with defects, the consequences of the geometric instability to tension are not known.
We exploit this analogy between 2D crystals and columnar materials to characterize the structural instabilities triggered by the stress from topological defects in the cross-sectional order of bundles. Specifically, we consider the instabilities driven by elementary 5-fold (positive) and 7-fold (negative) disclinations in otherwise hexagonally ordered bundles, characterized respectively, by the removal or insertion of a 60 • wedge of crystalline material [29]. Just as in 2D crystals, disclinations and dislocations, can be characterized via a Volterra construction corresponding to the mismatch of lattice rotation and displacement around a closed loop encircling the defect. Disclinations are quantified by a topological charge, s, measuring the angular turning of lattice directions around the defect, which must be integer multiples of 2π/6 so that the array remains 6-fold at all points except defect cores. The analogous problem, the shape of hexagonally ordered membranes with 5-and 7fold disclinations, was studied by Seung and Nelson, in the context of the Föppl-von Kármán (FvK) theory of crystalline sheets [6]. Creating a 5-fold (7-fold) disclination requires removing (adding) a wedge of crystalline material from the sheet, which stretches (compresses) distances azimuthally around the defect, yielding tensile (compressive) stresses along the hoop direction. These stresses are mechanically balanced by compressive (tensile) stresses along the radial lines extending from the disclination. The fact that thin elastic sheets are unstable to compression now justifies the shapes favored by disclinations: Five-fold defects favor conical shapeswith positive Gaussian curvature-buckled along the radial lines extending from the defect; while 7-fold defects favor saddle shapes-with negative Gaussian curvaturebuckled along the azimuthal hoops surrounding the defect (shown visually in Fig. 1(a)). The buckling behavior of a crystalline sheet of radius R is governed by a single dimensionless number, the Föppl-von Kármán (FvK) number, γ s ≡ Y R 2 /B, which characterizes the crystal's relative resistance to in-plane stretching vs. out-of-plane bending, described by the elastic moduli Y and B, respectively. Disclinated crystals remain flat for small γ s , but become unstable above a threshold value (which is somewhat higher for 7-fold than 5-fold), at which point in-plane stretching exceed the cost of bending to yield a buckled 3D shape.
Complementary to this mechanical perspective is the geometrical one, in which topological defects redefine the metric of the 2D surface in which the crystal is embedded [4]. Consistent with the Gauss-Bonnet theorem [35], the deficit or excess angle associated with disclinations can be accommodated without far-field strain, provided that it is balanced by the integrated Gaussian curvature of the sheet. This is illustrated in Fig. 1(b), where the relative lengths of the radial (dashed black) and circumferential (solid purple) paths along the surface depends The differing geometries can be characterized by comparing the radial distance between two disks R (dashed black line), and the circumferential distance from an outer disk to itself C (solid purple line). For a flat surface C = 2πR (center), but for a positively curved surface C < 2πR (left), while a for a negatively curved surface C > 2πR (right). (c) Equivalent geometries exist for flexible columnar materials when we consider the distance of separation perpendicular to the filaments. A pattern of twist (left) reduces the circumferential distance between filaments (solid purple), similar to a dome. Alternatively, a pattern of splay (right) reduces the radial (dashed black) distance between filaments, similar to a saddle. (d) Based on these geometries, conjectured structures are shown for a bundle with a 5-fold disclination (left), and 7-fold disclination (right). on the curvature. Similarly, in a columnar material, variations in filament orientation can be linked to geometrical constraints on their spacing [30,32], which in turn has a precise connection to the metric geometry. This can be seen in the simplified depictions of filament bundles in Fig. 1(c), which also shows the radial (dashed black) and circumferential (solid purple) distances between filaments. The result is a unique equivalence between a pattern of filament orientation and a corresponding surface with a Gaussian curvature of where t x and t y are respectively the x and y components of filament tilt in the plane normal to their mean orientation,ẑ. Given the generic instability of columnar structures to the internal tensile stresses generated by disclinations, it is reasonable to expect that sufficiently flexible bundles will buckle into non-parallel 3D shapes. However, it is a priori unclear exactly what 3D patterns of orientation will be trigger by such defect-generated stresses, nor what parameters control their relaxed shapes (as γ s does for crystalline sheets). Beyond their mechanical and geometrical correspondence as 2D ordered materials, crystalline membranes and columnar bundles have obvious and profound differences. Specifically, bundles are fully three-dimensional structures, i.e. their full degrees of freedom are not reducible to a 2D manifold. And it remains to be understood how their buckling behavior relates to their well-studied 2D membrane counterparts. For example, assuming the pattern of tilt to be axisymmetric in response to a centered disclinations, one might anticipate that 5-and 7-fold defects generate the respective 3D double-twist (left) and undulating splay (right) tilt patterns shown in Fig. 1(d) [32]. As we describe below, the distinction between such shape modes gives rise to vastly different elastic energies. This reveals a fundamental asymmetry between the ability of bundles to realize the analogs of positive-or negative-curvature metric geometries; and as a consequence, we observe a profoundly asymmetric response to these two elementary defect types.
In this paper, we employ a combination of continuum elastic theory and discrete-filament simulations of a minimal model of cohesive bundles to study the shape transitions driven by single elementary 5-or 7-fold disclinations. Based on an axisymmetric model of infinite-length bundles, we show that the buckling behavior of bundles is governed by a parameter that we a call the filamentaryvon Karman (fil-vK) number. Analogous to the FvK number for membranes, the fil-vK measures the dimensionless ratio of inter-columnar distortions (imposed by defects) to the cost of the lengthwise bending of filaments. However, unlike the case of crystalline sheets, we show that 5-fold defects lead to shape buckling of all bundles, without a threshold for size or filament flexibility. In contrast, we show that 7-fold defects are characterized by a finite fil-vK number instability threshold for axisymmetric splay undulations (i.e. finite bundle diameter or filament stiffness). We show that this dramatic asymmetry derives from the existence of a uniquely soft torsional mode available to 5-fold defects that generates an equivalent positive curvature without lengthwise variation in strains. Whereas for 7-fold defects, no such shape mode exists that provides "negative-curvature" without breaking longitudinal (i.e. lengthwise) symmetry, and thus amplifying inter-column strains. Indeed, we find that this additional frustration of negative-curvature underlies a much more profound symmetry-breaking for 7-fold bundles that can be described with the simplistic assumption of axisymmetry. Comparing our continuum analysis to discrete-filament simulations we find that for 5-fold defective bundles composed of very flexible filaments, there is a spectrum of metastable torsionally-wrinkled modes. While for 7-fold defective bundles, we find significantly more complex non-axisymmetric modes that compromise the drive for "negative-curvature" with the preference for uniform inter-filament strains along the length. We argue that these low-symmetry "counter-twisted" tilt patterns allow for a significant reduction of the threshold diameter for buckling 7-fold defective bundles, yet a finite threshold must remain in the infinite length limit due to the breaking of lengthwise symmetry. Nevertheless, these results are consistent with the qualitative distinctions captured by the asymmetric shape modes. Finally, we conclude with a discussion of fil-vK number values and their ramifications for various experimental systems of cohesive filament bundles.

II. AXISYMMETRIC SHAPE INSTABILITIES
In this section we explore a continuum elastic description of a columnar bundle possessing a 5-or 7-fold disclination in its cross section. We then analyze how these defects trigger axisymmetric shape instabilities. While we will show in Sec. III, that the assumption of axisymmetry ultimately fails for 7-fold defective bundles due to frustration between negative-curvature geometry and longitudinal symmetry, the analysis of the axisymmetric shape modes most clearly illustrates the mechanical and geometric principles that underly a profoundly asymmetric response to positive vs. negative disclinations. Furthermore, it highlights the critical combination of the material parameters that govern the shape response of bundles to defects. Our analysis examines the case of a hexagonally ordered columnar material, which can be considered to be a generic example of a filament bundle.

A. Continuum elasticity of defective bundles
Consider an initially cylindrical bundles with a radius R and length L → ∞ along theẑ axis. Here, as in refs. [30,36], the stress-free reference state is described as a 2D hexagonal array of parallel filaments. Deformations are measured relative to the reference configuration by the displacements u ⊥ (x) of local filaments at a point x in the bundle. As there is no cost for lengthwise displacements in such a material, u ⊥ (x) is 2D and perpendicular toẑ. The elastic energy for deformation is where λ and µ are the Lamé coefficients deriving from inter-filament cohesive forces [33,37], and related to the 2D Young's modulus, Y = 4µ(λ + µ)/(λ + 2µ), and Poisson ratio, ν = λ/(λ + 2µ) of the filament array. The strain tensor has components in the xy plane [38] where t i is the in-plane component of the filament tangent unit vector. In the limit of small tilt (|∂ z u ⊥ | 1), the tangent is Additionally, we consider the elastic cost of lengthwise gradients of t(x) associated with filament curvature κ = |(t · ∇)t| [39] where the value of the Frank constant, K, is proportional to the intrinsic bending modulus of filaments, B. The ratio of intra-to inter-filament elastic moduli defines a length scale λ 2 b ≡ K/Y , typically associated with the penetration depth of bending deformation in a columnar material [34]. Comparing this size scale to the lateral size of the bundle defines the dimensionless fil-vK number, Analogous to the FvK number for thin membranes γ s , γ assesses the relative costs of inter-filament vs. intrafilament deformations in the bundle, and as we show below, is critical for regulating the buckling of unstable bundles.
To explore the connection between the defect-induced instabilities of crystalline membranes and columnar bundles, we first consider the Euler-Lagrange equations of E = E elas + E bend , for the case of a tilt pattern t(x) that is fixed along the length (δt(x) = 0). These are simply conditions of in-plane force balance, with a stress of σ ij = δ ij λu ii + 2µu ij . Like the case of a membrane with a fixed topography, the relaxation of the displacement u ⊥ for bundles can be derived from the conditions of in-plane force balance augmented with a compatibility equation. This equation enforces the stress contributions from both the in-plane components of t(x) as well as singularities in the displacement fields associated with topological defects in the 2D crystalline order, where s(x) = α s α δ (2) (x − x α ) is the areal density of a topological disclination charge (s α is the charge and x α is the position of the disclination α), and K eff is the Gaussian curvature of a surface that approximates the inter-columnar metric of the 2D bundle cross section [32]. Illustrated visually, Fig. 1(b) shows a 2D surface that approximates the inter-filament distances found in Fig. 1(c). Hence, both topological defects (s(x)) and filament tilt patterns for which K eff = 0, act as sources for far-field inter-filament stresses. This relation has been previously used to show, for example, that positively charged (5fold) disclinations become stable for bundles with fixed and sufficiently large double-twist (i.e. the tilt pattern shown on the left in Fig. 1(d)). However, there is a complication in determining the buckling modes for negatively charged (7-fold) disclinations: this unusual case of double-twist is the only tilt pattern that yields a constant strain along the bundle's length; but this pattern alone creates the wrong effective curvature (positive rather than negative). Alternatively, while a locally-splayed geometry where K eff < 0, may partially neutralized a negatively charged (7-fold) disclinations-the right image in Fig. 1(d)-this tilt pattern unfavorably breaks lengthwise symmetry (i.e. because ∂ z u ⊥ = t ⊥ ). In Sec. II B, this axisymmetric splay pattern is used to determine the generic γ-dependence of buckling 7-fold defective bundles. Though, in Sec. III, our discrete filament simulations reveal that frustration between negative curvature and longitudinal symmetry leads to a far more complex and lower energy tilt pattern that preempts the axisymmetric instability for 7-fold, but not 5-fold, defective bundles. We will show that the breaking of lengthwise symmetry cannot be avoided for 7-fold defective bundle, and thus the qualitative and profound distinctions between 5-and 7-fold defective bundles predicted by the continuum analysis, are still borne out by simulations with unconstrained shapes.
Our purpose is to understand the equilibrium patterns of displacement (and thus orientation) that result from fixed topological defect structure. Therefore, to accurately determine mechanical equilibrium through the Euler-Lagrange conditions, we must consider lengthwise variation of u ⊥ and the associated variation of t ⊥ ∂ z u ⊥ . These more general equilibrium conditions used in the instability analysis below take the form, Relative to the case of a fixed tilt pattern in eqn (7), this force balance introduces two additional terms. The first term, ∂ z t i σ ij , couples stresses in consecutive "layers" of the bundle, and it is the analog of the "Young-Laplace" contribution to the normal force (i.e. in the first FvK equation) from in-plane stresses in curved membranes. The last term, proportional to K, derives from torques generated by bending of the filaments, which are expressed here as in-plane forces.
B. Linear stability of parallel, defective bundles Given this foundation of a continuum elastic model of columnar materials, we now employ linear stability to determine the buckling patterns caused by centered disclinations. The results will reveal a fundamental difference in the charge (7-vs 5-fold) dependence of deformation, owing to symmetry breaking in the lengthwise direction, an aspect unique to columnar materials.
To begin, we consider the stability of an initially parallel bundle (t 0 =ẑ) possessing a centered disclination, whose equilibrium stress σ 0 ij satisfy eqns (8) and (9), and σ rφ = 0, where r and φ are respectively the crosssectional in-plane radial and polar angle coordinates, and s = ±π/3 is the topological disclination charge (where ± refers respectively to 5-and 7-fold defects). From this base displacement field, u 0 , generated by the defect (associated with σ 0 ij ), we apply non-parallel displacements of δu(x), such that u(x) = u 0 (x) + δu(x). In particular, we consider deformations that are periodic along z and axisymmetric in the plane, These two periodic shape modes we refer to as, splayand torsional wrinkles respectively, for k = 0. We can now analyze the instability of a parallel bundle to splay (δu r = 0) or torsional (δu φ = 0) shape modes derived from the existence of solutions to the force balance equations. Naturally, we consider the limit of small deflections from the initial parallel state; or in other words, the solution to eqn (9), to linear order in δu(x). The (linearized) force balance equations can be recast in a simple and surprisingly familiar form (see Appendix B), where is the radial part of the 2D Laplacian, and the exact forms of V α (r) and α are given in eqs. (B13) and (B15). The "potential", V α (r) ∝ −s(kR) 2 (ln r + C 0 ), derives from the coupling of tilt to the defect-induced stress, while the eigenvalue, α ∝ (kR) 4 /γ, derives from forces induced by bending. Thus, the stability of defective bundles to axisymmetric splay-or torsionally-wrinkled shapes are formally equivalent to finding (zero angular momentum) bound states of a 2D hydrogen atom "energy" − α , whose central "charge" is sk 2 . The boundary conditions are determined by the condition of finite stress at the center of the bundle, or δu(0) = 0, and vanishing stresses at the outer bundle surface, which take the form for radial stresses and azimuthal stresses While superficially similar in form, the distinct boundary conditions underly profound difference between splay and torsional deformations of 2D columnar materials. Torsional modes with "zero kinetic energy", that is δu φ ∝ r, generate no shear stress at the bundle surface, while the same is not true corresponding to radial modes, which indicate that splay ground states acquire "kinetic energy". As a result, for splay modes the existence of "bound states", where α ≥ 0, occurs only for finite k where V r (r) is sufficiently strong, while for torsional modes bound states exist for all wave vectors down to k → 0.
The results of this ground state analysis is shown in Fig. 2 (for ν = 1/3, chosen to make comparisons to our discrete model later on). Here γ c (k) shows the critical fil-vK number, above which 5-and 7-fold defective bundles are unstable to torsional and splay instabilities at a wave vector k. The distinct wave vector dependence of these instabilities follows from simple energetic arguments: Consider first, a torsional mode, δu φ ≈ τ 0 r cos(kz), where τ 0 is a constant. Because this mode is a purely rigid rotation aroundẑ, to linear order elastic strains vanish (i.e. ∂ i u j + ∂ j u i = 0), and the only elastic strains are generated by tilt δt ≈ τ 0 krφ, leading to a mean strain δu φφ ≈ −k 2 τ 2 0 r 2 that is compressive in the hoop direction. This leads to a relaxation of defect induced stress of dA σ 0 φφ δu φφ > Y sk 2 R 4 τ 2 0 , which is dominated by large r where σ 0 φφ (r)/s > 0. Combining this with the bending cost gives an energy density ε tor for torsional-wrinkling (relative to the parallel state) This shows that the relaxation of tensile strain generated by 5-fold defects (s = +π/3) exceeds the bending cost for modes, k < k c ≈ λ −1 b , or γ c ∼ (kR) 2 , with the scaling shown in Fig. 2(a). Hence, bundles of any size or stiffness are unstable to long-wavelength (k → 0) torsionalwrinkles. In infinite bundles, such modes correspond to uniform helical twist δt Ωrφ, studied previously as an ansatz for elastic energy ground states in the presence of 5-fold disclinations [30].
As L → ∞, this lack of a threshold for the shape instability driven by 5-fold defects is unlike the analogous problem of conical buckling of crystalline sheets. This difference derives from the fact that in the latter case, the elastic energy released by conical buckling is proportional to the square of the sheet curvature (i.e. the Gaussian curvature), just as the (positive) energy cost of bending. Hence, membrane stiffness must fall below a critical value for shape buckling. However, for cohesive bundles, the specific shape mode driven by positive disclinations is a soft mode, generating elastic elastic costs only at O(Ω 4 ), which for small twists, are always overwhelmed by the elastic energy released by twist, proportional to sΩ 2 . (a) A 5-fold disclination yields a stability range that decreases to zero in the long-wavelength limit of k → 0, i.e. homogeneous pitch. (b) Alternatively, a 7-fold disclination is unstable to splay undulations. The minimum γ unstable mode (black dot) occurs for kR = 21.8, and γc = 13, 685. The inset shows a schematic comparison between the axisymmetric instability and the γc corresponding to the non-axisymmetric shape modes explored by the discrete filament simulations (green dot), analyzed in Sec. III B and extrapolated to the L → ∞ limit.
Turning to the case of 7-fold defects (s = −π/3), the energetic analysis proceeds along similar lines for a shape mode ansatz of radial splay, which we approximate by the linear profile δu r ≈ ρ 0 r cos(kz), where ρ 0 is a constant. However, unlike torsional wrinkles, splay leads not only to tilt-induced radial strains (∝ (kρ 0 r) 2 ), but also to linear in-plane strains, due to area dilation ∇ ⊥ · δu ∼ ρ 0 cos(kz), generating an elastic cost (per unit volume) of ≈ Y ρ 2 0 . Thus, the cost of splay modes has the form, (16) where the cost at k → 0 derives from the well-known coupling between splay and lengthwise density variations in columnar systems. The radial splay modes that relax elastic energy (due to the collapse of the radial tension generated by 7-fold defects) are unlike twist in that they are not soft modes. The balance between relaxing bending and elastic energy selects an optimal wrinkling wavelength k c ≈ |s| 1/2 λ −1 b , with a net relaxation proportional to −Y s 2 γ. Thus, only when the fil-vK number exceeds a threshold value, will the 7-fold defect drive (finite k) splay-wrinkling of the bundle. This result for 7-fold defects is shown in Fig. 2(b). Unlike the torsional wrinkling in the presence of a 5-fold defect, which becomes unstable at long-wavelengths, now there is a range of long-wavelengths (small k), for which no unstable solution exists at any γ. The minimum unstable value of γ c , occurs at a mode kR = 21.8, for which γ c = 13, 685, setting an upper limit threshold fil-vK number for the splay instability driven by a 7-fold disclination. For large k, the stability line again follows the scaling of γ c ∼ (kR) 2 .
The stability analysis to axisymmetric shape modes illustrates a profound asymmetry between the response to 5-vs. 7-fold defects in the infinite length limit. Bundles with centered 5-fold defects are always unstable; while for 7-fold defects, parallel bundles are stable up to a finite γ (proportional to their lateral area), beyond which they become unstable to lengthwise shape modes at a finite wavelength. This is in sharp contrast to crystalline membranes, where the modest asymmetry in the bending cost of respectively conical vs. saddle-like shapes driven by 5-and 7-fold defects lead to only slight difference in the critical FvK number for buckling. For bundles, this dramatic asymmetry can be attributed to the fact that "double-twist" generates the equivalent of positive curvature geometries (K eff > 0), but requires no strain variation along the bundle's length. In fact, it can be shown that the uniform double-twist pattern is the only texture that does not break lengthwise symmetry [40]. Therefore it is the only soft mode available for deformation (i.e. available to any γ = 0) because it doesn't lead to the immense strains caused by area dilation. Consequently, generating negative equivalent curvature (K eff < 0) through axisymmetric splay leads to an elastic cost for such modes that does not vanish in the k → 0 limit.
In Sec. III B we find that 7-fold defective bundles are still unstable to modes that (must) break lengthwise symmetry; however there exists a more exotic tilt pattern of lower energy that allows for buckling at a lower (though necessarily still nonzero) γ. As a consequence, the asymmetric splay mode that was described analytically in the current section is in fact preempted by this nonaxisymmetric shape mode, shown schematically in the inset to Fig. 2(b). Although this instability is triggered earlier (γ c ≈ 1, 500 and kR ≈ 3), we will show below that the negative-curvature tilt pattern breaks lengthwise symmetry as required, and therefore imposes a finite value for the critical fil-vK number. This behavior is in stark contrast to 5-fold defective bundles that are unstable for all γ.
As alluded to above, the linear stability analysis of axisymmetric modes does not necessarily capture the true symmetries of the most stable deformation pattern, nor the far-from threshold buckled configuration. Additionally, we have ignored the effects of bundle ends by taking the L/R → ∞ limit. In the following sections, we lift these constraints and compare our analytic results to those from simulations of a finite-length discrete-filament model of cohesive bundles.

III. DISCRETE MODEL OF COHESIVE FILAMENT BUNDLES
Here we introduce a bead-spring model of cohesive filaments. Our purpose is to determine the elastic energy  (17) and (18). The distance of closest contact (red line) represents the true separation between filaments i and j.
ground states of bundles containing disclinations, without constraining the symmetry of their deformed shapes, as was done in section II. This model treats individual filaments as semi-flexible and cohesive "featureless" tubes, that incur no elastic cost for lengthwise sliding of neighboring filaments, but do generate costs for lateral deformations that strain inter-filament distances.
A bundle contains N f filaments, indexed by i = 1 . . . N f , with each filament discretized into N v vertices, or "beads", indexed by n = 1 . . . N v . Vertex positions along a single filament are located at the position, x i,n , and we define i,n as the length of the line segment between vertices n and n+1 on filament i. The local tangent at n is defined asT i,n = (x i,n+1 − x i,n )/ i,n , from which the cost of intrafilament bending is defined In the limit that N v → ∞, this energy asymptotically approaches the standard elastic energy for a semi-flexible, worm-like chain. The elastic cost of cohesive interactions between neighboring filaments i and j are modeled as generic Hookean springs, where ∆ n,ij represents the distance of closest contact from vertex n on filament i, to a point along filament j. This distance intersects j at a right angle as shown in Fig. 3 [41]. In our discretized model, filaments are composed of line segments anchored to jointed vertices, where ∆ n, ij is calculated between a vertex and its neighboring segment, rather than between vertices [42]. For sufficiently large N v this model allows for frictionless sliding between neighbor filaments (particularly when they are straight). Assembling eqns (17) and (18), the total free energy of our discrete filament model is where the final sum is over all the nearest neighbor filaments j, to filament i. In Appendix A, it is shown that in the limit of N v → ∞ and N f → ∞, the elasticity of this model approaches the continuum limit described by eqns (2)and (5) where 0 is the initial intra-filament vertex spacing [43].
It is straightforward to including an intra-filament stretching cost that favors a constant i,n = 0 and acts to maintain inextensibility of filaments. This would be necessary, for example, for a physically accurate description of the interplay between 3D shape and the loss of cohesive contact at the ends of fixed length filaments. However, the present goal is to explore and analyze the buckling behavior of defective bundles, extrapolating to the L → ∞ limit where the effects of boundary interactions are presumably negligible. Hence, filament lengths are not fixed, and rather the vertical z coordinates of vertices are fixed at equally spaced layers (with a separation of 0 ). This has the effect of achieving smaller inter-filament elastic stretching at bundle ends than would occur in the fixedlength case. Strictly speaking, in this model the volume of the bundle and the length of the filaments are no longer conserved. However, this situation is arguably more relevant of certain biological or supramolecular fibers that self-assemble by adjusting length and radius simultaneously. Regardless, we discuss how the non-generic treatment of the bundle ends influences the buckling behavior for finite bundle length.

A. 5-Fold Disclinations
Here, we consider the shape transition of bundles with a central 5-fold disclination, introduced through a fixed topology of inter-filament elastic bonds. We focus on the case of a high-aspect ratio with L/R = 8, as we find that finite-length end effects play a relatively small role in the their buckling behavior. We consider bundles with with N f = 306 filaments (radius R ≈ 10∆ 0 ) and a vertex spacing of 0 = 0.2∆ 0 . Energy minimization was performed using the GSL conjugate gradient package for C [44]. Below we explore the structure and energetics of stable and meta-stable states beginning with highly flexible filaments with γ = 25, 000. Then, after energy minimization, γ is reduced by increasing B, and the energy is minimized again. This process is repeated over 50 steps (in even logarithmic decrements of γ) until γ = 0.25.
Results for 5-fold disclinations are shown in Fig. 4, where we plot the mean filament twist angle, θ , defined as the mean angle of all filaments with the centerline (i.e. cos θ i,n =T i,n ·ẑ) vs. γ. Data point colors represent different initial configurations generated by applying an azimuthal displacement pattern to all filament vertices of the form 0.2r cos(nπz/L), where n is an integer that counts the number of times the handedness of the helical twist changes along the bundle length. This procedure allows us to bias the lengthwise symmetry of distinct equilibrium shapes. The case of n = 0 produces a homogeneously twisted bundle, where all filaments possess an identical pitch, which is the lowest energy state for all tested values of γ. Note that bundles with 5-fold disclinations are always unstable to homogeneous twist (n = 0) for all values of γ, consistent with the L → ∞ linear-stability results described above.
In addition to uniform pitch states, we see in Fig. 4 that metastable oscillating twist states are mechanically stable for sufficiently large values of γ. These torsionally wrinkled shapes are characterized by an alternating direction of twist along the z axis. Our discrete-filament simulations find that for a given n-wrinkled shape, there is a value of fil-vK, designated γ * (n), below which, the bundle becomes unstable to a lower-n structure. These points are highlighted as the large dots in Fig. 4(a). As γ is decreased, eventually γ * (n) is reached, and the bundle becomes unstable and undergoes a large transition to a new lower energy and lower n of torsional wrinkles. To highlight this trend, Fig. 5(a) shows the total energy vs. γ. For large γ, the transitions between alternating handedness of twist become sharp kink-like boundaries, consistent with bending being concentrated over length scales proportional to λ b = R/γ −1/2 .
The region of stability of n-wrinkled bundles is consistent with the linear-stability analysis of the continuum model, by assuming that torsional oscillations are commensurate with the finite length of the bundle, or k n = πn/L. On one hand, γ c (k), predicted by the continuum model, defines the point at which the straight bundle becomes unstable to torsional wrinkling at wave vector k, while γ * (n), measured from simulation results, is the smallest value that an n-wrinkled bundle is observed to be stable. These two thresholds always satisfy γ * (n) > γ c (k n ). It is not clear what limits the ability to resolve mechanical equilibrium of the n-wrinkled state all the way down to the parallel state (i.e. θ = 0); presumably this derives from the combination of the inherent precision limit of our discrete-filament model and the vanishing of energetic barriers between nearly unstable and stable modes (with lower n). Notwithstanding the loss of stability as γ approaches the limit of stability for the n-wrinkled mode, we estimate the value of this threshold by fitting the γ-dependence of an n-wrinkled mode to the region near γ γ * (n) where θ 0 is the maximum twist angle far from the transition point (22.1 • for a 5-fold disclination), and ζ is a value that regulates the speed of the transition. This particular form has two motivations: first, the expectation that near to the stability threshold θ ∼ |γ − γ c | −1/2 , a characteristic of a supercritical bifurcation; and second, the predicted γ-dependence for equilibrium uniform n = 0 twist of 5-fold defective bundles is expected to have the form of eqn (20), with θ 0 = 2 √ 3/9 radians, ζ = 32/3 and γ c = 0 [30]. This fit is shown to agree well with the n = 0 results in Fig. 4 (dashed line). The values of the γ c = 0 for n ≥ 1 extracted from fits to eqn (20) are shown in  While the total energy of a wrinkled bundles is always found to be increasing with n for a given γ, we find evidence that, in the large-γ limit, oscillating twist structures for high-n tend towards a lower elastic energy than lower-n structures. Fig. 5(b) shows the elastic contribution to the total energy for large values of the fil-vK number, which shows that E elas from relatively high-n states (e.g. n = 6) tends to decrease faster with γ than lower-n structures. Extrapolating this to even larger values of γ suggests that for sufficiently flexible filaments, highly-wrinkled bundles (n → ∞) could become the lowest-energy shape equilibria (even lower than the uniform twist state) in the γ → ∞ limit, where the cost of filament bending is negligible in comparison to the elastic cost.
Recalling the analogous case of crystalline membranes, Seung and Nelson argued that in the asymptoticly flexible limit of γ s → ∞, the far-field elastic stress of a 5-fold disclination can be completely screened. This stress focusing is achieved by a nearly isometric conical shape that concentrates Gaussian curvature to the disclination position [6]. For the torsionally wrinkled bundles, the tilt pattern that concentrates K eff to the bundle center is not the uniformly twisted one, but one where filaments tilt rapidly from theẑ direction within the core of the bundle, and adopt a constant θ in the outer bulk of the cross section. Fig. 5(c) shows the twist angle θ as a function of a filament's radial distance from the centerline. While this tilt pattern is possible within a given cross section of the bundle, it requires deviations from the constant helical pitch and introduces shear deformations that grow along the bundle length. Nevertheless we observe that for large γ, bundle shapes tend towards a similar "curvature focusing" geometry, made possible the torsional wrinkles. Fig. 4(d), shows triangulated surfaces with inter-vertex distances equal to the inter-filament distances in the discrete filament model mid-way between the alternating wrinkles. These surfaces show the progressive focusing of curvature towards the central defect as the number of wrinkles increases with γ, similar to crystalline membranes (see Appendix C for how to calculate the Gaussian curvature). If true, such a feature would be important for stabilizing nontrivial topologies in crystalline columnar materials.

B. 7-fold disclinations
In this section, we analyze the shape-transitions of bundles with centered 7-fold disclinations. Unlike the torsional wrinkles of 5-fold defects, our discrete-filament simulations show that 7-fold defects lead to buckled shapes that significantly break axisymmetry, in a manner unlike the splay-undulation ansatz analyzed in Section II B. Nevertheless, despite the different optimal tilt pattern, discrete-filament simulations do reveal that 7-fold defects favor shapes that break the lengthwise symmetry of inter-filament strains, and, as required, only above a critical γ. However, the tilt pattern is such that area dilation, and therefore costly in-plane strains, are minimized.
Given the complexity of the optimal buckled shapes, we first focus on the limit of infinitely rigid filaments (γ → 0) with a finite length, L. Bundles of N f = 106 filaments (R ≈ 5∆ 0 ), of different lengths are shown in Fig. 7(a-b), with L = R and L = 2R respectively, showing a deformation pattern with two generic features. Most obvious is the tilt of the bundle's centerline with respect to the z axis. Superimposed on this near-uniform tilt is a more subtle pattern of tilt-variation within the bundle's cross section. This pattern is more easily illustrated via the projections of filament tilt in a plane perpendicular to the bundle centerline, as shown in Fig. 7(cd). The "double vortex" pattern viewed in this perspective reveals the surprising emergence of twist driven by 7-fold defects, far from the radial splay pattern assumed on the grounds of axisymmetry. This pattern, which we call the counter twist tilt pattern, is composed of two double-twisting domains of opposite handedness.
To show that this tilt pattern effectively screens the defect-induced stresses, we analyze the distribution of equivalent Gaussian curvature, K eff , using the discretefilament analysis from Appendix C. Fig. 7(c-d) shows that regions of K eff < 0 are predominantly focused at the central 7-fold disclination. One can understand how this counter twist tilt pattern is the metric equivalent The colored lines represent the distance of closest contact between neighboring filaments. The blue circumferential and red radial paths are relatively unaffected by tilt, while the yellow radial path is shortened due to filaments tilting. (f) Filaments mapped to an equivalent surface with identical colored paths, but with distances now shortened by curvature rather than filament tilt.
to negative Gaussian curvature through the schematic Figs. 7(e)and (f). Here we highlight distinct interfilament paths (colored) in the cross section. In this metric analogy, we are interested in the accumulated distance of closest approach between neighbors spanned along the path, rather than the path length in the planar crosssectional cut. This distance is dependent on, and can only be shortened by, filaments tilting into the direction of their separation. In Fig. 7(e), the circumferential (blue) path encompasses the 7-fold disclination, and is relatively unperturbed by filament tilt. However, the radial (yellow) path along the mid-line separating the distinct double-twist domains, is effectively shorted by tilting. Hence, this tilt pattern allows the bundle to shorten distances along (certain) radial directions while keeping the circumferential distance the same relative to a planar geometry, a hallmark of negative curvature geometries. An equivalent negatively curved surface, shown in Fig. 7(f), shows the same filament tangents and distances of closest approach between them, but now lengths are adjusted by out-of-plane surface curvature rather than in-plane filament tilt. Both patterns are effective at relaxing the compressive hoop stresses generated by 7-fold defects, by expanding the circumference while keep radially separated filaments close to their preferred separation distance.
For this rigid filament case, we note that the tilt pattern is highly sensitive to length. For example, in Fig. 7(c-d) the tilt variation and K eff decreases significantly from the L = R case to L = 2R. To understand the length dependence, we analyze this pattern using two quantities: the first measures the tilt of the centerline, θ , or the mean angle of filaments with respect to the z axis; and the second measures the prominence of the double vortex tilt pattern seen in Fig. 7(c-d), δθ , defined as the mean tilt angle of filaments away from the centerline. From these parameters we can estimate how the 7-fold bundle energy varies with, θ and δθ , as well as L and R. (For a derivation of the comparison between our continuum elastic and discrete bead-spring models, see Appendix A.) We begin with a simple ansatz for counter twist which has two double-twist patterns centered on x = 0 and y = ±R/2, and a mean orientation of θ x. Assuming that θ δθ, we may use eqn (1) (and Appendix C) to computed the effective negative curvature of the pattern, giving K eff ≈ −3 θ δθ/R 2 . Hence, the negative curvature geometry relaxes the elastic energy over the bulk of the bundle by an amount δE relax ≈ −Y V |s| θ δθ, where Y ≈ −1 0 . However, generating this tilt pattern requires two additional costs in rigid filament bundles. First, the tilt variation leads to in-plane strains (dilation and shear) that grow along the bundle length as u ij ≈ δθz/R. This cost leads to an elastic penalty that grows rapidly with length δE elas ≈ Y V (δθ) 2 (L/R) 2 .
FIG. 8. Plots of (a) the mean tilt angle, and (b) the mean tilt angle deviation, vs. bundle aspect ratio L/R. Dotted lines show the power laws predicted from eqn (22), which follows the limit R/∆0 1. Inset plots show the same data plotted vs. the dimensionless combinations given in eqn (22). Note that the double-asymptotic limit of L/R R/∆0 1 underlies eqn (22), highlighted by the yellow (lighter) data points in these insets.
Finally, the mean tilt of the bundle axis introduces a stretching cost at the ends of the bundle through the tangential "slip" of neighbor filaments. For small tilts, the length of the "slipping" regions are slip ≈ ∆ 0 θ [45], over which the inter filament distance is stretched by an amount ∆ − ∆ 0 ≈ ∆ 0 θ 2 , leading to an elastic cost for tilt at the bundle ends of δE ends ≈ Y R 2 ∆ 0 θ 5 .
Combining these energetic terms and minimizing with respect to mean-tilt and tilt-variation, we find While this scaling suggests that both angles vanish as L → ∞, δθ (i.e. the double-vortex tilt pattern) is expected to decrease far more rapidly with bundle length. Fig. 8(b) indeed shows this decrease with L/R, and is in further agreement with the weaker power-law dependence on R/∆ 0 . Fig. 8(a) shows that the mean tilt value, θ , falls with both L/R and R/∆ 0 , in numerical agreement with the power law of eqn. (22) in the asymptotic regime L/R R/∆ 0 1. This scaling suggest that such tilt patterns for shorter bundles of rigid filaments better screen defect stresses than longer ones.
We turn now to the case of finite flexibility (i.e. γ = 0), and map out the γ-dependence of the 7-fold disclination buckling transition with the discrete filament bundle model. Results are found in a manner similar to the case of 5-fold bundles, but with γ increased in a stepwise manner between energy minimizations in order to investigate the buckling transition. Final results are shown in Fig. 9(a) for a bundle with N f = 428 (R ≈ 10∆ 0 ) and three different aspect ratios L/R, where we measure the degree of shape buckling by δθ , the mean variation of filament tangents with respect to centerline.
For short bundles, we find a gradual increase of δθ with γ, consistent with the intuitive notion that filament flexibility simply reduces the cost of the counter twist tilt pattern previously shown to be favorable for γ → 0. As bundle length increases to the large aspect ratio limit, L/R 1, we find that this gradual increase sharpens. For small γ, consistent with the scaling above, δθ tends to zero with increasing L/R. However, we find that the large-γ buckling persists and tends toward an L-independent behavior as L/R increases. These trends are consistent with the emergence of a finite buckling threshold as L → ∞, as predicted by the linear stability of the axisymmetric model in Section II B. The behavior of the critical fil-vK number for this emergent doublevortex instability is schematically illustrated in the inset of Fig. 2(b). Compared to the rotationally symmetric constrained ansatz from the continuum elasticity stability analysis, this tilt pattern is extrapolated at the L → ∞ limit to occur at a lower γ c ≈ 1, 500 and kR ≈ 3. As discussed, this finite-γ threshold for buckling by 7-fold defects can be attribute to the elastic costs of breaking lengthwise symmetry of inter-filament spacing, a necessary consequence of the K eff < 0 tilt pattern. Fig. 9(a) shows the tendency with increased L/R towards something like a second-order transition occurring around γ ≈ 2, 000. This value is notably an order of magnitude smaller than the threshold predicted by the continuum theory (from Fig. 2(b)). In large part, we expect this discrepancy in the buckling threshold is due to the fact that 7-fold bundles adopt shapes that are far from the axisymmetric ansatz of radial splay, and presumably, relax elastic costs far more efficiently. The specific pattern of buckling is comparable to the counter twist tilt pattern exhibited by finite-L rigid filament bundles shown in Fig. 7. Additionally, we find that for large L, complex patterns are (at least for intermediate γ) consistent with a varying pattern of local counter twist along the bundle's length, where the span of one counter twist zone is much smaller that L. We can expect that the length of these counter twist zones, h, is set by a balance between inter-and intra-filament elasticity. We then balance over a vertical length h the inter-filament elastic cost, ∼ Y (δθ) 2 (h/R) 2 , with the bending cost, ∼ K(δθ/h) 2 , to find an optimal size at This length is notably the same characteristic scale of optimal splay undulations. Beyond this length scale, the buildup of in-plane strains must be relaxed by a boundary layer that allows filament directions to reorient (e.g. through a reversal of the counter twist tilt pattern).
The full 3D structures for various values of γ can be seen in Fig. 9(b); while Fig. 9(c) shows maps of cross-sectional mean Gaussian curvature K eff at distinct heights. For intermediate γ, the midsection of the bundle appears to alternate in the handedness of the double vortex, seen in Fig. 9(c). Although the exact tilt pattern is not as simple as the ideal case found in Fig. 7(c), it clearly is not the splay ansatz proposed in Fig.1(d), and can be interpreted as a non-symmetric double vortex pattern. However, as γ is increased still further, the buckled structures lose their alternating double vortex character and instead become more disorganized, possibly an indication that larger effects of inter-filament "friction" at high curvatures in our model prevent full equilibration. Nevertheless, we note that these "crumpled" shapes continue to show increasing negative effective curvature for larger γ, as well as a further decreasing distance between adjacent "crumples", consistent with the expected decrease in the span of counter twist domains h ∼ γ −1/4 . All in all, these shapes are far from the proposed radial undulations in Section II B, a detail that will be discussed in the next section.

IV. DISCUSSION
In the classical theories of disclinations, dating back to the geometric constructions of Volterra [46] up to modern elastic theory treatments [47,48], positive vs. negative disclinations are symmetric with respect to their energetic costs in linear elastic theory. This asymmetry is only weakly broken in the 2D crystals, due to the unequal costs of conical-vs. saddle-like bending, and leads to few percent shift of the critical FvK number for buckling [6]. Alternatively, disclination driven buckling in columnar structures appears to be in a distinct class, where the geometric packing constraints select one of the two signs of defects as particularly low energy, as shown schematically in the single defect self-energy plotted in Fig. 10. Five-fold defects drive torsional buckling of bundles of any diameter and filament stiffness (i.e. any nonzero γ), while for L → ∞, 7-fold defects only drive shape transitions above a threshold bundle diameter (threshold γ). Simply put, this implies that positive disclinations should always be more abundant than negative disclinations in columnar/filamentous bundles. Here, we discuss the geometric origins of this asymmetry between defect signs, as well as the susceptibility to defect-induced instabilities for various material systems.
Simulated elastic energy ground states for 7-fold defects reveal a surprising symmetry breaking feature, where buckling breaks both the lengthwise and axial symmetry of the initially parallel bundle (e.g. Fig. 7). What is perhaps most surprising, is the emergence of the common tilt pattern of local double-twist, in response to both positive-and negative-charged disclinations. Doubletwist is a cylindrically symmetric pattern of filament tilt, in which filament orientations twist around a central axis. A single double-twist domain has the tilt pattern t ⊥ Ωrφ and a metric equivalent curvature K eff = 3Ω 2 , which is consistent with the ability of this structure to screen positive disclination stresses. Why then, should 7-fold defects also drive the formation of double-twisted tilt patterns if they desire negative equivalent curvature?
The answer to this question derives from the interplay between the geometry of 2D tilt patterns in a given cross section and the lengthwise variation of inter-filament distances in the bundle. To illustrate this, we define the strain variation, v ij ≡ ∂ z u ij , to measure changes in intercolumn strain along the bundle. From eqn (3), we have where ∂ z t κn, with κ and n being the respective local curvature and normal of a filament's Frenet-Serret frame. This tensor shows that certain tilt patterns in a 2D cross section require the buildup of elastic strains up-or downstream from the section when v ij = 0. Furthermore, it can be shown that uniform double-twist is, in fact, the only non-parallel tilt pattern which does not required lengthwise variations of inter-filament distances [40]. Because double-twist has strictly anti-symmetric in-plane gradients (i.e. ∂ i t j = ij Ω), it is easily verified that all components of v ij vanish for a single, uniform-Ω domain. Again, the lack of lengthwise strain variation for uniform double-twist is linked to the absence of a threshold fil-vK for torsional buckling in 5-fold defective bundles. For the case of 7-fold defective bundles, it is now easy to see why the axisymmetric radially splay ansatz is not preferred. Although it does generate a favorable local distribution of K eff < 0, this tilt pattern unfavorably requires v ii ∇ ⊥ · t ⊥ = 0 everywhere in the cross section. Therefore, this axisymmetric geometry is presumably a poor compromise between the preference for negative K eff and minimal elastic strain variation along the bundle. Alternatively, the counter-twist observed in Section III B, composed of two oppositely rotating double-twisted domains, focuses a region of highly-negative curvature between the opposing domains (near to the disclination position). While simultaneously, a minimal lengthwise strain build-up along the bundle cross section is maintained. This is illustrated in Fig. 11, which shows the map of double-twist and splay-density (v ii ), for a 2D cross section of a simulated 7-fold defective bundle of rigid filaments found in Section III B. Unlike radial splay, the counter-twist tilt pattern effectively expels splay from most of the bundle's cross section, while generating a sufficient measure of K eff < 0 (Fig. 7(c)) near the bundle's center to screen the in-plane stresses generated by the 7-fold defect.
The above arguments suggest that the unique geometry of double-twist make this tilt pattern a potentially inexpensive "building block" for more complex 3D buckled cohesive bundles, well beyond the seemingly ideal centered 5-fold defects. The generic emergence of doubletwist can further be illustrated by the shape equilibria of bundles possessing off-centered 5-fold disclinations, shown in Fig. 12. Here we see that as the defect position, r d , increases away from the bundle center, the centerline of the bundle buckles helically with a curvature that increases with γ. The equilibrium curvature of the centerline, κ c , depends non-monotonically on r d . Superficially, this buckling pattern would seem to imply a qualitatively different responses of non-axisymmetry defect distributions. On the contrary, this "writhing" bundle equilibria are in fact themselves regions of a uniform double-twist domain, but with the center of rotation located away from the bundle center and eventually outside of the bundle's cross section. This behavior is likely the result of non-linear (large-tilt) corrections to the metric geometry of uniform double-twist described previously [31]. Here K eff is concentrated at at the center of rotation and tapers to an effectively flat geometry for large radial distances from this axis (compared to helical pitch). Thus, the shifting center of double-twist rotation is likely driven by the polarized stress distribution created by off-center defects. The bundles ultimately straighten for larger r d due to the elastic screening of defect stress by the free boundary [30,36].
Next, we consider the relevance of defect-induced shape buckling in cohesive filament bundles realized for a range of physical systems. In particular, we estimate the range of accessible fil-vK numbers based on measured or predicted values of their resistance to bending, their elasticity of inter-filament cohesion, and their observed values of bundle radius, R. The buckling behavior to defects is critically sensitive to γ ≡ (R/λ b ) 2 , which is the ratio of bundle size relative to the material dependent length scale, λ b = K/Y , itself a measure of the ratio of intrafilament to inter-filament elasticity. The bend elastic cost   [59,60] to shorter-range and stronger binding due to polyvalent counterion binding [58,61]. Here, cohesion between polymer microfiber arrays is modeled by capillary bridging [62,63].
can be estimated simply from the bending stiffness of filaments B, and their diameters d, as K ≈ B/d 2 . Far less well characterized is the elasticity of inter filament cohesion, which in turn, controls the elasticity of the 2D filament array. Following refs [64,65], we estimate Y by assuming that the elastic stiffness of inter-filament cohesion (strictly speaking the curvature of the inter-filament binding potential) can be estimate as /σ 2 , where is the cohesive energy per unit length between neighbor filaments, and σ is a length scale characterizing the cohesive range. Neglecting numerical prefactors accounting for lattice geometry, this estimate then gives Y ≈ /σ 2 .
In Table I, we compare values of γ for vastly different classes of cohesive bundles: carbon nanotubes; biological filaments (DNA and microtubules); and polymeric microfibers assembled by surface (or capillary) interactions. Significantly, we find that the materials can reach a γ that spans up to 15 orders of magnitude, up to γ ∼ 10 6 . In general, this suggests that almost all cohesive filament bundle assemblies are susceptible to buckling by 5-fold disclinations, but in general, most conditions fall below the range where they are perturbed by 7-fold disclinations (i.e. γ 2000). For example, nanotube bundles that condensed from solutions [51] grow to relative narrow widths (∼ 100 nm) and reach only γ ∼ 10 −6 − 10 −5 .
Alternatively, much larger (∼ 5 µm) "yarns" generated by spinning from nanotube "forests" extend far into the large fil-vK range, potentially exceeding the critical γ for the 7-fold disclination-induced shape buckling shown in Fig. 9.
Microtubule bundles formed by polyvalent counterion condensation have been observed to exhibit longwavelength undulations, interpreted as 3D writhing configurations [58], not unlike the helical buckling of offcenter defect patterns of Fig. 12. For these systems there is also some evidence of "irregular" packing in the 2D cross section, though at present there has been no attempt to quantify defect distributions in these experiments. It has been proposed that a non-equilibrium process of cohesive condensation of filament bundles may often lead to the trapping of topological defects into the cross-sectional order [66], which may have consequences on their structure and assembly. From the values of interaction and mechanics of microtubules, in particular for relatively brittle binding by polyvalent counterion bridging, we modestly estimate that γ ∼ 0.1. This small allowance for flexibility may drive a detectable degree of helical buckling in the presence of positive disclinations, as suggested by the low-γ range of Fig. 5.
Finally, we conclude with a simple estimate for macro-scopic and elastic filaments, held together by cohesive surface contact (i.e. the range of attractions is small compared to the filament diameter). In this case, both filament bending stiffness, B ≈ E f d 4 , and inter-filament contact stiffness, Y ≈ E f , are proportional to the elastic modulus, E f , of the filaments. Hence, their ratio is expected to be independent of E f , and therefore roughly in this strong cohesive contact regime, we expect that the fil-vK becomes independent of material parameters and simply dependent on the number of filaments within the cross section, γ ≈ N f .

V. CONCLUSION
We have determined the disclination-induced buckling instabilities for cohesive filament bundles, and resolved their dependence on the size and mechanics of the assemblies. These results point to the ratio of intra-filament to inter-filament elastic stiffness as the key materialdependent quantity that regulates the emergent shapes of defective bundles. Compared to flexible membranes, the equilibrium shapes of bundles exhibit a far greater nontrivial dependence on the types of defects. This results from the interplay between the local metric geometry of 2D bundle cross sections and the lengthwise variations of filament spacing within the bundle. We find that bundles are vastly more susceptible to shape deformation by 5-fold defects (positive disclinations) than 7-folds. While bundles possessing other cross-sectional symmetries (e.g. 4-fold) may require anisotropic elastic costs as well as distinct topological charges of defects, we expect this basic conclusion to hold independent of lattice symmetry, as the geometric principles underlying metric coupling to tilt and uniform longitudinal spacing are generic.
The principles developed in this study suggest new strategies for engineering the 3D shapes of bundles through their controlled 2D packing. We envision various synthetic strategies that may be exploited to template columnar or filamentous assemblies with controlled topological defects in their cross-sectional order. This feat would be much in the same spirit of kirigami engineering of point-like disclinations used to engineer the 3D (and self-folding) shapes of sheets. This could be achieved, for example, by i) the capillary cohesion of nanofabricated 2D arrays of high aspect ratio flexible pillars [62], ii) the controlled defect formation in columnar assemblies grown epitaxially from a lithographically templated surfaces [67,68], and iii) so-called DNA-origami techniques to engineer dsDNA bundles [69] with programmed 2D cross sections, containing intentionally placed disclinations and dislocations. For example, as shown in Fig. 12, controlling the placement of a 5-fold defect and lead to spontaneously writhing and twisting 3D configurations, dramatically reshaping the responses of the assembly to a range of stimuli (e.g. mechanical, electronic, photonic, etc.). In systems where the cohesion between filaments and their elastic stiffness may be externally tunable (say in temperature-or field-responsive materials), we envision that the 3D shapes of such bundles could be adjustable, driven to wind and unwind responsively.
To a first approximation it can be anticipated that multiple elementary disclinations of the same charge lead to effects equivalent to single higher charged defects, namely, added drive for positive or negative curvature bundle textures. However, given the highly non-linear dependence of buckling to disclination type, the response to even a single edge dislocation, i.e. a 5-7 "dipole", is likely to be far more complex for bundles than for their geometrical analogs of 2D membranes [6]. Previous work provides some insight into how dislocations drive certain tilt patterns [70]. First off, dislocations have an extra degree of motion over disclinations, namely their orientation. In much the same way 5-fold disclinations are expected to be more prevalent in filament bundles, dislocations are expected to be oriented such that the 5-fold end is turned towards the center of twist (corresponding to the removal of a partial row extending radially to the boundary), this generate stresses that can be relaxed by positive curvature tilt patterns (e.g. twist). Furthermore, the strength of the tilt-dislocation stress coupling is strongly position dependent, with a preference for radially-aligned dislocations situated at R/ √ 3 from the bundle center. Alternatively, misoriented dislocations (e.g. with the negative disclination end closer to the center) cannot be relaxed by twist, and may instead lead to non-axisymmetric and longitudinally asymmetric deformations modes. Lastly, unlike disclinations, whose strength is quantized in multiples of π/3, the net strength of dislocations to drive shape instabilities depends on the ratio of Burger's vector (a microscopic dimension proportional to lattice spacing) to bundle size, b/R. Therefore, we would anticipate that the drive to buckle the bundle into a 3D shape may be more sensitively controlled by the collective effects (number, orientation, and locations) of multiple dislocations arranged in a bundle as compared to one or a few isolated dislocations.
While such experimental directions remain to be explored, the principles governing the response to elementary 5-and 7-fold disclinations lay the ground work for engineering custom-made defect configurations in 2D filament arrays, which could be used to "program" specified 3D shapes of bundles. However, beyond the isolated disclinations discussed here, the geometric and mechanical principles that govern the collective responses to multi-defect patterns remain open.
NSF through Award No. PHY 16-07611) where some of this work was completed. Numerical computations were performed on the UMass Shared Cluster at the Massachusetts Green High Performance Computing Center.

Appendix A: Continuum and discrete model parameters
Here we describe the correspondence between parameters of the discrete filament model of Sec. III and the elastic moduli of the continuum theory of columnar bundles presented in Sec. II. We begin with the squared curvature of the ith filament with tangent T i (s) where lim n →0 is the continuum limit. We can convert from the discrete model of bending energy in eqn (17), to the continuum model in eqn (5) by summing over all "bending bonds" on filament i, and then summing over all filaments, yielding the total bending energy where dA N f i=1 ρ −1 0 , and ρ 0 = ∆ −2 0 / √ 3 is the stressfree areal density of filaments giving, For the discrete model of the inter-filament elasticity, from eqn (18) we have the elastic energy of the nth layer, where ∆ 0 is the preferred local spacing between filaments. Structurally, the bundle model is essentially N v stacks of the hexagonal bead-spring model of ref. [6], with a geometrically non-linear coupling between local tilt and lattice strain. Thus, summing over these layers we have from which we have the elastic constants, corresponding to 2D Youngs modulus, and 2D Poisson ratio ν = 1/3. From these the fil-vK number may be estimated from the discrete model parameters as The boundary conditions on the sides of the bundle are simply σ 1 rr = σ 1 rθ = 0, or specifically λρ(R)/R = (λ + 2µ)∂ r ρ(R) = 0 (B9) ∂ r τ (R) − τ (R)/R = 0. (B10) And finally, we have the boundary conditions for the derivatives of the displacements at the ends of the bundles, but we will neglect these by assuming that solutions are periodic, and of the form ρ(x) = δu r (r) cos(kz); τ (x) = δu φ cos(kz).
To compare to a finite length bundle, we might consider wavelengths that are commensurate with the bundle length, k = 2πn/L, though to be clear, these purely sinusoidal deformations will not allow us to match the free end boundary conditions from eqns (B3) and (B4). Presumably, a boundary layer is required to match the purely periodic solutions we consider, to the free end calculations. Therefore, we work under the assumption that the length scale of this boundary layer will vanish proportional to K/Y , and hence can be ignored for large aspect ratio bundles (L/R 1), and large bundle fil-vK number (γ 1). To proceed, we rewrite the equations in dimensionless variables, by measuring all lengths in units of bundle width R, and stresses in units of Y . Doing this, and recalling the definition of the 2D Poisson ratio, ν = λ/(λ + 2µ), we rewrite eqns (B7) and (B8) in the form of eqn (12) with the effective "potentials" In this way, we have recast the linear stability calculation in terms of an eigenvalue problem, with the boundary conditions νδu r (R) + ∂ r δu r (R) = 0 (B16) ∂ r δu φ (1) − δu φ (1)/R = 0. (B17) For the linear stability calculation, we are interested in the ground state solution, i.e. the smallest values of r or φ that are consistent with our boundary conditions. This will correspond to the first instability-the lowest value of γ-at which a given periodic mode k becomes unstable. In the main text and Fig. 2, we solve these equations for the most unstable wavenumber, k, given centered 5-fold (s = +2π/6) and 7-fold (s = −2π/6) disclinations.
Appendix C: Discrete approximation of equivalent Gaussian curvature For an arbitrary discrete surface composed of vertices, edges, and triangular faces, the Gaussian curvature can be defined in terms of the deficit of internal angle at a single vertex, as where ψ α is the internal vertex angle for face α, A is the summed area of all the faces attributed to the vertex, and f is the number of triangular faces adjoining the vertex [71,72]. We can generate an equivalent mesh for a bundle by finding the points by intersection of all filaments with a plane, here chosen to be the z plane; followed by a Delaunay triangulation. While this mesh is itself flat, with zero Gaussian curvature, the relevant distance between filaments is the distance of closest contact which is generally out of plane (i.e. similar to ∆ n,ij , but defined within a cross-sectional plane rather than at vertex n).
Explicitly, if the in-plane spacing between two filaments, i and j, is ∆ = x j − x i , then the out of plane distance of closest contact from filament i to filament j, is defined as Therefore, by computing the corrected side length, |∆ ⊥ | of a given triangle, we calculate the internal angles, ψ α . From this we compute the effective Gaussian curvature, K eff , at given filament (and at a given z) by summing the internal angles in formula eqn (C1) for the local triangular neighbor array. In essence, this method works by distorting the dimensions of every triangle based on the relevant distances of closest contact, then stitching them back together in a manner that preserves the original network topology, but requiring out-of-plane orientations. The total effective curvature of the surface, K total eff , is simply the sum of the equivalent Gaussian curvature of each individual filament at that z plane.
Although it is possible to always calculate the equivalent Gaussian curvature for a given pattern of filament tilt, it is not always possible to find an example of the equivalent surface that is embeddable in R 3 , or one that is strictly unique, possessing the same metric data [73]. The reconstructed surfaces in Fig. 5(d) were found by energetically relaxing the positions of a bead-spring model of vertices, with a spring constant k s and a spring rest length of ∆ ⊥ (which is unique for each vertex). Energy minimization was performed on an initially axisymmetric conical surface, and proceeded until the total energy of all the springs fell below 0.001k s .